【Colab】「ニューラルネットワーク、多様体、トポロジー」を実装する

Christopher Olah氏のブログ記事[「ニューラルネットワーク、多様体、トポロジー」](http://colah.github.io/posts/2014-03-NN-Manifolds-Topology/)にインスパイアされて、Olah氏の記事に記載されていたいくつかの動画を実装してみました。


# 「ニューラルネットワーク、多様体、トポロジー」記事概要

Google Brain、OpenAI、AnthropicAIなど機械学習の最前線で活躍するChiritopher Olah氏が2014年に投降したブログ記事「ニューラルネットワーク、多様体、トポロジー」が面白いです。
深層学習をはじめとするニューラルネットワークによる機械学習がある層から次の層へと進むにつれて「何をやっているのか」、結果として「なぜ優れたパターン認識性能を生み出せるのか」を視覚的に理解しようとした記事になります。（記事自体は多様体学習にまで踏み込んだより深い内容ですが、この記事では扱いません。というかあつかえませんでしｔ）

例えば、下記のように青色の線と赤色の線を識別したい「二値分類問題」をニューラルネットで解く場合に、インプットされた入力値をどのように分離可能な状態へと変換していくかの過程が、視覚的に直観的に理解できるような動画がいくつかアップされていています。

![687474703a2f2f636f6c61682e6769746875622e696f2f706f7374732f323031342d30332d4e4e2d4d616e69666f6c64732d546f706f6c6f67792f696d672f73706972616c2e312d322e322d322d322d322d322d322e676966.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/107669/438fc7af-b36d-fcdc-e822-18311f50e36e.gif)

記事のソースとなるコードはないみたいなので、記事の内容から「恐らくこういうことをやっているんだろうな」という推定ベースで再現してみました。

より詳細な記事内容については、@KojiOhkiさんが翻訳して下さっていたので[そちら](https://qiita.com/KojiOhki/items/af2241027b00f892d2bd)も参考にしました。

https://qiita.com/KojiOhki/items/af2241027b00f892d2bd

# 実装
Google Colab（python3系）上で実装しました。ニューラルネット構築に際しては、モジュールを使わずnumpyだけで書いてみました。

#### モジュールのインポート
```python
import random
import numpy                   as np
import matplotlib.pyplot as plt
from matplotlib.animation import PillowWriter
from matplotlib import animation, rc
from IPython.display import HTML
```
僕個人のドライブへのマウントなども設定していますが、そちらは省略しました。


#### 活性化関数ほか関数群の定義
ニューラルネットは「線形（アフィン）変換」→「活性化関数による写像」→「線形変換」→・・・と続いていくので、線形変換の関数と、（いくつかの）活性化関数を用意しました。また、誤差逆伝播による学習をしていくので、微分関数も用意しています。

最後の方で出力層のための活性化関数と誤差関数も用意しています。今回は2値問題なので、出力層への活性関数としてsigmoid関数、誤差関数として交差エントロピーと利用しました。交差エントロピーとsigmoid関数の合成関数はその微分値が綺麗な形になるため、合成関数の微分関数として記述しています。

```python
''' アフィン変換 '''
def affine(x, W, b):
    return np.dot(x, W) + b

# 勾配
def affineBack(du, x, W, b):
    dx = np.dot(du, W.T)
    dW = np.dot(x.T, du)
    db = np.sum(du, axis=0)    
    return dx, dW, db
    
''' ReLU関数 '''
def relu(x):
    return np.maximum(0.0, x)

# 勾配
def reluBack(u):
    return np.where(u > 0.0, 1.0, 0.0)
    
''' tanh関数 '''
def tanh(x):
    return ( np.exp(x) - np.exp(-x) ) / ( np.exp(x) + np.exp(-x) )

def tanhBack(u):
    return  (1.0 - tanh(u)**2.0)

''' シグモイド関数 '''
def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))

''' 誤差関数：交差エントロピー'''
# 誤差(交差エントロピー）＋活性化関数(sigmoid)勾配
def sigmoidCrossEntropyErrorBack(y, t):
    return (y - t) / y.shape[0]

def crossEntropyError(y, t):
    delta = 1e-7  #log 0対策
    return - np.sum(t * np.log(y + delta)) / y.shape[0]
```

#### ニューラルネットクラスの定義
次にNeuralNetworkクラスを定義します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/107669/e4672e94-8903-98a1-6f23-e2ab41e12b74.png)


プーリングなどはないシンプルな全結合ニューラルネットワークになります。
元記事では画像ごとに隠れ層の数が変更されているみたいですが、ここでは隠れ層は3層でfixしてみました。各層のノード数は変更可能な形で実装していますが、後ほどハイパーパラメータとしてノード数は固定します。

```python
class NeuralNetwork:
    ''' 初期値 '''
    def __init__(self, ISIZE, H1SIZE, H2SIZE, H3SIZE, OSIZE, LEARNRATE):
        ''' 活性化関数・およびその微分関数の決定 '''
        self.activH1     = tanh
        self.activH2     = tanh
        self.activH3     = tanh
        self.activH1Back = tanhBack
        self.activH2Back = tanhBack
        self.activH3Back = tanhBack
        self.activO      = sigmoid  
        self.activOBack  = sigmoidCrossEntropyErrorBack   

        ''' 各層のユニット数 '''
        self.iSize  = ISIZE
        self.h1Size = H1SIZE
        self.h2Size = H2SIZE
        self.h3Size = H3SIZE
        self.oSize  = OSIZE 

        ''' 重み・バイアス初期値 '''
        self.W1 = np.random.rand(self.iSize , self.h1Size) * 2.0 - 1.0
        self.W2 = np.random.rand(self.h1Size, self.h2Size) * 2.0 - 1.0
        self.W3 = np.random.rand(self.h2Size, self.h3Size) * 2.0 - 1.0
        self.W4 = np.random.rand(self.h3Size, self.oSize ) * 2.0 - 1.0

        self.b1 = np.zeros(self.h1Size) + 0.01
        self.b2 = np.zeros(self.h2Size) + 0.01
        self.b3 = np.zeros(self.h3Size) + 0.01
        self.b4 = np.zeros(self.oSize ) + 0.01

    ''' 順伝播：入力信号　⇒　隠れ層　⇒　出力 '''
    def forward(self, Input):
        u1 = affine(Input, self.W1, self.b1)
        z1 = self.activH1(u1)
        u2 = affine(z1, self.W2, self.b2)
        z2 = self.activH2(u2)
        u3 = affine(z2, self.W3, self.b3)
        z3 = self.activH3(u3)
        u4 = affine(z3, self.W4, self.b4)
        y  = self.activO(u4)        
        
        return u1, z1, u2, z2, u3, z3, u4, y

    ''' 逆伝播：SGD '''
    def backward(self, Input, CorrectD):
        u1, z1, u2, z2, u3, z3, u4, y = self.forward(Input)
  
        du4           = self.activOBack(y, CorrectD)
        dz3, dW4, db4 = affineBack(du4, z3, self.W4, self.b4)
        du3           = dz3 * self.activH3Back(u3)
        dz2, dW3, db3 = affineBack(du3, z2, self.W3, self.b3)
        du2           = dz2 * self.activH2Back(u2)
        dz1, dW2, db2 = affineBack(du2, z1, self.W2, self.b2)
        du1           = dz1 * self.activH1Back(u1)
        dx,  dW1, db1 = affineBack(du1, Input, self.W1, self.b1)        
        
        # 重み更新：右辺の重みはt-1期
        self.W1   -= LEARNRATE * dW1
        self.b1   -= LEARNRATE * db1
        self.W2   -= LEARNRATE * dW2
        self.b2   -= LEARNRATE * db2
        self.W3   -= LEARNRATE * dW3
        self.b3   -= LEARNRATE * db3
        self.W4   -= LEARNRATE * dW4
        self.b4   -= LEARNRATE * db4
```
各層の活性化関数はtanh関数にしました。（後ほど活性化関数をいじってみたいと思います。）

あまり拡張性のある形で記述できなかったのは少し悔しいですが、今回の検討範囲ではこれで十分かなと思っています。重みパラメータの初期値なども適当です。対象とする問題が簡単なものになるので、慣性項やAdaGradなど学習効率改善のためのパラメータはセットしていません。

#### 学習のための関数
学習を実行する関数を用意します。
```python
def learning(Data):
  # 学習
  for i in range(ITER):
    # 学習：バッチサイズ
    index = np.random.choice(range(Data.shape[0]), BATCH)
    input = Data[index, :][:,[0,1]]
    t     = Data[index, :][:,[2,3]]

    NN.backward(input, t)

    # 誤差経過観察
    if i % 10000 == 0:
      u1, z1, u2, z2, u3, z3, u4, out  = NN.forward(input)

      error = crossEntropyError(out, t)
      print("%d th: %f "%(i, error))
```
バッチ毎にデータからランダムでサンプルをピックアップして学習します。定期的に誤差がどの程度減っているかを可視化するようにもしています。


#### グローバル変数（ハイパーパラメータ含む）の定義
```python
N         = 5000        # 各要素データ数
GSIZE     = 50          # 曲面可視化のためのグリッド点
ISIZE     = 2           # 入力サイズ（x, y座標）
H1SIZE    = 20          # 隠れ層のユニット数
H2SIZE    = 10          # 隠れ層のユニット数
H3SIZE    = 5           # 隠れ層のユニット数
OSIZE     = 2           # 出力サイズ（redバイナリ、blueバイナリ）
BATCH     = 10          # バッチサイズ
ITER      = 50000       # 学習の回数
LEARNRATE = 0.05        # 学習係数
FRAME     = 30          # アニメーションフレーム数（1層ごと）
```

データサンプル数として、赤色5,000個、青色5,000個の合計10,000個を用意しました。
ニューラルネットによって、入力したデータの座標空間がどのように「歪められていくか」を見せるためにグリッド点も用意します。（元記事では格子線だったのですが、matplotlibで2次元の格子線を描く方法が分からずやむなくグリッド点にしています）

入力値は各データ点のx,y座標です。隠れ層のユニットサイズはそれぞれ、第1層:20、第2層:10、第3層:5と設定しました。出力は入力値が赤である確率と青である確率の2つの値を返すようにしています。

バッチサイズは10で、学習は5万回としました。色々試してみると数千回で十分誤差が下がっているようですが、安全を見て5万回としています。最後の変数は動画作成時のための変数です。


#### データの生成
データ生成にあたっては、元記事の数字から「恐らくこの関数だろう」という値を生成しました。3パターンのデータを生成しています。
1. 二つの三角関数
1. 二つの螺旋関数
1. 円とドーナツ


##### 1. 二つの三角関数
$x \in \[-1, 1\]$に対して、
```math
\begin{align}
y_{red} = - 0.5  \cos(3x) + 0.5 \\
y_{red} = - 0.5  \cos(3x) - 0.5
\end{align}
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/107669/a0ceeb3a-270e-d2ce-09df-72c8c8b1b6f1.png)


##### 2. 二つの螺旋関数
$\theta \in \[-4.0, 0.75\]$に対して、
```math
\begin{align}
x_{red} &=& - 0.65 \times e^{0.4  \theta} \times \cos{\theta} \\
y_{red} &=& - 1.0  \times e^{0.4  \theta} \times \sin{\theta} \\
x_{blue} &=& 0.65 \times e^{0.4  \theta} \times \cos{\theta} \\
y_{blue} &=& 1.0  \times e^{0.4  \theta} \times \sin{\theta} 
\end{align}
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/107669/ad6d2ce1-162a-0fea-e19f-b4dbecd94acd.png)


##### 2. 円内部とドーナツ
```math
\begin{align}
円 &=& \lbrace x \mid d(x,0)<1/3 \rbrace  \\
ドーナツ &=& \lbrace x \mid 2/3<d(x,0)<1 \rbrace 
\end{align}
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/107669/d4f2d13f-a183-b53c-905c-7eb39306dbeb.png)


それぞれを生成するコードはこちらになります。
```python
def dataCreate(id):
  # データ構造：[x, y, red, blue]
  # red列は赤なら1,青なら0

  # 二つの三角関数
  if id == 1:
    x1 = np.linspace(-1.0, 1.0, N)
    x2 = x1
    y1 = - np.cos(x1*3.0) / 2.0  + 0.5
    y2 = - np.cos(x2*3.0) / 2.0  - 0.5
  # 螺旋関数
  elif id == 2:
    theta = np.linspace(-4.0, 0.75, N)
    a     = 0.5
    b     = 0.4

    x1 = - 1.3 * a * np.exp(b * theta) * np.cos(theta)
    y1 = - 2.0 * a * np.exp(b * theta) * np.sin(theta)
    x2 =   1.3 * a * np.exp(b * theta) * np.cos(theta)
    y2 =   2.0 * a * np.exp(b * theta) * np.sin(theta)
  # 円とドーナツ
  elif id == 3:
    theta1  = 2.0 * np.pi * np.random.rand(N)
    radius1 = np.sqrt(np.random.rand(N)) / 3.0
    x1 = radius1 * np.cos(theta1)
    y1 = radius1 * np.sin(theta1)

    theta2  = 2.0 * np.pi * np.random.rand(N)
    radius2 = np.sqrt(np.random.rand(N)) / 3.0 +  2.0 / 3.0 
    x2 = radius2 * np.cos(theta2)
    y2 = radius2 * np.sin(theta2)

  col1 = np.zeros(N)
  col2 = np.ones(N)

  data1 = np.concatenate([x1.reshape(N,1), y1.reshape(N,1), col1.reshape(N,1), col2.reshape(N,1)], 1)
  data2 = np.concatenate([x2.reshape(N,1), y2.reshape(N,1), col2.reshape(N,1), col1.reshape(N,1)], 1)
  Data  = np.concatenate([data1, data2], 0)

  if id == 3:
    plt.scatter(x1, y1, c='red')
    plt.scatter(x2, y2, c='blue')
  else:
    plt.plot(x1, y1, c='red')
    plt.plot(x2, y2, c='blue')
  plt.xlim(-1.0,1.0)
  plt.ylim(-1.0,1.0)    

  gg     = np.linspace(-1.0, 1.0, GSIZE)
  xx, yy = np.meshgrid(gg,gg)
  grid   = np.concatenate([np.ravel(xx).reshape(GSIZE**2,1), np.ravel(yy).reshape(GSIZE**2,1)], 1)

  return Data, grid
```

#### 動画作成
動画の作成方法ですが、下記ステップになります。
1. まずはサンプルデータを元に学習実施
1. 学習完了後、すべての入力データ（およびグリッド点）を使って、各中間層（線形変換層、活性化層も分ける）毎にノードを生成
1. 各層を2次元空間にマッピング。ノードが3以上ある場合は、一番上の2つのノードの値を利用（2つのノードへと射影する）
1. 層と層の間を滑らかにつなぐために、線形補間
1. 全体を動画化（入力層→第1層線形変換層→活性化層→第2層線形変換層→・・→出力層）

以下が動画作成のコードになります。
```python
def aniCreator(Data, id):
  # 入力層から隠れ層・・・出力層までを可視化
  # ノードが3以上の場合は、はじめの2つのノードに射影
  fig = plt.figure()
  ims = []

  # 学習後パラメータを使って隠れ層を構築
  u1, z1, u2, z2, u3, z3, u4, out  = NN.forward(Data[:,[0,1]])

  # 格子点 ⇒ 隠れ層
  u1G, z1G, u2G, z2G, u3G, z3G, u4G, outG  = NN.forward(grid)

  for i in range(FRAME * 8):
    if i <= FRAME - 1:
      layer1  = Data
      layerG1 = grid
      layer2  = Data
      layerG2 = grid
    elif i <= FRAME*2 - 1:
      layer1  = Data
      layerG1 = grid
      layer2  = u1
      layerG2 = u1G
    elif i <= FRAME*3 - 1:
      layer1  = u1
      layerG1 = u1G
      layer2  = z1
      layerG2 = z1G
    elif i <= FRAME*4 - 1:
      layer1  = z1
      layerG1 = z1G
      layer2  = u2
      layerG2 = u2G
    elif i <= FRAME*5 - 1:
      layer1  = u2
      layerG1 = u2G
      layer2  = z2
      layerG2 = z2G
    elif i <= FRAME*6 - 1:
      layer1  = z2
      layerG1 = z2G
      layer2  = u3
      layerG2 = u3G
    elif i <= FRAME*7 - 1:
      layer1  = u3
      layerG1 = u3G
      layer2  = z3
      layerG2 = z3G
    else:
      layer1  = z3
      layerG1 = z3G
      layer2  = z3
      layerG2 = z3G

    j = i % FRAME
    col1x = (FRAME-j)/FRAME * layer1[np.where(Data[:,2] == 0), 0] + j/FRAME * layer2[np.where(Data[:,2] == 0), 0]
    col1y = (FRAME-j)/FRAME * layer1[np.where(Data[:,2] == 0), 1] + j/FRAME * layer2[np.where(Data[:,2] == 0), 1]
    col2x = (FRAME-j)/FRAME * layer1[np.where(Data[:,2] == 1), 0] + j/FRAME * layer2[np.where(Data[:,2] == 1), 0]
    col2y = (FRAME-j)/FRAME * layer1[np.where(Data[:,2] == 1), 1] + j/FRAME * layer2[np.where(Data[:,2] == 1), 1]

    gridx = (FRAME-j)/FRAME * layerG1[:,0] + j/FRAME * layerG2[:,0]
    gridy = (FRAME-j)/FRAME * layerG1[:,1] + j/FRAME * layerG2[:,1]

    im1 = plt.scatter(gridx, gridy, s=2, c='darkgray', edgecolors='darkgray') 
    im2 = plt.scatter(col1x, col1y, s=5, c='red', edgecolors='red')              
    im3 = plt.scatter(col2x, col2y, s=5, c='blue', edgecolors='blue')             
    ims.append([im1, im2, im3])                  # グラフを配列 ims に追加

  # 描画設定
  plt.xlim(-1.5,1.5)
  plt.ylim(-1.5,1.5)

  # 動画生成
  anim = animation.ArtistAnimation(fig, ims, interval=50)
  if id == 1:
    anim.save("animation1.gif", writer="pillow")
  elif id == 2:
    anim.save("animation2.gif", writer="pillow")
  else:
    anim.save("animation3.gif", writer="pillow")

  # Google Colaboratoryの場合必要
  rc('animation', html='jshtml')
  plt.close()
  return anim
```

各層から次の層へと線形補間する際の分岐が冗長な形になってしまいました。もっと効率的な書き方はできると思いますが思いつきませんでした。


# 結果
#### 1. 二つの三角関数
まずは2つの三角関数を学習してみました。元記事では静止画しかありませんでしたが、折角なのでこちらも動画化しました。

![animation1.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/107669/4c89ae2d-1581-71bd-bed8-76aa084e79a0.gif)

中間層が3層、各層のユニット数が十分あるため簡単に分類できています。出力層に近づくにつれ、赤部分と青部分がこれ以上ないくらい離れていくのが見えます。


#### 2. 二つの螺旋関数
元記事では各層のノード数の最大値が2つの場合とそうでない場合を載せていますが、こちらではノード数が十分ある場合のみ再現しました。再現といっても、パラメータの数がそもそも違うのと、線形分離する「答え」がたくさんあるため、厳密に元記事を再現しているわけではありません。ただ、入力層の空間を捻じ曲げながら赤と青を分離しようとしている点は再現できているかなと思います。

![animation2.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/107669/fc0119e6-7472-3b42-2007-4bb17852bbd2.gif)


#### 3. 円とドーナツ
最後に円とドーナツですが、こちらについても元記事とは異なりここではノード数が3つ以上あるため、2次元では絶対に線形分離できない問題に対して、高次元空間へと射影して分離できています。

![animation3.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/107669/d54da46e-0ad3-2acc-269e-6700e8035138.gif)



# まとめと考察
Colah氏の示唆に富む記事に載っていた面白い動画を再現してみました。一通り再現してみて改めて、「シンプルなデータセットと簡単なモデルに落とし込んだうえでニューラルネットの挙動を可視化するとすごく示唆にとんでいるな」と思いました。

線形変換層は元のデータ空間を平行移動させたり回転させたり、拡大縮小することはできても、点同士の相対的な位置関係は変わりません。そこで活性化関数で空間を捻じ曲げることで位置関係を壊すことができます。また、一旦活性関数で捻じ曲げた空間を再度拡大縮小したり回転させることで更なる捻じ曲げに対する準備ができるようになります。このような「こねくり回し」をたくさん繰り返すことで、複雑に絡み合ったデータを特徴量毎に分類することができるようになりそうだなと思いました。厳密な理解ではないですが、層を増やして「ディープにする」意味もイメージが湧きます。

今後は、活性化関数を変更してみたり、もう少し難易度が高い問題を扱ってみたりしながらニューラルネットワークの挙動を理解したいと思います。