NN 実装のための 行列演算・画像演算 活用　全結合Affin処理の例

cuBLAS と NPP を使用して、GPUで動くＮＮの自前実装を行ってます。用いた関数群を紹介します。

#背景
* 深層学習の理解を深めようと自前で実装を行ってます。
* GPU速い。[15K￥位のgeforce GTX 1050tiでも、CPUに比べ10～100倍以上の実行速度があることを体感](http://qiita.com/t-tkd3a/items/fbd1d97fae51a880ad1f)しました。
* GPUで演算を行わせるため、CUDA Toolkit の行列演算(cuBLAS)と 画像演算(NPP)のライブラリを活用します。

本記事では、行列演算と画像演算の活用しての、NN の全結合のAffin処理 を 公開します。

#参考文献
* [cuBLAS CUDA ToolKit](http://docs.nvidia.com/cuda/cublas/#axzz4hCMGJZPb) 
* [About BLAS - BLASの簡単な使い方](http://azalea.s35.xrea.com/blas/about.html)
* [NPP CUDA Toolkit](http://docs.nvidia.com/cuda/npp/#axzz4hCMGJZPb)
* [IPP インテル(R) インテグレート パフォーマンス・プリミティブ リファレンス・マニュアル(PDF)](http://www.intel.co.jp/content/dam/www/public/ijkk/jp/ja/documents/developer/IPP_2_Final_i.pdf)

#cuBLAS

##行列の積演算 図示

![cr-01-1.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/3ee88c01-afd3-c9e4-7ef2-638bff061935.png)
↓の図は↑の式を示します。　
![cr-01-2.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/e27e9313-e2c4-1b56-8113-8bd287bd7f84.png)
この図示方法、行数 列数の関係が一目で分かるのでいいですね。

## 行毎 および 列毎 の総和
[gemv 関数](http://docs.nvidia.com/cuda/cublas/#cublas-lt-t-gt-gemv) を用いて 行毎 または 列毎の 総和を算出する方法
gemv 関数は 行列と ベクトルの 積 を算出する関数
　$y = αAx + βy$
　　$α,β$ はスカラ値、$x,y$はベクトル値、$A$は行列
　　引数で 行列$A$の 転置 の指示が可能

![cr-02_A.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/73651097-cc0a-17b4-c382-46951d66384e.png)

↑の行列を例に、gemv関数で 行毎 または 列毎 の総和の算出方法を示します。

#### 行毎の総和

![cr-02_A行.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/3ee11e56-0c66-644b-bd17-a82da7317fa0.png)

* ベクトル$x$ の値は 全て1.0 
* ベクトル$y$ は、 行列$A$ の 行毎の総和 となる。

#### 列毎の総和
![cr-02_A列.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/cabedf6b-48b6-5c12-7e3c-0793a5375b7a.png)

* 行列$A$ を転置。
* ベクトル$x$ の値は 全て1.0 
* ベクトル$y$ は 行列$A$ の 列毎の総和 となる。

補足
> * 入力画像のPixel が 整数型 の場合には、 NPP の SumWindowRow関数 , SumWindowColumn関数 が利用可能。
*  入力画像が float型 の場合に適した関数は見つからず、gemv関数で行ってますが、 適した関数があるのでしょうか？


## Broadcast
[ger関数](http://docs.nvidia.com/cuda/cublas/#cublas-lt-t-gt-ger) を用いて 行 または 列 のBroadcast を行う方法
ger関数は 列ベクトル と 行ベクトル の積の 行列 を算出する関数。
　$A=αxy^T+A$
　　$α$ はスカラ値、$x,y$ はベクトル、$A$は行列。


### 列ベクトル の Broadcast

![cr-05.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/908c7cd9-024e-cdb0-28bb-c405618146c4.png)

* 列ベクトル$x$ が コピー元。
* 行ベクトル$y$ は 全て1.0
* 結果行列 は $x$ をBroadcast したものとなる。

### 行ベクトル の Broadcast

![cr-04.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/c05d8256-134c-02bf-4748-c34089657396.png)

* 行ベクトル$y$ が コピー元。
* 列ベクトル$x$ は 全て1.0
* 結果行列 は $y$ をBroadcast したものとなる。

補足
> もっと適した関数があるのでしょうか？


# NPP 
画像演算用のライブラリです。
cuBLASと補い合って使え便利です。

抜粋したものを下記します。
・ 下表内の"画像"は、画像内の各Pixel値 を示す。
・ 画像どうしの演算は、Width Height Channel bit数 が等しい前提で。
・ 各Pixel値 を float で扱える関数を抽出。

## 算術演算
|関数識別名|概要|備考|
|:-------|:-----|:-----|
|Add|画像 ＋ 画像|
|AddC|画像 ＋ スカラ値|
|Sub|画像 － 画像|
|SubC|画像 － スカラ値|
|Mul|画像 × 画像|アダマール積|
|MulC|画像 × スカラ値|
|Div|画像 ÷ 画像|
|DivC|画像 ÷ スカラ値|
|Abs|画像の絶対値|
|Sqr|画像$^2$|
|Sqrt|$\sqrt{画像}$ |
|Ln|自然対数|
|Exp|指数関数|
|AddSquare|$B=A^2+B$|$A,B$は画像|
|AddProduct|$C=A×B+C$|$A,B,C$は画像　バイアスありのアダマール積 |
|AddWeighted|$B=αA+(1-α)B$| $A,B$は画像　$α$はスカラ値|
|AbsDiff|Abs( $A-B$ )| $A,B$は画像|
|AbsDiffC|Abs( $A-α$ )| $A$は画像　$α$はスカラ値|

気づき
>* cuBLAS には アダマール積がなく、GPUのカーネル関数を書くのかと思いましたが、NPPのMulはまんまアダマール積ですね...。
* GPUのカーネル関数を書かなくても、NPPの関数の組み合わせれば、NN用の演算は実現できそうですね。

## 比較演算
|関数識別名|概要|備考|
|:-------|:-----|:-----|
|Threshold|各Pixelと閾値を比較演算し、真ならば、閾値に置き換える|比較演算は「より大きい」と「より小さい」|
|threshold_Val|各Pixelと閾値を比較演算し、真ならば、引数値に置き換える|比較演算は「より大きい」と「より小さい」|
|Threshold_LTValGTVal|各Pixelと２つの閾値(上限,下限)を比較演算し、真ならば、上限側,下限側のそれぞれの引数値に置き換える||

気づき
> * Threshold関数で ReLUができますね！
* LeakyReLUも算術演算と組み合わせでできますね。

## NPP 関数のサンプルコード
NPPのMUL関数のサンプルコードを示します。
内容

* float型 100×100 Pixel の2画像のMUL  
* 画像A 全て 2.0 の値 を初期値とする
* 画像B 全て 4.0 の値 を初期値とする
* MUL演算の結果はBに格納されてくる
* Bの値を表示し 値が 8.0 になっている事を確認する

```cpp:nppiMul_Sample.cpp
#include <cuda_runtime.h>
#include <npp.h>

int nppMul_example()
{
	int matMean = (100*100);
	size_t sz = sizeof(float) * matMean;
	float *hA , *hB;
	float *dA , *dB;

	hA = (float *)calloc(matMean,sizeof(float));
	if (hA == NULL){fprintf(stderr,"[%s L:%d]: calloc Error.\n",__FILE__,__LINE__);return -1;}

	hB = (float *)calloc(matMean,sizeof(float));
	if (hB == NULL){fprintf(stderr,"[%s L:%d]: calloc Error.\n",__FILE__,__LINE__);return -1;}

	for(int i=0;i<matMean;i++)
	{
		hA[i] = 2.0;
		hB[i] = 4.0;
	}

    if (cudaMalloc((void **)&dA, sz ) != cudaSuccess)
    {
        fprintf(stderr, "!!!! device memory allocation error (allocate A)\n");
        return -1;
    }

    if (cudaMalloc((void **)&dB, sz ) != cudaSuccess)
    {
        fprintf(stderr, "!!!! device memory allocation error (allocate B)\n");
        return -1;
    }

    // ホストのメモリからビデオメモリに配列をコピー
    cudaMemcpy( dA, hA, sz, cudaMemcpyHostToDevice);
    cudaMemcpy( dB, hB, sz, cudaMemcpyHostToDevice);

    // Mul
    NppiSize roi = {100,100};

    NppStatus nppRet;
    nppRet = nppiMul_32f_C1IR (
        dA,		//const Npp32f ∗ pSrc,
        100*sizeof(float),		//int nSrcStep,
        dB,		//Npp32f ∗ pSrcDst,
        100*sizeof(float),		//int nSrcDstStep,
        roi		//NppiSize oSizeROI
    );

	fprintf(stdout,"nppRet:%d\n",nppRet);

	// ビデオメモリからホストのメモリにデータをコピー
	cudaMemcpy( hA, dA, sz, cudaMemcpyDeviceToHost);
	cudaMemcpy( hB, dB, sz, cudaMemcpyDeviceToHost);

    // Aの確認　先頭の 10値を表示 
	printf("hA: ");
	for(int i=0;i<10;i++) { printf ("%2.2f ",hA[i]); }
	printf("\n");

    // B(結果)の確認　 
	printf("hB: ");
	for(int i=0;i<10000;i++)
	{
		if (i%100==0){printf("\n");}
		printf ("[%02d]%2.2f ",i%100,hB[i]);
	}
	printf("\n");

    // B(結果)は意図通り 全て 8.0 になっている。 

    cudaFree(dA);
    cudaFree(dB);
    free(hA);
    free(hB);

    return 0;
}

// 省略
```
補足
>* cuBLAS と同様に メモリ転送と NPP の目的の関数を呼び出すのみです。お手軽！
* GPU内のスレッド等は隠蔽されてます。


# 全結合NN の Affin の 行列演算
## 順伝播
　　![cr-06.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/92407302-20a6-80c8-fdd3-69334a2b7060.png)

* $M$ は前層からの順伝播の値 $m$個のパーセプトロン 
* $N$ は今層での Affin 後の値 $n$個のパーセプトロン
* 個々のパーセプトロンの式は、[以前の記事](http://qiita.com/t-tkd3a/items/9bf50f2e10e6a15b6ed5)にも記してあります。参考に...。

### オンライン 
　　$N=(Weight)M + Bias$

　　![cr-07.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/3a19d528-f859-219a-6163-e2e530db6719.png)

### ミニバッチ
　　$N=(Weight)M + Bias$
　　![cr-08.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/50dec098-11a4-4929-d29d-594262517d7b.png)


* cuBLASの[gemm関数](http://docs.nvidia.com/cuda/cublas/#cublas-lt-t-gt-gemm) $C=αAB+βC$ ( $A,B,C$ は行列、$α,β$はスカラ値)を用いる。
* $α=1.0$　　$β=1.0$
* $C$(上図$N$）は、初期値として$Bias$を ミニバッチ数の Broadcast を行っておく。

### ミニバッチ用行列 の 行優先メモリイメージ 
cuBLASでは メモリを列優先として扱うが、C言語では行優先で扱うため、上記行列(列優先）の行優先の表記も示します。

行優先表記
　　![cr-08_C.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/44a38d6b-b4be-6c69-eed3-7f13ef302928.png)

前記事での[(行優先・列優先)実装の話](http://qiita.com/t-tkd3a/items/098f51c2cc40a5bddc38#%E5%AE%9F%E8%A3%85%E3%81%AE%E8%A9%B1)も参考に...。

## 逆伝播

　　![cr-09.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/a97169a7-1fad-e1f8-69d6-66a1d27110c2.png)

* $dN$ は後続層からの逆伝播の入力値  
* $dM$ $dW$ $dB$ を算出する

### dM
　　$dM=(Weight)^TdN$

　　![cr-10.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/8f0e0434-f350-524f-a74a-7d211ca71530.png)

* $Weight$ は 転置 
* cuBLASの[gemm関数](http://docs.nvidia.com/cuda/cublas/#cublas-lt-t-gt-gemm) $C=αAB+βC$ ( $A,B,C$ は行列、$α,β$はスカラ値)を用いる。
* $α=1.0$　　$β=0.0$


### dM 行優先メモリイメージ
　　行優先表記
　　![cr-10_c.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/cf2bd0e8-ee07-c605-16d0-56365fc97de2.png)

### dW
　　$dW=dN(M)^T$


　　![cr-11.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/ff6c3959-551f-ac97-521e-266bca72053a.png)

* $M$ は 転置 
* cuBLASの[gemm関数](http://docs.nvidia.com/cuda/cublas/#cublas-lt-t-gt-gemm) $C=αAB+βC$ ( $A,B,C$ は行列、$α,β$はスカラ値)を用いる。
* $α=1.0$　　$β=0.0$


### dW 行優先メモリイメージ
　　行優先表記
　　![cr-11_c.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/74d51d9a-352c-f27b-8b37-f791b8c12381.png)



### dB
　　$dB=\sum^{ミニバッチ数}_{i=1} dN$

　　![cr-12.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/541b7f5d-e39b-83d9-5675-a7bf50646458.png)

* [gemv 関数](http://docs.nvidia.com/cuda/cublas/#cublas-lt-t-gt-gemv) を用いて 行列$dN$ の 行毎の総和 を算出する


### dB 行優先メモリイメージ
　　行優先表記
　　![cr-12_c.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/f69fa4b1-d3e1-97dc-62f0-6e1abb874f4c.png)