KerasのConv2Dの行列式演算

KerasでいうところのConv2Dがどのような演算をやっているかどういう風に理解してますか。
よくモデルの図解では直方体のデータ変形の例で示されますよね。
じゃあこれがどんな演算かっていうと初心者向け解説だと、畳み込みや特徴量抽出の説明をしてそれで終わってしまうことがままあります。
そうすると自分は、いや、畳み込みわかるけどConv2Dのパラメータ数とかどうなってるのか分らないよ、と思うわけです。

ということでConv2Dの自分なりの行列計算の理解をまとめたいと思います。もしかしたら詳しい参考書とかには書かれてるのかもしれませんが自分の読んだ本には行列計算は載っていませんでしたので……。
畳み込みや行列自体については解説しませんので畳み込みの原理は理解されたうえでお読みください。
もしも書いてる内容が間違っている場合は指摘していただけると幸いです。

```cnn.py
inputs = Input(shape=(28,28,1))                        # shape =(28,28,1)
x = Conv2D(32, (3, 3), activation='relu')(inputs)      # shape (28,28,1)=>(26,26,32)
x = Conv2D(64, (3, 3), activation='relu')(x)           # shape (26,26,32)=>(24,24,64)
...
```

#Dense演算におけるベクトル変形
さて、Conv2Dの解説をやると言っておきながら最初にニューラルネットワーク(NN)のDense演算でのベクトル次元の変形を行列式演算で説明します。
理由は先にやっているとConv2Dの行列演算の理解が進みやすいからです。
よくあるKerasでのMNISTのNNモデルを例にとって説明してみましょう。

```nn.py
inputs = Input(shape=(28,28))                         # shape =(28,28)
x1 = Flatten()(inputs)                                # shape (28,28)=>(784)
x2 = Dense(512, activation='relu')(x1)                # shape (784)=>(512)
x3 = Dense(512, activation='relu')(x2)                # shape (512)=>(512)
y = Dense(10, activation='softmax')(x3)               # shape (512)=>(10)
```
$(28,28)$のデータをFlattenによって$(784)$次元のベクトルにし、ベクトル長を$(784)→(512)→(512)→(10)$と変形させます。
この演算は行列式の数式でいうと

```math
    \left(
    \begin{array}{ccc}
      x2_{1} \\
      x2_{2} \\
      \vdots \\
      x2_{512}
    \end{array}
    \right)
   = \left(
    \begin{array}{ccc}
      A_{1,1} & A_{2,1} & \cdots & A_{784,1} \\
      A_{1,2} & A_{2,2} & \cdots & A_{784,2} \\
      \vdots  & \vdots  & \ddots & \vdots    \\
      A_{1,512} & A_{2,512} & \cdots & A_{784,512} 
    \end{array}
  \right)
    \left(
    \begin{array}{ccc}
      x1_{1} \\
      x1_{2} \\
      \vdots \\
      x1_{784}
    \end{array}
    \right)
     +
    \left(
    \begin{array}{ccc}
      B_{1} \\
      B_{2} \\
      \vdots \\
      B_{512}
    \end{array}
    \right)
```
```math
    \left(
    \begin{array}{ccc}
      x3_{1} \\
      x3_{2} \\
      \vdots \\
      x3_{512}
    \end{array}
    \right)
   = \left(
    \begin{array}{ccc}
      C_{1,1} & C_{2,1} & \cdots & C_{512,1} \\
      C_{1,2} & C_{2,2} & \cdots & C_{512,2} \\
      \vdots  & \vdots  & \ddots & \vdots    \\
      C_{1,512} & C_{2,512} & \cdots & C_{512,512} 
    \end{array}
  \right)
    \left(
    \begin{array}{ccc}
      x2_{1} \\
      x2_{2} \\
      \vdots \\
      x2_{512}
    \end{array}
    \right)
     +
    \left(
    \begin{array}{ccc}
      D_{1} \\
      D_{2} \\
      \vdots \\
      D_{512}
    \end{array}
    \right)
```
```math
    \left(
    \begin{array}{ccc}
      y_{1} \\
      y_{2} \\
      \vdots \\
      y_{10}
    \end{array}
    \right)
   = \left(
    \begin{array}{ccc}
      E_{1,1} & E_{2,1} & \cdots & E_{512,1} \\
      E_{1,2} & E_{2,2} & \cdots & E_{512,2} \\
      \vdots  & \vdots  & \ddots & \vdots    \\
      E_{1,10} & E_{2,10} & \cdots & E_{512,10} 
    \end{array}
  \right)
    \left(
    \begin{array}{ccc}
      x3_{1} \\
      x3_{2} \\
      \vdots \\
      x3_{512}
    \end{array}
    \right)
     +
    \left(
    \begin{array}{ccc}
      F_{1} \\
      F_{2} \\
      \vdots \\
      F_{10}
    \end{array}
    \right)
```
という計算を行います。(活性化関数の計算は省略しています)
この$A～F$の行列の値がモデルの重みであり、NN学習では損失を小さくするようなモデル重みを学習します。
ところで一番上の行列式は以下のようにまとめることが可能です。

```math
    \left(
    \begin{array}{ccc}
      x2_{1} \\
      x2_{2} \\
      \vdots \\
      x2_{512}
    \end{array}
    \right)
   = \left(
    \begin{array}{ccc}
      A_{1,1} & A_{2,1} & \cdots & A_{784,1} & B_{1} \\
      A_{1,2} & A_{2,2} & \cdots & A_{784,2} & B_{2} \\
      \vdots  & \vdots  & \ddots & \vdots    & \vdots \\
      A_{1,512} & A_{2,512} & \cdots & A_{784,512} & B_{512}
    \end{array}
  \right)
    \left(
    \begin{array}{ccc}
      x1_{1} \\
      x1_{2} \\
      \vdots \\
      x1_{784} \\
      1
    \end{array}
    \right)
```
```math
　　x2 = A'\cdot x1'
```
つまり入力次元$M$次元、出力次元$N$次元のベクトル変形には$(M+1)\cdot N$のパラメータ数が必要です。

#Conv2D演算における直方体変形
ではKerasにおける記述にて以下のようなConv2Dの演算によって直方体データの形が$(26,26,32)→(24,24,64)$に変形する行列演算を考えてみましょう。

```cnn.py
inputs = Input(shape=(28,28,1))                        # shape =(28,28,1)
x1 = Conv2D(32, (3, 3), activation='relu')(inputs)     # shape (28,28,1)=>(26,26,32)
x2 = Conv2D(64, (3, 3), activation='relu')(x1)         # shape (26,26,32)=>(24,24,64)
...
```
このとき`x2 = Conv2D(64, (3, 3), activation='relu')(x1)`に相当するConv2D演算は

```math
    \left(
    \begin{array}{ccc}
      x2_{channel1}[24×24] \\
      x2_{channel2}[24×24] \\
      \vdots \\
      x2_{channel64}[24×24]
    \end{array}
    \right)
   = \left(
    \begin{array}{ccc}
      A_{1,1}[3×3]\ast & A_{2,1}[3×3]\ast & \cdots & A_{32,1}[3×3]\ast \\
      A_{1,2}[3×3]\ast & A_{2,2}[3×3]\ast & \cdots & A_{32,2}[3×3]\ast \\
      \vdots  & \vdots  & \ddots & \vdots    \\
      A_{1,64}[3×3]\ast & A_{2,64}[3×3]\ast & \cdots & A_{32,64}[3×3]\ast 
    \end{array}
  \right)
    \left(
    \begin{array}{ccc}
      x1_{channel1}[26×26] \\
      x1_{channel2}[26×26] \\
      \vdots \\
      x1_{channel32}[26×26]
    \end{array}
    \right)
     +
    \left(
    \begin{array}{ccc}
      B_{1} \\
      B_{2} \\
      \vdots \\
      B_{64}
    \end{array}
    \right) \cdot E [24×24]
```
という数式で表すことができます。
また出力の1チャンネル目に限れば以下のように畳み込みの線形和で記述できます。

```math
x2_{channel1}[24×24]　= A_{1,1}[3×3] \ast x1_{channel1}[26×26] + A_{2,1}[3×3] \ast x1_{channel2}[26×26] + \cdots + A_{32,1}[3×3] \ast x1_{channel32}[26×26] + B_{1} \cdot E [24×24]
```

ここで$\ast$は畳み込み演算子、$A_{i,j}[3×3]$は$(3,3)$サイズの畳み込みフィルター、$E[24×24]$はすべての要素の値が$1$の$(24,24)$サイズの行列です。

###DenseとConv2Dの比較
Denseのベクトル変形とConv2Dにおける直方体変形を比較してみましょう。
Denseのベクトル次元数はConv2Dではチャンネル数に相当することが分かります。
また変換行列$A$における各要素$A_{i,j}$はDenseの場合は一つの数値ですが、Conv2Dの場合はこれは$(3,3)$サイズの畳み込みフィルターになります。
ここで注意していただきたいのはConv2Dにおける畳み込みフィルターの総数(パラメータ数)は入力チャンネル数$32$と出力チャンネル数$64$の積に依存し変換画像の大きさには依存しないということです。
(無論、入力画像が大きくなれば**計算所要時間は増えますが**)
つまり、畳み込みフィルターサイズを$(L,L)$、入力チャンネルを$M$、出力チャンネルを$N$とすると直方体変形には$(L\cdot L\cdot M\cdot N + N)$のパラメータ数が必要です。

最後に実際にVGG16におけるパラメータ数を確認してみましょう。
$(224,224,3)→(224,224,64)$の変形では3*3*3*64+64=1792ですし、
$(224,224,64)→(224,224,64)$の変形では3*3*64*64+64=36928となっていて
Conv2Dのパラメータ数は画像サイズとは関係ないことが分かります。

```
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
input_1 (InputLayer)         (None, 224, 224, 3)       0         
_________________________________________________________________
block1_conv1 (Conv2D)        (None, 224, 224, 64)      1792      
_________________________________________________________________
block1_conv2 (Conv2D)        (None, 224, 224, 64)      36928     
_________________________________________________________________
block1_pool (MaxPooling2D)   (None, 112, 112, 64)      0         
_________________________________________________________________
……
```
#考察
- モデル図解での直方体データから直方体データへの変形がConv2Dやpoolingである。直方体自体がConv2Dではないんですよね。当たり前かもしれませんけど。
- ところが、下のVGG16モデルの図のようにさも直方体本体がConv2Dみたいな図の書かれ方をされることがあります。なんでだろう？自分の理解だと直方体間の隙間がConv2Dに相当するのだが。
<img src="https://qiita-image-store.s3.amazonaws.com/0/299928/0091c4ac-0751-afd8-2560-c4e102386b2e.png"  width=60%>

- Kerasだと`x2 = Conv2D(`**`64`**`, (3, 3), activation='relu')(x1) `における64に相当する出力チャンネル数を指してフィルター数と呼びます。([Convolutionalレイヤー - Keras Documentation](https://keras.io/ja/layers/convolutional/))
しかし、行列計算上における畳み込みのフィルターの総数は入力チャンネル、出力チャンネル数の積ですから、今回の解釈に関してこれでよいのかっていう疑問はあります。
- 先に述べたようにConv2Dのパラメータ数は入力チャンネル、出力チャンネル数の積に依存します。すべてのフィルター数を2倍にするとパラメータ数は4倍に、すべてのフィルター数を4倍にするとパラメータ数は16倍になるわけですから、モデル作成においてはフィルター数の増加するより並列に畳み込み層を追加するほうがパラメータ数の増加を抑えられます。
- Conv2Dのパラメータ数は入力チャンネル、出力チャンネル数の積に依存し、画像サイズに依存しません。これは入力画像サイズの異なるモデルの再学習においてフィルター数が同じであるならばConv2Dのモデル重みは入力画像サイズが異なってもパラメータ数は等しいため使いまわせるのだろうと推測されます。ただし、全結合におけるベクトル次元は画像サイズに依存するので、入力画像サイズが異なる場合、全結合部分のモデル重みは使いまわせません。



#まとめ
直方体データ変形といっても実際にどんな演算すれば三次元の直方体データを変形できるのかっていうのは意外にイメージしにくいのではないかと思います。おそらく三次元配列を計算する機会があまりありませんから。
今回はConv2D演算を行列式にて表記してみました。
データを直方体で表したときConv2Dは1,2次元目(縦横方向)に関して畳み込みを行い、3次元目(チャンネル方向)には全結合を行っているのに感覚的に近いかと思いました。

#おまけ:SeparableConv2Dはどうなるの？
どうせならKerasのSeparableConv2Dも行列式で書いてみましょう。
SeparableConv2DとはDepthwise convolutionとPointwise convolution((1×1)畳み込み)をまとめて行うのをこう呼ぶらしいです。SeparableConv2Dの形は下記記事を参考にしました。
[tf.nn.conv2d, tf.nn.depthwise_conv2d, tf.nn.separable_conv2dのチャネル数と意味](https://qiita.com/YusukeSuzuki@github/items/0764d15b9d0b97ec1e16)
[Kerasの作者@fcholletさんのCVPR'17論文XceptionとGoogleのMobileNets論文を読んだ](https://qiita.com/yu4u/items/34cd33b944d8bdca142d)
[MobileNets: CNNのサイズ・計算コストの削減手法＿翻訳・要約](https://qiita.com/HiromuMasuda0228/items/7dd0b764804d2aa199e4)

この時、`x2 = SeparableConv2D(64, (3, 3), depth_multiplier=3, activation='relu')(x1)`で示されるdepth_multiplier=3のSeparableConv2Dの行列式演算は以下のように書けます。
畳み込みフィルターの数はdepth_multiplierと入力チャンネル数の積に比例するのでConv2Dより少ないパラメータでConv2Dを表現することになります。

```math
    \left(
    \begin{array}{ccc}
      x2_{channel1}[24×24] \\
      x2_{channel2}[24×24] \\
      \vdots \\
      x2_{channel64}[24×24]
    \end{array}
    \right)
   =
   \left(
    \begin{array}{ccc}
      Pw_{1,1} & Pw_{2,1} & \cdots & Pw_{96,1} & B_{1}\\
      Pw_{1,2} & Pw_{2,2} & \cdots & Pw_{96,2} & B_{2}\\
      \vdots  & \vdots  & \ddots & \vdots  & \vdots  \\
      Pw_{1,64} & Pw_{2,64} & \cdots & Pw_{96,64} & B_{64}
    \end{array}
  \right) 
  \left(
    \begin{array}{ccc}
      Dw_{1,1}[3×3]\ast & 0 & \cdots & 0 & 0\\
      Dw_{1,2}[3×3]\ast & 0 & \cdots & 0 & 0\\
      Dw_{1,3}[3×3]\ast & 0 & \cdots & 0 & 0\\
      0 & Dw_{2,4}[3×3]\ast & \cdots & 0 & 0\\
      0 & Dw_{2,5}[3×3]\ast & \cdots & 0 & 0\\
      0 & Dw_{2,6}[3×3]\ast & \cdots & 0 & 0\\
      \vdots  & \vdots  & \ddots & \vdots  & \vdots  \\
      0 & 0 & \cdots & Dw_{32,94}[3×3]\ast & 0 \\
      0 & 0 & \cdots & Dw_{32,95}[3×3]\ast & 0 \\
      0 & 0 & \cdots & Dw_{32,96}[3×3]\ast & 0 \\ 
      0 & 0 & \cdots & 0 & 1
    \end{array}
  \right)
    \left(
    \begin{array}{ccc}
      x1_{channel1}[26×26] \\
      x1_{channel2}[26×26] \\
      \vdots \\
      x1_{channel32}[26×26] \\
      E[24×24]
    \end{array}
    \right)

```
ここで畳み込みフィルターサイズを$(L,L)$、入力チャンネルを$M$、出力チャンネルを$N$、depth_multiplierを$O$とするとSeparableConv2D直方体変形には$(M \cdot O + 1) \cdot N +((L \cdot L)\cdot M \cdot O)$のパラメータ数が必要です。

実際にモデルのパラメータ数を確認すると
$(26,26,32)→(24,24,64)$のSeparableConv2Dの計算でパラメータ数は(32*3+1)\*64+((3*3)*32*3)=7072で合っているのが分かります。

```sep_cnn.py
inputs = Input(shape=(28,28,1))                                                     # shape =(28,28,1)
x1 = SeparableConv2D(32, (3, 3), depth_multiplier=3, activation='relu')(inputs)     # shape (28,28,1)=>(26,26,32)
x2 = SeparableConv2D(64, (3, 3), depth_multiplier=3, activation='relu')(x1)         # shape (26,26,32)=>(24,24,64)
...
...
_________________________________________________________________
Layer (type)                 Output Shape              Param #
=================================================================
input_1 (InputLayer)         (None, 28, 28, 1)         0
_________________________________________________________________
separable_conv2d_1 (Separabl (None, 26, 26, 32)        155
_________________________________________________________________
separable_conv2d_2 (Separabl (None, 24, 24, 64)        7072
_________________________________________________________________
...
```