計算するって何だろう（for ディープラーニング）

## はじめに
ディープラーニングを学習する過程で、計算についてのメンタルモデル——計算とはこういうものだというイメージ——が拡張されるのを感じました。

本記事は、その新しいメンタルモデルを共有するために書かれました。

## 計算とは
世界で一番易しい計算、といえば以下でしょう。

$$
1 + 1 =\,?
$$

この易しい計算を難しく分解すると以下となります。

- 構造：$y = f(x_0, x_1) = x_0 + x_1$
- 要素：$x_0 = 1, x_1 = 1, y = unknown$

そして、計算とは、構造と要素（一部が未知数）が与えられた状態で、**構造に「フィット」する要素の組を求めること**と定義できます。

上述の計算は、順方向——左辺から右辺——へと、一度だけ行われます。

世の中の大多数の人々にとって、計算とは、順方向へと一度だけ行われるものでしょう。しかし、そうではない計算を考えることもできます。

## 双方向に、何度でも（１）
例えば、以下のような計算を考えます。

$$
?\,+\,?\,= 2
$$

- 構造：$y = f(x_0, x_1) = x_0 + x_1$
- 要素：$x_0 = unknown, x_1 = unknown, y = 2$

構造は先ほどの「世界で一番易しい計算」と同じですが、要素が2点異なります。

- 左辺が未知数であること
- 要素の組が一意に定まりそうにないこと

以上の計算に対し、以下の3つのアクションを組み合わせて、計算を行うものとします。

- 初期値の更新
- 順方向への計算
- 逆方向へのフィードバック

例えば、こんなアプローチが考えられるでしょう。

1. 未知数をランダムに初期化する
1. 順方向への計算を行う
1. 計算結果と真の値の誤差が十分小さければ終了する
1. 逆方向へとフィードバックする（未知数をランダムに初期化する）
1. 2に戻る

計算は、双方向に何度も行われます。誤差は、以下の誤差関数eで表すものとします。

$$
e = \frac{1}{2}\{y - (x_0 + x_1)\}^2
$$

二乗しているのは正負の違いを吸収するため、二分の一倍しているのは、のちの都合上です。

以下がコードの例です。


```python
import random

def randint():
    return random.random() * 10 - 5

def calc():
    x0, x1, y = randint(), randint(), 2
    delta = 0.0001
    counter = 0
    while True:
        counter += 1
        # forward
        out = x0 + x1
        
        # backward
        error = 0.5 * (y - out) ** 2
        if error < delta: break
        dout = randint(), randint()
        x0, x1 = dout[0], dout[1]
    message = '{}回の計算で、x0が{:.3f}、x1が{:.3f}と判明。計算結果は{:.3f}です。'
    print(message.format(counter, x0, x1, x0 + x1))

for _ in range(10): calc()
```

    463回の計算で、x0が-1.724、x1が3.728と判明。計算結果は2.004です。
    1回の計算で、x0が4.247、x1が-2.259と判明。計算結果は1.987です。
    99回の計算で、x0が0.465、x1が1.529と判明。計算結果は1.993です。
    230回の計算で、x0が-1.126、x1が3.137と判明。計算結果は2.010です。
    983回の計算で、x0が0.227、x1が1.784と判明。計算結果は2.011です。
    143回の計算で、x0が2.904、x1が-0.902と判明。計算結果は2.003です。
    22回の計算で、x0が1.947、x1が0.066と判明。計算結果は2.013です。
    112回の計算で、x0が2.865、x1が-0.872と判明。計算結果は1.992です。
    1338回の計算で、x0が-2.327、x1が4.338と判明。計算結果は2.011です。
    83回の計算で、x0が3.545、x1が-1.535と判明。計算結果は2.010です。


確かに（近似値レベルではありますが）計算できており、アプローチが有効であることがわかります。

しかし、以下の点を改善できないものでしょうか。

- 計算回数が多い
- 計算回数がバラつく

## 双方向に、何度でも（２）
結論的には、以下アプローチにより、計算回数を少なく、かつ、安定させることができます。

1. 未知数をランダムに初期化する
1. 順方向への計算を行う
1. 計算結果と真の値の誤差が十分小さければ終了する
1. 逆方向へとフィードバックする（要素を誤差関数の偏微分結果で更新する）
1. 2に戻る

微分が出てくるのは、それが誤差を小さくするための合理的な方法だからです。というのも、誤差関数eを微分することで、勾配がわかるので、要素を増やせば良いのか、減らせば良いのかわかるのです。

数式は以下のとおり。

$$
e = \frac{1}{2}\{y - (x_0 + x_1)\}^2
$$

$$
\frac{\partial e}{\partial x_0} = \frac{\partial e}{\partial x_1} = -y + x_0 + x_1
$$

以下がコードの例です。


```python
import random

def randint():
    return random.random() * 10 - 5

def calc():
    x0, x1, y = randint(), randint(), 2
    alpha, delta = 0.1, 0.0001
    counter = 0
    while True:
        counter += 1
        # forward
        out = x0 + x1
        
        # backward
        error = 0.5 * (y - out) ** 2
        if error < delta: break
        dout = (-y + x0 + x1, -y + x0 + x1)
        x0, x1 = x0 - alpha * dout[0], x1 - alpha * dout[1]
    message = '{}回の計算で、x0が{:.3f}、x1が{:.3f}と判明。計算結果は{:.3f}です。'
    print(message.format(counter, x0, x1, x0 + x1))

for _ in range(10): calc()
```

    19回の計算で、x0が-1.544、x1が3.556と判明。計算結果は2.012です。
    30回の計算で、x0が0.149、x1が1.838と判明。計算結果は1.988です。
    30回の計算で、x0が2.167、x1が-0.179と判明。計算結果は1.987です。
    17回の計算で、x0が2.103、x1が-0.089と判明。計算結果は2.014です。
    24回の計算で、x0が3.955、x1が-1.968と判明。計算結果は1.988です。
    11回の計算で、x0が4.346、x1が-2.335と判明。計算結果は2.012です。
    29回の計算で、x0が-0.809、x1が2.797と判明。計算結果は1.988です。
    21回の計算で、x0が1.667、x1が0.347と判明。計算結果は2.014です。
    14回の計算で、x0が4.051、x1が-2.065と判明。計算結果は1.986です。
    28回の計算で、x0が0.947、x1が1.039と判明。計算結果は1.986です。


計算回数が少なく、かつ、安定していることがわかります。

乗算についても、見てみましょう。数式は以下のとおり。

$$
e = \frac{1}{2}\{y - (x_0 \cdot x_1)\}^2
$$

$$
\frac{\partial e}{\partial x_0} = -y \cdot x_1 + {x_1}^2 \cdot x_0
$$

$$
\frac{\partial e}{\partial x_1} = -y \cdot x_0 + {x_0}^2 \cdot x_1
$$

以下がコードの例です。


```python
import random

def randint():
    return random.random() * 10 - 5

def calc():
    x0, x1, y = randint(), randint(), 2
    alpha, delta = 0.1, 0.0001
    counter = 0
    while True:
        counter += 1
        # forward
        out = x0 * x1
        
        # backward
        error = 0.5 * (y - out) ** 2
        if error < delta: break
        dout = (-y * x1 + (x1 ** 2) * x0, -y * x0 + (x0 ** 2) * x1)
        x0, x1 = x0 - alpha * dout[0], x1 - alpha * dout[1]
    message = '{}回の計算で、x0が{:.3f}、x1が{:.3f}と判明。計算結果は{:.3f}です。'
    print(message.format(counter, x0, x1, x0 * x1))

for _ in range(10): calc()
```

    19回の計算で、x0が1.504、x1が1.320と判明。計算結果は1.986です。
    21回の計算で、x0が1.338、x1が1.488と判明。計算結果は1.990です。
    8回の計算で、x0が2.364、x1が0.841と判明。計算結果は1.988です。
    17回の計算で、x0が-1.661、x1が-1.198と判明。計算結果は1.991です。
    18回の計算で、x0が1.238、x1が1.607と判明。計算結果は1.989です。
    17回の計算で、x0が1.266、x1が1.571と判明。計算結果は1.988です。
    13回の計算で、x0が1.468、x1が1.357と判明。計算結果は1.991です。
    10回の計算で、x0が1.486、x1が1.353と判明。計算結果は2.010です。
    19回の計算で、x0が1.305、x1が1.524と判明。計算結果は1.989です。
    3回の計算で、x0が0.616、x1が3.265と判明。計算結果は2.010です。


加算と同様、計算回数が少なく、かつ、安定していることがわかります。

## ディープラーニング
ディープラーニングで、レイヤーを重ねてネットワークを定義するとき、計算の構造が定まります。また、計算の要素のうち、左辺の入力値と、右辺の出力値も定まります（教師あり学習）。しかし、その他の要素（重みやバイアス）については、未知数のままです。これをランダムに初期化し、双方向に何度も計算を行い、構造にフィットする要素の組を求める、ということが行われます。構造にフィットする要素の組は複数あるので、そのうちの1つを求めているイメージです。

[Keras Blog](https://blog.keras.io/building-autoencoders-in-keras.html)を参考に、-100〜100の実数xを入力、および、出力するオートエンコーダ（自己符号化器）を実装してみます。xをエンコードし、すぐにデコードします。活性化関数は用いません。数式で表すと、以下の構造にフィットする重み（Weight）とバイアス（Bias）の組を求める問題です。とても単純な例です。

$$
w_{dec} \cdot (w_{enc} \cdot x + b_{enc}) + b_{dec} = x
$$

以下に、コードと実行結果を示します。


```python
from keras.layers import Input, Dense
from keras.models import Model

# this is the size of our encoded representations
encoding_dim = 1

# this is our input placeholder
input_img = Input(shape=(1,))
# "encoded" is the encoded representation of the input
encoded = Dense(encoding_dim, activation=None)(input_img)
# "decoded" is the lossy reconstruction of the input
decoded = Dense(1, activation=None)(encoded)

# this model maps an input to its reconstruction
autoencoder = Model(input_img, decoded)

# this model maps an input to its encoded representation
encoder = Model(input_img, encoded)

# create a placeholder for an encoded (32-dimensional) input
encoded_input = Input(shape=(encoding_dim,))
# retrieve the last layer of the autoencoder model
decoder_layer = autoencoder.layers[-1]
# create the decoder model
decoder = Model(encoded_input, decoder_layer(encoded_input))

print(autoencoder.summary())
autoencoder.compile(optimizer='nadam', loss='mean_squared_error')

def load_data():
    data = np.arange(-100, 100, 0.01)
    np.random.shuffle(data)
    border = len(data) * 9 // 10
    return data[: border], data[border :]

import numpy as np
x_train, x_test = load_data()

print(x_train.shape)
print(x_test.shape)

autoencoder.fit(x_train, x_train,
                epochs=30,
                batch_size=256,
                shuffle=True,
                validation_data=(x_test, x_test))

# encode and decode some digits
# note that we take them from the *test* set
encoded_imgs = encoder.predict(x_test)
decoded_imgs = decoder.predict(encoded_imgs)

n = 5  # how many digits we will display

print('-------------')
for i in range(n):
    # display original
    print(x_test[i])
    
    # display encoded representation
    print(encoded_imgs[i][0])
    
    # display reconstruction
    print(decoded_imgs[i][0])
    
    print('-------------')
    
autoencoder.get_weights()
```

    _________________________________________________________________
    Layer (type)                 Output Shape              Param #   
    =================================================================
    input_1 (InputLayer)         (None, 1)                 0         
    _________________________________________________________________
    dense_1 (Dense)              (None, 1)                 2         
    _________________________________________________________________
    dense_2 (Dense)              (None, 1)                 2         
    =================================================================
    Total params: 4
    Trainable params: 4
    Non-trainable params: 0
    _________________________________________________________________
    None
    (18000,)
    (2000,)
    Train on 18000 samples, validate on 2000 samples
    Epoch 1/30
    18000/18000 [==============================] - 0s 13us/step - loss: 2567.6244 - val_loss: 1672.4072
    Epoch 2/30
    18000/18000 [==============================] - 0s 4us/step - loss: 1124.5867 - val_loss: 711.7286
    Epoch 3/30
    18000/18000 [==============================] - 0s 4us/step - loss: 471.2158 - val_loss: 291.2063
    〜中略〜
    Epoch 28/30
    18000/18000 [==============================] - 0s 4us/step - loss: 1.7792e-08 - val_loss: 1.7726e-08
    Epoch 29/30
    18000/18000 [==============================] - 0s 4us/step - loss: 1.7062e-08 - val_loss: 1.6747e-08
    Epoch 30/30
    18000/18000 [==============================] - 0s 4us/step - loss: 1.5183e-08 - val_loss: 1.4879e-08
    -------------
    32.99000000006805
    -35.798763
    32.990067
    -------------
    -33.85999999996616
    36.74323
    -33.860077
    -------------
    -95.41999999999766
    103.54479
    -95.4202
    -------------
    -86.44999999999307
    93.811035
    -86.45018
    -------------
    -27.51999999996292
    29.863403
    -27.52006
    -------------

    [array([[-1.0851457]], dtype=float32),
     array([0.00019342], dtype=float32),
     array([[-0.92153716]], dtype=float32),
     array([0.00017642], dtype=float32)]



5回実行し、得られた重みとバイアスは、以下のとおりです。

|$w_{enc}$|$b_{enc}$|$w_{dec}$|$b_{dec}$|
|:-:|:-:|:-:|:-:|
|-0.52114826|0.0008159|-1.9188374|0.0015671|
|1.6570777|0.01260246|0.6034708|-0.00760303|
|-1.3857455|-0.01825085|-0.72163296|-0.01317422|
|-0.6569005|-0.00015455|-1.5222994|-0.00023566|
|-1.0851457|0.00019342|-0.92153716|0.00017642]|

構造にフィットする要素の組が複数ある、ということがわかります。

## まとめ

どうやら、ディープラーニングを理解するためには、計算についてのメンタルモデルを、以下のとおりに拡張する必要がありそうです。

まず、大前提として、計算について、以下のとおりに定義します。

- 計算とは、構造と要素（一部が未知数）が与えられた状態で、構造に「フィット」する要素の組を求めることである

その上で、メンタルモデルを以下のとおりに拡張します。

- 計算とは、順方向に一度だけ行われるものだ
- 構造にフィットする要素の組は、一意に定まるものだ

↓

- 計算は、双方向に何度も行われることがある
  - フィッティング途上の要素を用いて順方向に計算
  - 逆方向へとフィードバック（要素を誤差関数の偏微分結果で更新）
- 構造にフィットする要素の組が、複数あることもある

最終的に等号が成立する点においては同じなのですが、過程が異なるのです。

以上を示して、本記事を終えたいと思います。