だいたい1時間で分かるDeepLearningの仕組み by Julia

本稿は、Julia Advent Calendar 2019 21日目の記事です

Advent Calendar 初参加ですが、よろしくお願いいたします



本稿の内容および対象読者

本稿では書籍『ゼロから作るDeepLearning〜Pythonで学ぶディープラーニングの理論と実装〜』の内容を噛み砕いて解説するとともに、プログラミング言語 Julia でニューラルネットワークを実装していきます

Juliaを利用する理由は以下の通りです

  • 素の計算速度がPythonより遥かに速い
  • Python(NumPy)より行列式の扱いが直感的で分かりやすい
  • 数式をそのままプログラムに落とし込みやすい
  • Pythonの便利なパッケージを都合よく利用することも可能

なお、本稿の全プログラムは https://github.com/amenoyoya/ml-algorithm に置いてあります

対象読者

各プログラムにはなるべく多くコメントを入れ、読めばある程度内容が分かるように作っているつもりですが、数学の基礎知識や Julia の基本文法について、一から解説するということは行っていません

そのため、対象読者としては以下のような方を想定しています

  • 高校レベルの数学(行列、微分の基礎的な知識)がある程度分かっている人
  • Julia 1.3.0 の基本文法が分かっている人
  • がっつりDeepLearningをやりたい(暇な)人
    • DeepLearningで何か作りたいものがある場合は、PyTorchやKeras等のライブラリの使い方を覚えたほうが手っ取り早い

環境構築

Environment

  • OS:
    • Ubuntu 18.04
    • Windows 10
  • Jupyter Notebook: 4.4.0
    • Anaconda: 4.5.11
      • Python: 3.7.4
  • Julia: 1.3.0

Installation on Ubuntu 18.04

# install in home directory
$ cd ~

# download julia-1.3.0
$ wget -O - https://julialang-s3.julialang.org/bin/linux/x64/1.3/julia-1.3.0-linux-x86_64.tar.gz | tar zxvf -

# export path
## /usr/local/bin/ に julia 実行ファイルのシンボリックを作成
$ sudo ln -s ~/julia-1.3.0/bin/julia /usr/local/bin/julia

# confirm version
$ julia -v
julia version 1.3.0

Installation on Windows 10

※ここでは Chocolatey パッケージマネージャを使った方法を採用します

Win + X |> A キー => 管理者権限 PowerShell 起動

# chocolateyパッケージマネージャを入れていない場合は導入する
> Set-ExecutionPolicy Bypass -Scope Process -Force; iex ((New-Object System.Net.WebClient).DownloadString('https://chocolatey.org/install.ps1'))

# install julia
> choco install -y julia

# confirm version
> julia -v
julia version 1.3.0

Julia でよく使うパッケージを導入

インストールするパッケージは以下の通りです

  • HTTP:
    • HTTP通信を行うためのパッケージ
  • DataFrames:

    • データフレームを扱うためのパッケージ

      $ julia -e 'using Pkg; Pkg.add("HTTP"); Pkg.add("DataFrames")'

Jupyter Notebook 導入

Jupyter Notebook
Pythonプログラムを実行しながらノート形式で保存できるソフトウェアです
プログラムの実行結果を確認しながら分析を進められるため、機械学習やデータ分析においてよく用いられます
Python以外の言語にも多数対応しており、Juliaでも利用することが可能です

※ 本稿では、Anaconda(Python3)環境は導入済みという想定でインストール作業を進めます

# install jupyter notebook
$ conda install jupyter

# Jupyter Notebook 用のJuliaカーネルのインストール
$ julia -e 'using Pkg; Pkg.add("IJulia"); Pkg.build("IJulia")'

# jupyter notebook のカーネルを確認
$ jupyter kernelspec list
Available kernels:
  julia-1.3    /home/user/.local/share/jupyter/kernels/julia-1.3   # <- Juliaが使えるようになっている
  python3      /home/user/miniconda3/share/jupyter/kernels/python3

# Jupyter Notebook 起動
$ jupyter notebook

# => localhost:8888 で Jupyter Notebook 起動
## Juliaを使うには New > Julia 1.3.0 を選択する

機械学習用パッケージ導入

機械学習、DeepLearning用のパッケージを導入します

とは言え、基本的にはPythonのパッケージをラッピングしたものばかりであるため、まずはPythonパッケージのインストールから始めます

# Pythonパッケージ ScikitLearn, Matplotlib をAnacondaでインストール
$ conda install scikit-learn matplotlib

# install Julia PyCall package
## PyCall: PythonプログラムをJuliaから呼び出すパッケージ
$ julia -e 'using Pkg; Pkg.add("PyCall")'

# install Julia ScikitLearn bundled package
$ julia -e 'using Pkg; Pkg.add("ScikitLearn")'

# install Julia Matplotlib bundled package
$ julia -e 'using Pkg; Pkg.add("PyPlot")'

# install Julia MachineLeaning Datasets package
## MLDatasets: MNIST等、機械学習用のデータセットをまとめたパッケージ
$ julia -e 'using Pkg; Pkg.add("MLDatasets")'

パーセプトロン

まず、ディープラーニング(ディープニューラルネットワーク)の起源とされる、パーセプトロンの仕組みを学びます

  • パーセプトロン
    • 1957年ローゼンブラットにより考案されたアルゴリズム
    • 複数の信号を受け取り、一つの信号を出力する
    • 出力する信号は「流す/流さない(1 or 0)」の二値
    • 各入力信号には重要度があり、それを重みで表現する

perceptron.png

Input: Julia 1.3.0
# 内積関数 dot を使えるようにする
using LinearAlgebra

# パーセプトロン
## 信号を受け取り、<入力信号> * <重み> + <バイアス> から出力(0 | 1)を計算するパーセプトロン関数を生成
make_perceptron(weight::Array{Float64,1}, bias::Float64) = begin
    # 入力信号x と 重みweight の内積に biasを加算した値が 0を超えているなら 1を返す
    return (x::Array{Int,1}) -> dot(x, weight) + bias > 0 ? 1 : 0
end

ANDパーセプトロン

パーセプトロンを使って、単純な演算装置を構築することができます

まずは、AND演算をパーセプトロンで実装してみましょう

AND演算とは以下のような演算です

  • 0 か 1 のいずれかの値を引数として、2つ受け取る
  • 2つの引数が両方共 1 の場合 1 を返す(それ以外であれば 0 を返す)
Input: Julia 1.3.0
"""
AND
0, 0 => 0
0, 1 => 0
1, 0 => 0
1, 1 => 1
"""

# ANDパーセプトロン
## 重み: 0.5, 0.5
## バイアス: -0.75
AND = make_perceptron([0.5, 0.5], -0.75)

for x in [[0, 0], [0, 1], [1, 0], [1, 1]]
    println("$(x) => $(AND(x))")
end
Output
[0, 0] => 0
[0, 1] => 0
[1, 0] => 0
[1, 1] => 1

パーセプトロンの数式の意味

パーセプトロンにおいて、$w_i$ は入力信号の重要度、$b$ はニューロンの発火のしやすさを表します

また、基本的に $\sum_{i=1}^n |w_i| = 1$ になり、$|b| \leqq 1$ となります

通常 $b < 0$ なので、$b$ の絶対値はニューロンの発火のしづらさを表します

この前提のもと、ANDパーセプトロンのパラメータの意味を考えてみましょう

  • $w_1, w_2$ ともに 0.5 であるため、2つの入力信号の重要度は変わらない
    • 言い換えれば、2つの入力値の入力順を変えても結果は変わらない
  • $b$ は -0.75 であり、$ 0.25 - 1 $ に等しい
    • つまり、25%の確率で発火し、75%の確率で発火しないことを意味する

次に、実際のAND演算子を考えてみます

  • 2つの入力値は、入力順が変わっても出力値は変わらない
    • $[0, 1] = [1, 0] = 0$
  • 入力値の4つの組み合わせ($[0, 0], [0, 1], [1, 0], [1, 1]$)のうち、出力値が1(発火)になるのは1つ($[1, 1]$)のみ
    • 発火の確率は $14 = 0.25$ だから、25%

このように、入力信号の重要度と発火のしやすさだけを考えれば、ビット演算子のニューロンを組み立てることが可能です

Input: Julia 1.3.0
"""
OR
0, 0 => 0
0, 1 => 1
1, 0 => 1
1, 1 => 1
"""

# ORパーセプトロン
## 重み: 0.5, 0.5 => 入力信号の重要度は変わらない
## バイアス: -0.25 => 25%の確率で発火しない([0,0]のみ), 75%の確率で発火([0,1], [1,0], [1,1])
OR = make_perceptron([0.5, 0.5], -0.25)

for x in [[0, 0], [0, 1], [1, 0], [1, 1]]
    println("$(x) => $(OR(x))")
end
Output
[0, 0] => 0
[0, 1] => 1
[1, 0] => 1
[1, 1] => 1
Input: Julia 1.3.0
"""
NAND: NOT AND
0, 0 => 1
0, 1 => 1
1, 0 => 1
1, 1 => 0
"""

# NANDパーセプトロン
## 重みとバイアスをそれぞれ符号反転すれば NOT 演算に等しくなる: NOT AND = -AND(x, b)
## 重み: -0.5, -0.5 => 入力信号の重要度は変わらない
## バイアス: 0.75 => 25%の確率で発火([1,1])
NAND = make_perceptron([-0.5, -0.5], 0.75)

for x in [[0, 0], [0, 1], [1, 0], [1, 1]]
    println("$(x) => $(NAND(x))")
end
Output
[0, 0] => 1
[0, 1] => 1
[1, 0] => 1
[1, 1] => 0
Input: Julia 1.3.0
"""
NOR: NOT OR
0, 0 => 1
0, 1 => 0
1, 0 => 0
1, 1 => 0
"""

# NORパーセプトロン
## 重みとバイアスをそれぞれ符号反転すれば NOT 演算に等しくなる: NOT OR = -OR(x, b)
## 重み: -0.5, -0.5 => 入力信号の重要度は変わらない
## バイアス: 0.25 => 75%の確率で発火([0,1], [1,0], [1,1])
NOR = make_perceptron([-0.5, -0.5], 0.25)

for x in [[0, 0], [0, 1], [1, 0], [1, 1]]
    println("$(x) => $(NOR(x))")
end
Output
[0, 0] => 1
[0, 1] => 0
[1, 0] => 0
[1, 1] => 0

分類問題としてのパーセプトロン

上記のように、パーセプトロンは極めて簡単な式で分類問題を解くことができます

すなわちパーセプトロンは「入力信号として2つのビット値をとり、その出力値を 0 or 1 に分類するアルゴリズム」であると言えます

図示すると、以下のような直線で分類を行っていると考えることができます

perceptron_and_or.png

パーセプトロンの限界

パーセプトロンによる分類は、直線を用いたものであるため、直線で区切ることができないデータを分類することはできません

ビット演算で言うと、XORなどがそれに当てはまります(下図のように、XORの出力値は直線で分類することができない)

perceptron_xor.png


多層パーセプトロン

XOR回路図

XORは、正確に言えば「単一のパーセプトロンで分類することができない」だけであり、複数のパーセプトロンを重ねればこの問題を解くことができます

すなわちXORの分類問題は、多層パーセプトロンで解くことが可能です

そのためには、まずはXORゲートをAND, OR, NAND, NORゲートの組み合わせに分解する必要があります(下図)

perceptron_xor.png

Input: Julia 1.3.0
"""
XOR: AND(NAND(x1,x2), OR(x1,x2))
0, 0 => 0
0, 1 => 1
1, 0 => 1
1, 1 => 0
"""

# XORパーセプトロン
## NAND, ORパーセプトロンを ANDパーセプトロンで結合する
XOR(x::Array{Int,1}) = AND([NAND(x), OR(x)])

for x in [[0, 0], [0, 1], [1, 0], [1, 1]]
    println("$(x) => $(XOR(x))")
end
Output
[0, 0] => 0
[0, 1] => 1
[1, 0] => 1
[1, 1] => 0

このように、パーセプトロンの層を積み重ねることでより複雑な分類問題を解くことが可能です

この「層を積み重ねる」という考えは、ディープラーニングの基本思想の1つであるため、極めて重要です

今回のXOR分類問題は、中間層を一つ追加するだけで解くことができたが、より多くの中間層を積み重ねればより複雑な問題にも応用可能となります

perceptron2.png


ニューラルネットワーク

それではいよいよ、ディープラーニングの基本アルゴリズムであるニューラルネットワークについて学びます

パーセプトロンからニューラルネットワークへ

ここで、パーセプトロンの仕組みを再び考えてみましょう

perceptron0.png

ここでバイアスについて、重みが $b$ で値が 1 固定の入力信号であると考えると、以下のように変形できます

perceptron_bias.png

これによって、パーセプトロンの式を $入力値×重み$ に統一することができます($a = \sum_{i=1}^n x_i*w_i$)

さらに、出力値 $y$ を得るための処理として、パーセプトロンでは $a \leqq 0 → 0$ | $a > 0 → 1$ という条件分岐を行っています

このように、入力信号の総和を出力信号に変換するための処理関数を活性化関数と呼び、$h(a)$ で表します

つまり、パーセプトロンにおける活性化関数は h(a) = 0 if a <= 0 else 1 と表現できます

このように、特定条件で0から1へ突然出力値が変わるような活性化関数をステップ関数と呼びます

perceptron_h.png

シグモイド関数

ステップ関数は微分することができないため、ニューラルネットワークにおいては、ステップ関数の代わりにシグモイド関数がよく使われてきました

$$ シグモイド関数: h(a) = \frac{1}{1 + e^{-a}} $$

"""

"""

# ステップ関数
step(a::Float64)::Float64 = (a <= 0 ? 0 : 1)

# シグモイド関数
sigmoid(a::Float64)::Float64 = 1 / (1 + exp(-a))

# プロット用|横軸
## collect(start:interval:end): start〜endまでinterval刻みの数値配列を生成
a = collect(-5.0:0.1:5.0)

# ステップ関数実行
## function.(array) で function(v) for v in array を実行
y = step.(a)
Input: Julia 1.3.0
# PyPlotパッケージ利用
using PyPlot

# ステップ関数プロット
plot(a, y)

step-graph.png

Input: Julia 1.3.0
# シグモイド関数プロット
plot(a, sigmoid.(a))

sigmoid-graph.png

ステップ関数とシグモイド関数

ステップ関数とシグモイド関数は、上記グラフを見ても分かるように、滑らかさ(連続性)に違いがあります

一方で、以下のような共通点もあります

  • 出力値は $0 \leqq y \leqq 1$ となる
    • 入力値が小さいほど 0 に近く、入力値が大きいほど 1 に近くなる
  • 非線形関数である
    • 一時直線では表現できない関数である

非線形関数

ニューラルネットワークでは、活性化関数に非線形関数を用いる必要があります

逆に言えば、活性化関数に線形関数を用いてはなりません

これは、線形関数を用いてしまうと、どれだけ層を深くしても、それと同じことを行う「隠れ層(中間層)のないネットワーク」が存在することになってしまうためです

例えば、活性化関数として $h(x) = cx$ という線形関数を採用し、$y(x) = h(h(h(x)))$ という3層ネットワークを構築したとしましょう

このネットワークは $y(x) = c * c * c * x$ という乗算処理を行いますが、$c$ は任意定数であることから、$y(x) = ax$($a = c^3$)という1層ネットワークに等しくなってしまいます

このように、線形関数を活性化関数として採用してしまうと、層を深くする意味がなくなってしまうため、必ず非線形関数を活性化関数とする必要があります

ReLU関数

シグモイド関数は古くから活性化関数としてよく使われてきましたが、最近ではReLU(Rectified Linear Unit)関数が採用されることが増えてきました

ReLU関数は、入力値が0以下なら0を出力し、入力値が0を超えたらその値をそのまま出力します

そのため、出力値が収束することはなく、より複雑な問題を解くことが可能です

Input: Julia 1.3.0
# ReLU関数
ReLU(x::Float64)::Float64 = (x <= 0 ? 0 : x)

# プロット用|入力値: -5.0, -4.9, ..., 4.9, 5.0
x = collect(-5.0:0.1:5.0)

# ReLU関数プロット
plot(x, ReLU.(x))

relu-graph.png


3層ニューラルネットワーク

パーセプトロン同様、ニューラルネットワークも複数の層を積み重ねることで複雑な計算を行うことができます

このように複数の層を重ねることで、「深く」計算できるようにしたニューラルネットワークをディープニューラルネットワークと呼び、これが俗にディープラーニングと呼ばれている技術です

3層ニューラルネットワークの実装

以下のような構造の3層ニューラルネットワークを実装してみましょう

neural_network_3.png

重みの記述を $w_{n p}^{(i)}$(n: 次層のニューロンの番号, p: 前層のニューロンの番号, i: 層の番号)とし、バイアスを加えると以下のようなネットワーク構造となります

neural_network_bias.png

さらに、各中間層(隠れ層)に活性化関数 $h()$ を導入し、層同士の信号伝達を表現すると以下のようになります

neural_network_layer.png

最後に、出力層に活性化関数 $σ()$ を導入すると、ニューラルネットワークの概念図が完成します

neural_network_y.png

ニューラルネットワーク1層のテンソル計算

$$ \begin{array}{ll} \begin{bmatrix} a_1^{(1)} \\ a_2^{(1)} \\ a_3^{(1)} \end{bmatrix} &= \begin{bmatrix} x_1 & x_2 \end{bmatrix} \begin{bmatrix} w_{11}^{(1)} & w_{21}^{(1)} & w_{31}^{(1)} \\ w_{12}^{(1)} & w_{22}^{(1)} & w_{32}^{(1)} \end{bmatrix} + \begin{bmatrix} b_1^{(1)} \\ b_2^{(1)} \\ b_3^{(1)} \end{bmatrix} \\ &= \begin{bmatrix} \sum_{p=1}^{2}{(x_1 * w_{1 p}^{(1)})} + b_1^{(1)} \\ \sum_{p=1}^{2}{(x_2 * w_{2 p}^{(1)})} + b_2^{(1)} \\ \sum_{p=1}^{2}{(x_3 * w_{3 p}^{(1)})} + b_3^{(1)} \end{bmatrix} \end{array} $$
Input: Julia 1.3.0
"""
3
"""

# ニューラルネットワーク構造体
mutable struct Network
    l::Int # ネットワーク層数
    b::Vector{Array{Float64,2}} # バイアス: [層1: Array{1x次層ニューロン数}, 層2: Array{1x次層ニューロン数}, ...]
    w::Vector{Array{Float64,2}} # 重み: [層1: Array{前層ニューロン数x次層ニューロン数}, 層2: Array{前層ニューロン数x次層ニューロン数}, ...]
end

# ニューラルネットワーク1層計算関数
## Network構造体, 層番号, 活性化関数, 層の入力信号 -> 層の出力信号 z
forward(network::Network, i::Int, h, x::Array{Float64,2})::Array{Float64,2} =
    h.(x * network.w[i] + network.b[i])

# ニューラルネットワーク計算関数
## Network構造体, 層番号, 中間層の活性化関数, 出力層の活性化関数, 入力信号 -> 出力信号 y
execute(network::Network, h, σ, x::Array{Float64,2}, i::Int=1)::Array{Float64,2} = 
    i == network.l ? forward(network, i, σ, x) : execute(network, h, σ, forward(network, i, h, x), i + 1)

# 中間層の活性化関数
## シグモイド関数
# シグモイド関数
sigmoid(a::Float64)::Float64 = 1 / (1 + exp(-a))

# 出力層の活性化関数
## 恒等関数: 入力値をそのまま出力
identity(x::Float64)::Float64 = x

# 3層ニューラルネットワーク構築
network = Network(
    3,
    [
        # i番目の層のバイアス: 次層の数×1 Array{Float64,2}
        [0.1 0.2 0.3], # 1 x 3
        [0.1 0.2], # 1 x 2
        [0.1 0.2] # 1 x 2
    ],
    [
        # i番目の層の重み: 前層の数×次層の数 Array{Float64,2}
        [0.1 0.3 0.5; 0.2 0.4 0.6], # 2 x 3
        [0.1 0.4; 0.2 0.5; 0.3 0.6], # 3 x 2
        [0.1 0.3; 0.2 0.4] # 2 x 2
    ]
)

# ニューラルネットワーク計算実行
## 中間層の活性化関数: sigmoid, 出力層の活性化関数: identity
## 入力値: x1 = 1.0, x2 = 0.5
execute(network, sigmoid, identity, [1.0 0.5])
Output
1×2 Array{Float64,2}:
 0.316827  0.696279

出力層の設計

ニューラルネットワークは、分類問題回帰問題の両方に用いることができます

ただし、分類問題と回帰問題のどちらに用いるかで、出力層の活性化関数を変更する必要があります

  • 分類問題
    • データがどのクラスに属するかを判定する問題
    • 活性化関数には、ソフトマックス関数のような、入力信号の総和を引数にとる関数がよく使われる
  • 回帰問題
    • 入力されたデータから数値を予測する問題
    • 活性化関数には、恒等関数のような、入力信号をそのまま引数にとる関数がよく使われる

neural_network_activate.png

ソフトマックス関数

分類問題の活性化関数として用いられるソフトマックス関数は以下のような式で表現されます

$$ y_k = \frac{exp(a_k)}{\sum_i^n exp(a_i)} $$

Input: Julia 1.3.0
# ソフトマックス関数
softmax(a::Vector{Float64})::Vector{Float64} = exp.(a) / sum(exp.(a))

# 動作確認
softmax([0.3, 2.9, 4.0])
Output
3-element Array{Float64,1}:
 0.018211273295547534
 0.2451918129350739  
 0.7365969137693785
オーバフロー対策

ソフトマックス関数は指数関数の計算を行うため、コンピュータ上で実装する場合、オーバフローに注意する必要があります

例えばJuliaの場合、$e^{1000}$ の結果は無限大を表す Inf となってしまいます

そのため、このような値を除算に用いると数値が不安定になってしまいます

このようなオーバフローを防ぐために、ソフトマックス関数を変形しましょう

任意定数 $C$ を分子と分母の両方に掛けると、以下のように変形できます

$$ \begin{array}{ll} y_k &=& \frac{exp(a_k)}{\sum_{i=1}^n exp(a_i)} \\ &=& \frac{C*exp(a_k)}{C*\sum_{i=1}^n exp(a_i)} \\ &=& \frac{exp(a_k + log(C))}{\sum_{i=1}^n exp(a_i + log(C))} \\ &=& \frac{exp(a_k + C')}{\sum_{i=1}^n exp(a_i + C')} \end{array} $$

上記より、$a_k$ から任意の定数を減算しても計算結果は変わりません

従って、入力信号のうち最大の値を各入力信号から減算すればオーバフローを回避することが可能です

Input: Julia 1.3.0
# オーバフロー対策済みソフトマックス関数
softmax(a::Vector{Float64})::Vector{Float64} = begin
    c = maximum(a)
    exp_a = exp.(a .- c)
    y = exp_a / sum(exp_a)
end

# 動作確認
y = softmax([0.3, 2.9, 4.0])
println(y)

# ソフトマックス関数は出力値の合計が 1 になる
println(sum(y))
Output
[0.01821127329554753, 0.24519181293507386, 0.7365969137693786]
1.0
ソフトマックス関数の特徴

ソフトマックス関数の出力値は、0〜1 の間の実数となり、その合計値は 1 になります

従って、ソフトマックス関数の出力値は、分類問題における各クラスに属する確率を表していると解釈することができます

ここで、ソフトマックス関数の出力値は入力値の大きさに単調比例する($exp(a_k)$ が単調増加関数であるため)ため、大小関係は関数処理前後で変わりません

そのため、一般にニューラルネットワークによる推論を行う場合、活性化関数(ソフトマックス関数)は省略されることが多いです

出力層のニューロン数

分類問題において、出力層のニューロン数は分類するクラス数に等しくなります

例えば、0〜9 の手書き数字を10個の数値クラス(0〜9)に分類する場合、出力層のニューロン数は10個必要になります

また、基本的に出力層のニューロンが出力する値が、そのままそのクラスに属する確率であると解釈されます


手書き数字の認識

続いて、学習済みのニューラルネットワークモデルを用いて、手書き数字の分類を行いましょう

この推論処理は、ニューラルネットワークの順方向伝搬と呼びます

MNISTデータセット

ここでは、MNISTという手書き数字の画像データセットを使用します

MNISTデータセットは0〜9までの数字画像から構成されており、訓練用画像60,000枚、テスト用画像10,000枚から成ります

一般的なMNISTデータセットの使い方では、訓練用画像を用いて学習を行い、学習したモデルでテスト画像に対してどれだけ正しく分類できるかを計測します

Input: Julia 1.3.0
# MLDatasetsパッケージのMNISTSデータセットを使う
using MLDatasets

# 訓練用画像データと教師データをロード
## train_x: 特徴量=<画像データ|28x28 グレースケール画像 60,000枚>{28x28x60000 Array{UInt8, 3}}
## train_y: 目的変数=<数値クラス|[0..9]の数値 60,000個>{60000 Array{Int, 1}}
train_x, train_y = MNIST.traindata()

# 一枚目の画像データ(28x28)をプロット
## imshowは縦軸と横軸が反転しているため転置行列を渡す
imshow(train_x[:, :, 1]')

# 一枚目の画像に描かれている数字
train_y[1]

5.png

ニューラルネットワーク構築

今回はすでに学習済みのモデルを使って推論処理のみを行ってみるため、そのモデルに合わせたニューラルネットワークを構築します

構築するニューラルネットワークは以下のようなものです

  • 3層ニューラルネットワーク
    • 入力層:
      • 手書き数字の画像データ(サイズ: 28x28)
      • ニューロン数: 28 * 28 = 784
      • 各ニューロンの入力値は 0.0〜1.0 の実数型である必要がある
    • 中間層(隠れ層)1:
      • ニューロン数: 50
      • 活性化関数: シグモイド関数
    • 中間層(隠れ層)2:
      • ニューロン数: 100
      • 活性化関数: シグモイド関数
    • 出力層:
      • ニューロン数: 10(0〜9の数字クラスに分類するため)
      • 活性化関数: ソフトマックス関数

なお、学習済みモデルは python-pickle 形式で保存されているものをGitHubからダウンロードして使用します

Input: Julia 1.3.0
"""
Python picklePython
"""

using HTTP
using PyCall
pickle = pyimport("pickle")

# ファイルにオブジェクト状態保存
save_pickle(filename, obj) = begin
    out = open(filename, "w")
    pickle.dump(obj, out)
    close(out)
end

# ファイルからオブジェクト状態復元
load_pickle(filename) = begin
    @pywith pybuiltin("open")(filename, "rb") as f begin
        return pickle.load(f)
    end
end

# インターネット上のファイルからオブジェクト状態復元
model = pickle.loads(
    HTTP.request(
        "GET", "https://github.com/oreilly-japan/deep-learning-from-scratch/blob/master/ch03/sample_weight.pkl?raw=true"
    ).body
)

"""

"""

# ニューラルネットワーク構造体
mutable struct Network
    l::Int # ネットワーク層数
    b::Vector{Array{Float64,2}} # バイアス: [層1: Array{1x次層ニューロン数}, 層2: Array{1x次層ニューロン数}, ...]
    w::Vector{Array{Float64,2}} # 重み: [層1: Array{前層ニューロン数x次層ニューロン数}, 層2: Array{前層ニューロン数x次層ニューロン数}, ...]
end

# 順方向伝搬関数
## Network構造体, 層番号, 活性化関数, 層の入力信号 -> 層の出力信号 z
forward(network::Network, i::Int, h, x::Array{Float64,2})::Array{Float64,2} =
    h(x * network.w[i] + network.b[i])

# 推論処理
## Network構造体, 層番号, 中間層の活性化関数, 出力層の活性化関数, 入力信号 -> 出力信号 y
predict(network::Network, h, σ, x::Array{Float64,2}, i::Int=1)::Array{Float64,2} = 
    i == network.l ? forward(network, i, σ, x) : predict(network, h, σ, forward(network, i, h, x), i + 1)

# 学習済みモデルのパラメータをセットして3層ニューラルネットワーク構築
init(model::Dict{Any,Any})::Network = Network(
    3,
    [
        # i番目の層のバイアス: Array{Float32,1} -> 1 x n Array{Float64,2}にキャスト
        Array{Float64,2}(reshape(model["b1"], (1, :))), # 1 x 50 Array{Float64,2}
        Array{Float64,2}(reshape(model["b2"], (1, :))), # 1 x 100 Array{Float64,2}
        Array{Float64,2}(reshape(model["b3"], (1, :))), # 1 x 10 Array{Float64,2}
    ],
    [
        # i番目の層の重み: Array{Float32,2} -> Array{Float64,2}にキャスト
        Array{Float64,2}(model["W1"]),
        Array{Float64,2}(model["W2"]),
        Array{Float64,2}(model["W3"]),
    ]
)

network = init(model)

# テスト用画像データとその解答をロード
test_x, test_y = MNIST.testdata()

# 入力値をニューラルネットワークモデル用に前処理
## 入力層は 1 x 784 Array{Float64,2} を受け付けるため、以下のように変換
x = Array{Float64,3}(reshape(test_x, (1, 784, :))) # => 1 x 784 x 10,000枚 Array{Flaot64,2}

# シグモイド関数
sigmoid(a::Array{Float64,2})::Array{Float64,2} = 1 ./ (1 .+ exp.(-a))

# ソフトマックス関数
softmax(a::Array{Float64,2})::Array{Float64,2} = begin
    c = maximum(a)
    exp_a = exp.(a .- c) # オーバフロー対策
    y = exp_a ./ sum(exp_a)
end

# 学習済みパラメータを適用したニューラルネットワークでテスト画像データ10,000枚に対して推論実行
## 中間層の活性化関数: シグモイド関数
## 出力層の活性化関数: ソフトマックス関数
predict_y = [predict(network, sigmoid, softmax, x[:, :, i]) for i in 1:10000]

# 各予測値のうち、最大値(最も確率の高い)のindex(数値クラス)を取得
## Juliaのindexは1から始まるため -1 する
predicted = [argmax(y)[2]-1 for y in predict_y]


# 解答データと比較し、予測値の正確性を計測
length(predicted[predicted .== test_y]) / 10000
Output
0.9352

このように、極めて単純な3層ニューラルネットワークでも、十分にパラメータを最適化されていれば、93.5%程度の精度を出すことができます

バッチ処理

一つのデータをニューラルネットワークに処理させる場合、1×入力層ニューロン数 の2次元テンソル(行列)を引数として渡していました

しかし通常、機械学習においては複数の(大量の)データを処理するため、今までのようにループ計算するのは余り効率が良くありません

これを解決するために、複数のデータを一つの塊(バッチ)として計算を行う処理をバッチ処理と呼びます

すなわち、1×入力層ニューロン数 の行列をデータ数分ループ計算させるのではなく、データ数×入力層ニューロン数 の行列を計算させる方が効率が良いということです

このアルゴリズムを実装してみましょう

# 入力値をニューラルネットワークモデル用に前処理したもの
# => 1 x 784 x 10,000 Array{Float64,3}
# これを 10,000 x 784 Array{Float64,2} に変換する

# まず、1 x 784 x 10,000 の次元の並び順を 10,000 x 784 x 1 に変換
## [1次元, 2次元, 3次元] => [3次元, 2次元, 1次元]
X = permutedims(x, [3, 2, 1])

# 10,000 x 784 x 1 Array{Float64,3} を 10,000 x 784 Array{Float64,2} に変換
X = Array{Float64,2}(reshape(X, (10000, 784)))

ここで順方向伝搬処理の流れをもう一度整理しましょう

データ1つに対して処理する場合は以下のような流れになっていました

入力層→中間層1
入力信号|(1×入力層ニューロン数)行列 * 重み|(入力層ニューロン数×中間層1ニューロン数)行列 = (1×中間層1ニューロン数)行列
(1×中間層1ニューロン数)行列 + バイアス|(1×中間層1ニューロン数)行列 = (1×中間層1ニューロン数)行列

これをバッチ処理させる場合、以下のようにバイアスの行数と不一致が起こるため、計算できません

入力層→中間層1
入力信号|(10,000×入力層ニューロン数)行列 * 重み|(入力層ニューロン数×中間層1ニューロン数)行列 = (10,000×中間層1ニューロン数)行列
(10,000×中間層1ニューロン数)行列 + バイアス|(1×中間層1ニューロン数)行列 = 行数不一致のため計算不可(10,000 \neq 1)

そのため、forward関数を修正して実行する必要があります

Input: Julia 1.3.0
"""

"""
# ニューラルネットワーク構造体
mutable struct Network
    l::Int # ネットワーク層数
    b::Vector{Array{Float64,2}} # バイアス: [層1: Array{1x次層ニューロン数}, 層2: Array{1x次層ニューロン数}, ...]
    w::Vector{Array{Float64,2}} # 重み: [層1: Array{前層ニューロン数x次層ニューロン数}, 層2: Array{前層ニューロン数x次層ニューロン数}, ...]
end

# 順方向伝搬関数
## Network構造体, 層番号, 活性化関数, 層の入力信号 -> 層の出力信号 z
forward(network::Network, i::Int, h, x::Array{Float64,2})::Array{Float64,2} =
    # バッチ処理の場合、前層の行数が1以上になるため、バイアスの加算をループ加算(.+)に変更
    h(x * network.w[i] .+ network.b[i])

# 推論処理
## Network構造体, 層番号, 中間層の活性化関数, 出力層の活性化関数, 入力信号 -> 出力信号 y
predict(network::Network, h, σ, x::Array{Float64,2}, i::Int=1)::Array{Float64,2} = 
    i == network.l ? forward(network, i, σ, x) : predict(network, h, σ, forward(network, i, h, x), i + 1)

# 学習済みモデルのパラメータをセットして3層ニューラルネットワーク構築
init(model::Dict{Any,Any})::Network = Network(
    3,
    [
        # i番目の層のバイアス: Array{Float32,1} -> 1 x n Array{Float64,2}にキャスト
        Array{Float64,2}(reshape(model["b1"], (1, :))), # 1 x 50 Array{Float64,2}
        Array{Float64,2}(reshape(model["b2"], (1, :))), # 1 x 100 Array{Float64,2}
        Array{Float64,2}(reshape(model["b3"], (1, :))), # 1 x 10 Array{Float64,2}
    ],
    [
        # i番目の層の重み: Array{Float32,2} -> Array{Float64,2}にキャスト
        Array{Float64,2}(model["W1"]),
        Array{Float64,2}(model["W2"]),
        Array{Float64,2}(model["W3"]),
    ]
)

network = init(model)

"""

"""
# バッチ対応版シグモイド関数
sigmoid(a::Array{Float64,2})::Array{Float64,2} = hcat(
    [1 ./ (1 .+ exp.(-a[row, :])) for row in 1:size(a, 1)]...
)'


# バッチ対応版ソフトマックス関数
softmax(a::Array{Float64,2})::Array{Float64,2} = begin
    y = []
    for row in 1:size(a, 1)
        c = maximum(a[row, :])
        exp_a = exp.(a[row, :] .- c) # オーバフロー対策
        push!(y, exp_a ./ sum(exp_a))
    end
    hcat(y...)'
end

"""

"""
# データ数×入力層ニューロン数 の入力信号を一回の順伝搬処理で推論できる
## [predict(network, sigmoid, softmax, x[:, :, i]) for i in 1:10000] より効率よく処理できる
predict_y = predict(network, sigmoid, softmax, X)

# 各行における10クラスのうち、最大値のindexを取得
## argmaxの dims=1 なら各列における行内の比較、dims=2 なら各行における列内の比較
indexes = argmax(predict_y; dims=2)

# CartesianIndexの2次元目の値のみ取得する
indexes = getindex.(indexes, [2])


# Juliaのindexは1から始まるため、-1 して数値クラスを取得
predicted = indexes .- 1

# 解答データと比較して、予測精度を計算
length(predicted[predicted .== test_y]) / 10000
Output
0.9352

備考

実際のところ、Python(NumPy)はforループ処理が極端に遅いため行列計算にすることで効率化されていましたが、Juliaのforループは十分に高速なため、ドット演算子(ベクトル演算子)とどちらが速いのか、筆者はまだ検証できていません

しかし、ベクトル演算子を使うことで無駄な一時変数のアロケーションを減らすことができ、結果的に計算速度を向上させることができるという話もあるため、極力ベクトル計算(テンソル計算)を行っておくに越したことはないと考えられます


データ駆動

機械学習においてはデータが命です

データからパターンを見つけ出し、データから答えを導く出すのが機械学習の本質とも言えます

ここで、入力データからデータの本質的なパターンを導き出す変換器を特徴量と呼びます

特徴量の抽出の歴史

最も単純な分類である線形分類においては、特徴量を抽出する必要はなく、入力データから有限回の学習により問題を解くことが可能でした

これは最初に行ったパーセプトロンによるAND, OR等の演算子の実装において証明されたように「パーセプトロンの収束原理」として知られています

一方、非線形分類問題は入力データからの自動的な学習が不可能であり、特徴量を抽出する必要が出てきました

当初は、特徴量の抽出からデータの学習(アルゴリズムの実装)まで全て人の手で行われており、このモデルはエキスパートシステムと呼ばれています

その後、特徴量はデータサイエンティストと呼ばれる人々によって抽出されるようになり、抽出された特徴量を機械学習アルゴリズム(SVM, kNNなど)を用いて自動的に学習することが可能になりました

そして近年話題となったディープニューラルネットワーク(ディープラーニング)においては、特徴量の抽出からデータの学習までを全てコンピュータによって行われるようになっています(画期的😊)

machine_learning_history.png

パラメータの最適化と損失関数

機械学習およびニューラルネットワークにおいては、データの学習とはすなわちパラメータの最適化です

例えばニューラルネットワークは、重みとバイアスというパラメータの値を少しずつ変化させ、計算結果が教師データと近い値になるように学習を行っています

この時、計算結果と教師データとの誤差を表現する関数を損失関数と呼び、この損失関数を最小化することが学習の目的となります

関数の最小値を求める場合、数学的には微分という手法が用いられます

微分値とは「ある瞬間における変化の量」と定義される値のことであり、わかりやすく言えば、グラフの接線の傾きのことです

簡単のため二次関数を例にすると下図の通り、微分値が0になる(接線の傾きが0になる)点が関数の最小値点であり、パラメータの最適点であると考えることができます

loss_function_optimization.png

ニューラルネットワークの学習においてはこのように、「損失関数」と「微分」がキーワードとなります

微分の実装

微分の数学的な定義は以下のようになります

$$ \frac{df(x)}{dx} = \lim_{h \to 0}{\frac{f(x+h) - f(x)}{h}} $$

この式は、関数 $f(x)$ において $x$ が極めて微小な量 $h$ だけ変化したときの $f(x)$ の変化量を表しています

ここで $x$ の極めて微小な変化量 $h$ をプログラム的にどう表すかが問題となります

数値微分

上記問題の最も単純な解決策は、$h$ に $10^{-50}$ などの極めて小さな値を入れてしまうことです

このような微分手法を数値微分と呼びます

この実装における注意点として、実装するプログラミング言語の丸め誤差を考慮しなければならないことと、数値微分で生じる誤差を極力減らすように工夫しなければならないことが挙げられます

丸め誤差とは、小数の小さな範囲で数値が省略されることで発生する誤差のことで、計算機イプシロンとして定義されています

ちなみに Julia においては 2.220446049250313e-16 という値になっています

理論上は、上記の値より大きな値を $h$ に設定すればよいが、計算量との兼ね合いから、一般的に $10^{-4}$ に設定すれば大抵の場合で上手く行くとされています

続いて、数値微分によって生じる誤差について考えてみましょう

数値微分では、「限りなく0に近い値」を「極めて微小な値」と置き換えて実装するため、その分の誤差が生じます

この誤差を極力減らすため、数値微分の実装では中心差分という差分をとることが多いです

ここで、微分の定義における $f(x+h) - f(x)$ のような差分を前方差分と呼びますが、中心差分では $x$ を中心としてその前後の差分 $f(x+h) - f(x-h)$ をとります

このように差分の平均値を計算することで、誤差を減らすことができます

ここまでの注意点を考慮して数値微分を実装すると以下のようになります

Input: Julia 1.3.0
# 数値微分関数
numeric_diff(f, x) = begin
    h = 1e-4 # 10^(-4)
    (f(x + h) - f(x - h)) / 2h # 中心差分から微分値を計算
end

"""
@test: 

: f(x) = 0.01x^2 + 0.1x
: f'(x) = 0.02x + 0.1
    => f'(5) = 0.2, f'(10) = 0.3
"""
# 例題の設定
f(x) = 0.01x^2 + 0.1x

# 数値微分: f'(5)
## => 0.1999999999990898 ≒ 0.2
println(numeric_diff(f, 5))

# 数値微分: f'(10)
## => 0.2999999999986347 ≒ 0.3
println(numeric_diff(f, 10))
Output
0.1999999999990898
0.2999999999986347

偏微分

上記の微分処理では、変数は $x$ のみだったため話は簡単でした

しかし通常、ニューラルネットワークにおいてはパラメータが複数あるため、1変数のみの微分では対応が難しいです

そこで、偏微分という手法が必要になってきます

ここでは簡単のため、2変数 $x_0, x_1$ で表される以下のような関数を考えてみましょう

$$ f(x_0, x_1) = x_0^2 + x_1^2 $$

Input: Julia 1.3.0
# 2変数関数
f(x0::Float64, x1::Float64) = x0^2 + x1^2

"""
@show: 
"""
# 軸準備
x0 = collect(-3.0:0.1:3.0) # x_0: -3.0, -2.9, -2.8, ..., 2.9, 3.0
x1 = collect(-3.0:0.1:3.0) # x_1: -3.0, -2.9, -2.8, ..., 2.9, 3.0

# X軸: [x0 × x0]: 61 x 61 Array{Float64,2}
X = hcat(fill(x0, length(x0))...)

# Y軸: [x1 × x1] の転置行列: 61 x 61 Array{Float64,2}
Y = hcat(fill(x1, length(x1))...)
Y = Y'

# Z軸: grid状にf(x0, x1)計算 = [f.(x0, x1) × f.(x0, x1')]
## => 61 x 61 Array{Float64, 2}
Z = [f(x, y) for x in x0, y in x1]

# グラフ描画用に行列を1次元配列に変換
X = reshape(X, (1, :))[1, :]
Y = reshape(Y, (1, :))[1, :]
Z = reshape(Z, (1, :))[1, :]

# 3Dグラフを描画
surf(X, Y, Z)

3d-graph.png

このように2つ以上の変数からなる関数に対する微分は、「どの変数に対する微分か」を区別する必要があります

このような微分を偏微分と呼びます

偏微分では、複数ある変数の中でターゲットとする変数を1つに絞り、他の変数は定数として扱います

Input: Julia 1.3.0
"""
@test: 

: f(x_0, x_1) = x_0^2 + x_1^2
: f/x_0 = 2x_0 (x_0)
    x_0 = 3, x_1 = 4  => f/x_0 = 2x_0 = 2 * 3 = 6
"""
# 2変数関数
f(x_0::Float64, x_1::Float64) = x_0^2 + x_1^2

# 1変数(x_0)関数として再定義(x_1 = 4 の定数)
g(x_0::Float64) = f(x_0, 4.0)

# 数値微分: x_0 = 3
## => 6.00000000000378 ≒ 6
println(numeric_diff(g, 3.0))
Output
6.00000000000378
Input: Julia 1.3.0
"""
@test: 

: f(x_0, x_1) = x_0^2 + x_1^2
: f/x_1 = 2x_1 (x_1)
    x_0 = 3, x_1 = 4  => f/x_1 = 2x_1 = 2 * 4 = 8
"""
# 2変数関数
f(x_0::Float64, x_1::Float64) = x_0^2 + x_1^2

# 1変数(x_1)関数として再定義(x_0 = 3 の定数)
g(x_1::Float64) = f(3.0, x_1)

# 数値微分: x_1 = 4
## => 7.999999999999119 ≒ 8
println(numeric_diff(g, 4.0))
Output
7.999999999999119

勾配

全ての変数に対する偏微分を行列形式でまとめたものを勾配と呼びます

上記の2変数関数を例にとると、勾配は以下のように表現されます

$$ \begin{array}{ll} Δf &=& (\frac{∂f}{∂x_0}, \frac{∂f}{∂x_1}) \\ &=& (2x_0, 2x_1) \end{array} $$
Input: Julia 1.3.0
# 数値微分による勾配の実装
## 関数 f(Float64配列), Float64配列 x -> Float64配列 勾配
numeric_gradient(f, x::Array{Float64,1})::Array{Float64,1} = begin
    h = 1e-4 # 10^(-4)
    grad = Array{Float64, 1}(undef, length(x)) # xと同じ長さの配列 [undef, undef, ...] を生成
    # 各変数ごとの数値微分を配列にまとめる
    for i in 1:length(x)
        # 指定indexの変数に対する中心差分を求める
        org = x[i]
        x[i] = org + h
        f1 = f(x) # f([..., x[i] + h, ...])
        x[i] = org - h
        f2 = f(x) # f([..., x[i] - h, ...])
        grad[i] = (f1 - f2) / 2h # i番目の変数に対する数値微分
        x[i] = org # x[i]の値をもとに戻す
    end
    return grad
end

"""
@test:
    2 f(x_0, x_1) = x_0^2 + x_1^2
     (3, 4), (0, 2), (3, 0) 
"""
# 多変数関数
f(x::Array{Float64,1}) = sum(x.^2)

# 点 (3, 4) における勾配
## (2x_0, 2x_1) = (6, 8)
println(numeric_gradient(f, [3.0, 4.0]))

# 点 (0, 2) における勾配
## (2x_0, 2x_1) = (0, 4)
println(numeric_gradient(f, [0.0, 2.0]))

# 点 (3, 0) における勾配
## (2x_0, 2x_1) = (6, 0)
println(numeric_gradient(f, [3.0, 0.0]))
Output
[6.00000000000378, 7.999999999999119]
[0.0, 4.000000000004]
[6.000000000012662, 0.0]

ここで、上記勾配のグラフを描いてみると面白いことが分かります

Input: Julia 1.3.0
# 軸準備
x0 = collect(-2.0 : 0.25 : 2.25)
x1 = collect(-2.0 : 0.25 : 2.25)

# X軸: [x0 × x0]: 18 x 18 Array{Float64,2}
X = hcat(fill(x0, length(x0))...)

# Y軸: [x1 × x1] の転置行列: 18 x 18 Array{Float64,2}
Y = hcat(fill(x1, length(x1))...)
Y = Y'

# Z軸: grid状に勾配(Δf([x0, x1]))計算
## => 18 x 18 Array{Float64, 2}
Z = [numeric_gradient(f, [x, y]) for x in x0, y in x1]

# グラフ描画用に行列を1次元配列に変換
X = reshape(X, (1, :))[1, :]
Y = reshape(Y, (1, :))[1, :]
Z = reshape(Z, (1, :))[1, :]

# quiver: x-y平面に座標(x, y)を基点として、成分(u, v)のベクトルを描く
U = [z[1] for z in Z]
V = [z[2] for z in Z]
quiver(X, Y, -U, -V)

quiver-graph.png

上記グラフのように、当該関数の勾配は、中心点(最適点)(0, 0) の方向を向いたベクトルで表され、中心点から離れるほど勾配の大きさ(傾き)は大きくなっています

すなわち、勾配を求めることによってどの方向にパラメータを移動すれば最適点に近づくのかを知ることができます

上記の例で考えると、例えばパラメータが (1, 1) という値のときの勾配を見れば、そのベクトルの向きが左下 (-, -) 方向であることが分かるため、次の試行パラメータとして (2, 2)(右上方向)のような非効率的なパラメータを選ぶことがなくなります

さらに言えば、勾配の量が大きいほど最適点から離れているということも分かるため、より効率的に試行パラメータを選定することができます

実際には、より複雑な関数においては、勾配の指し示す方向が必ずしも最適点ではないことも多いですが、少なくともその点において関数を減少させる方向を指し示しているのは間違いないです

この性質を利用して、損失関数を減少させ、最適パラメータを探索する方法を勾配降下法と呼びます

勾配降下法を数式で表すと以下のようになります

$$ \begin{array}{ll} \vec{x} = \vec{x} - η \frac{∂f}{∂ \vec{x}} \\ \\ \vec{x}: ベクトル x \\ η: 更新量(学習率) \\ \frac{∂f}{∂ \vec{x}}: 勾配 \end{array} $$

ここで、学習率の値は、前もって任意の定数に決めておく必要があります

このようなパラメータを、機械学習ではハイパーパラメータと呼び、一般的に大きすぎても小さすぎても「良い学習」は行うことができないとされています

そのためニューラルネットワークの学習では、学習率の値を変更しながら、正しく学習できているか確認作業を行う必要があります

Input: Julia 1.3.0
# 勾配降下法の実装
## 損失関数 f, ベクトル x, 学習率 lr, 試行回数 ste_num -> 最適化ベクトル
gradient_descent(f, x::Array{Float64,1}, lr::Float64=0.01, step_num::Int=100)::Array{Float64,1} = begin
    step_num > 0 ? gradient_descent(f, x - lr * numeric_gradient(f, x), lr, step_num - 1) : x
end

"""
@test:  f(x_0, x_1) = x_0^2 + x_1^2 
"""
# 初期パラメータ (-3, 4), 学習率 0.1 で勾配降下法を100回試行
gradient_descent(f, [-3.0, 4.0], 0.1, 100)
Output
2-element Array{Float64,1}:
 -6.111107928998789e-10
  8.148143905314271e-10
Input: Julia 1.3.0
# 学習率が大きすぎる場合
gradient_descent(f, [-3.0, 4.0], 10.0, 100)
Output
2-element Array{Float64,1}:
 -2.5898374737328363e13
 -1.2952486168965398e12
Input: Julia 1.3.0
# 学習率が小さすぎる場合
gradient_descent(f, [-3.0, 4.0], 1e-10, 100)
Output
2-element Array{Float64,1}:
 -2.999999939999995 
  3.9999999199999934

上記のように、学習率が大きすぎると大きく発散してしまい、逆に小さすぎる場合はほとんど学習が進行しません

この学習率のようなハイパーパラメータは、重みやバイアスといった自動的に調整されるパラメータと異なり、人の手で調整しなければならないパラメータです

そのため、一般的には最適なハイパーパラメータを根気よく探すという作業が必要になります


損失関数

2乗和誤差

損失関数として用いられる関数はいくつもあるが、最も有名なものの一つに2乗和誤差があります

2乗和誤差は以下の数式で表されます

$$ \begin{array}{ll} E = \frac{1}{2} \sum_{k}{} (y_k - t_k)^2 \\ \\ k: データの次元数 \\ y_k: ニューラルネットワークの出力(予測値) \\ t_k: 教師データ(正解値) \end{array} $$
Input: Julia 1.3.0
# 2乗和誤差
## Float64行列 予測値, Float64行列 正解値 -> Float64 誤差
mean_squared_error(y::Array{Float64,2}, t::Array{Float64,2})::Float64 = 0.5 * sum((y - t).^2)

"""
@test: 2
- : 2

- : 09 0.1, 0.05, 0.6, 0.0, 0.05, 0.1, 0.0, 0.1, 0.0, 0.0
    => 20.62
- : 09 0.1, 0.05, 0.1, 0.0, 0.05, 0.1, 0.0, 0.6, 0.0, 0.0
    => 70.62
"""
# 「2」を正解ラベルとする行列(教師データ)
## 以下のようにN個のラベルを 1×N Array{0|1, 2} で表現したものを one-hot表現と呼ぶ
t = [0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]

# 予測値: ラベル「2」の確率が最も高い場合
y_2 = [0.1 0.05 0.6 0.0 0.05 0.1 0.0 0.1 0.0 0.0]
## 2乗和誤差 => 小さい(正解に近い)
println(mean_squared_error(y_2, t))

# 予測値: ラベル「7」の確率が最も高い場合
y_7 = [0.1 0.05 0.1 0.0 0.05 0.1 0.0 0.6 0.0 0.0]
## 2乗和誤差 => 大きい(正解から遠い)
println(mean_squared_error(y_7, t))
Output
0.09750000000000003
0.5974999999999999

交差エントロピー誤差

分類問題の損失関数としては、交差エントロピー誤差がよく用いられます

これは、交差エントロピー誤差の計算量が2乗和誤差に比べて遥かに少なく、それでいながら十分に正確な誤差を表現を表すことができるからです

交差エントロピーは以下の数式で表されます

$$ \begin{array}{ll} E = - \sum_{k} t_k \times \ln{y_k} \\ \\ k: データの次元数 \\ y_k: ニューラルネットワークの出力(予測値) \\ t_k: 教師データ(正解値) \ln: 底がe(自然定数)の対数 \end{array} $$

上式は、教師データが one-hot表現(正解ラベルのインデックスのみが 1 の行列表現)である場合、実質正解ラベルの計算しか行いません

そのため、正解ラベルの確率に対する対数(逆指数)表現となります(確率が高いほど加速度的に誤差が小さくなる)

Input: Julia 1.3.0
# 交差エントロピー誤差
## ln(0) = -Inf が発生するのを防ぐため、予測値に微小な値(10^-7)を加算して計算する
cross_entropy_error(y::Array{Float64,2}, t::Array{Float64,2})::Float64 = -sum(t .* log.(y .+ 1e-7))

"""
@test: 
- : 2

- : 09 0.1, 0.05, 0.6, 0.0, 0.05, 0.1, 0.0, 0.1, 0.0, 0.0
    => 20.62
- : 09 0.1, 0.05, 0.1, 0.0, 0.05, 0.1, 0.0, 0.6, 0.0, 0.0
    => 70.62
"""
# 「2」を正解ラベルとする行列(教師データ)
## 以下のようにN個のラベルを 1×N Array{0|1, 2} で表現したものを one-hot表現と呼ぶ
t = [0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]

# 予測値: ラベル「2」の確率が最も高い場合
y_2 = [0.1 0.05 0.6 0.0 0.05 0.1 0.0 0.1 0.0 0.0]
## 交差エントロピー誤差 => 小さい(正解に近い)
println(cross_entropy_error(y_2, t))

# 予測値: ラベル「7」の確率が最も高い場合
y_7 = [0.1 0.05 0.1 0.0 0.05 0.1 0.0 0.6 0.0 0.0]
## 交差エントロピー誤差 => 大きい(正解から遠い)
println(cross_entropy_error(y_7, t))
Output
0.510825457099338
2.302584092994546

ミニバッチ学習

機械学習においては、処理しなければならないデータ数が極めて多いため、効率的に学習を行うためには、複数のデータをひとまとめにしたバッチ単位で学習を行う必要があることはすでに述べた通りです

しかし、膨大な量の訓練データを一つのバッチとして処理しようとすると、今度はメモリ不足になる可能性があります

そこで、訓練データの一部を小さな塊(ミニバッチ)として無作為に選び出し、それらに対して学習処理を行うことをミニバッチ学習と呼びます

ミニバッチ処理の実装

"""
: 
"""
# ミニバッチサイズ
batch_size = 100

# 添字用配列
## 1〜60,000 のランダムな数値からなる1次元配列を batch_size の長さで生成
indexes = rand(1 : size(train_x)[3], batch_size)

# ミニバッチ処理: 対象データ群から batch_size 個のデータを抜き出し
train_x[:, :, indexes]

学習アルゴリズムの実装

これまでに実装した「損失関数」「ミニバッチ」「勾配」「勾配降下法」をまとめることで、ニューラルネットワークの学習アルゴリズムを実装することができます

確率的勾配降下法

ニューラルネットワークの学習手順は以下のようなものが基本となります

  1. ミニバッチ
    • 訓練データからランダムに一部のデータを選び出す(ミニバッチ)
    • 一回の学習においては、このミニバッチの損失関数の値を減少させることを目的とする
  2. 勾配
    • ミニバッチの損失関数を減らすために、各重みパラメータの勾配を算出する
    • 勾配は、損失関数の値を最も減らす方向を示す
  3. パラメータの更新
    • 重みパラメータを勾配方向に微小量だけ更新する
    • 1.に戻って同様の手順を繰り返す

ここで、使用する訓練データをミニバッチとして無作為に選び出していることから、このような学習方法を確率的勾配降下法(Stochastic Gradient Descent)と呼びます

ディープラーニングの多くのフレームワークでは、確率的勾配降下法の頭文字をとってSGDという名前の関数で実装されているのが一般的です

2層ニューラルネットワークの実装

今回は、手書き数字の学習を行うためのニューラルネットワークとして、2層ニューラルネットワーク(隠れ層1つのニューラルネットワーク)を実装することにしましょう

ネットワーク設計は以下の通りです

  • 入力層:
    • 手書き数字の画像データ(サイズ: 28x28)
    • ニューロン数: 28 * 28 = 784
    • 各ニューロンの入力値は 0.0〜1.0 の実数型である必要がある
  • 中間層(隠れ層)1:
    • ニューロン数: 100
    • 活性化関数: シグモイド関数
  • 出力層:
    • ニューロン数: 10(0〜9の数字クラスに分類するため)
    • 活性化関数: ソフトマックス関数
    • 損失関数: 交差エントロピー関数
Input: Julia 1.3.0
"""
2
"""
# Network構造体を継承して2層ニューラルネットワーク実装
TwoLayerNetwork(weight_init_std::Float64=0.01) = Network(2,
    [
        zeros(Float64, 1, 100), # 1x100-Array{Float64,2} バイアス_1: 0行列
        zeros(Float64, 1, 10),  # 1x10-Array{Float64,2} バイアス_2: 0行列
    ],
    [
        rand(UInt8, 784, 100) * weight_init_std, # 784x100-Array{Float64,2} 重み_1: 任意整数 * weight_init_std の乱数行列
        rand(UInt8, 100, 10) * weight_init_std, # 100x10-Array{Float64,2} 重み_2: 任意整数 * weight_init_std の乱数行列
    ]
)

# 推論処理
## Network構造体, 入力信号 -> 出力信号 y
predict(network::Network, x::Array{Float64,2})::Array{Float64,2} = predict(network, sigmoid, softmax, x)

# 推論処理+損失関数
## Network構造体, 入力信号, 教師データ -> 交差エントロピー誤差
loss(network::Network, x::Array{Float64,2}, t::Array{Float64,2})::Float64 = cross_entropy_error(predict(network, x), t)

# 各パラメータの勾配計算
## Network構造体, 入力信号, 教師データ -> 各パラメータの勾配行列をまとめた辞書
numeric_gradient(network::Network, x::Array{Float64,2}, t::Array{Float64,2})::Dict{AbstractString, Array{Float64,2}} = begin
    loss_func = w -> loss(network, x, t)
    Dict(
        "B1" => numeric_gradient(loss_func, network.b[1]),
        "B2" => numeric_gradient(loss_func, network.b[2]),
        "W1" => numeric_gradient(loss_func, network.w[1]),
        "W2" => numeric_gradient(loss_func, network.w[2]),
    )
end

"""
@test: 2
"""
# ダミー入力データ: 0.0〜1.0 の784サイズデータ 100枚分
x = rand(Float64, 100, 784)

# 推論実行
net = TwoLayerNetwork()
y = predict(net, x)

# 各行ごとに列の合計値が1になっているか確認(softmax関数の特性の確認)
[sum(y[row, :]) for row in 1:size(y, 1)]
Output
100-element Array{Float64,1}:
 1.0
 1.0
 ⋮
 1.0

素直に数値微分で勾配降下法を実装すると、上記のようなコードになります

しかし、このコードには重大な問題があります

すなわち、実行時間がかかりすぎるという問題です

以下のコードで、数値微分による勾配降下法にかかる時間を計測してみましょう

Input: Julia 1.3.0
# 訓練用画像データと教師データをロード
## train_x: 特徴量=<画像データ|28x28 グレースケール画像 60,000枚>{28x28x60000 Array{UInt8, 3}}
## train_y: 目的変数=<数値クラス|[0..9]の数値 60,000個>{60000 Array{Int, 1}}
train_x, train_y = MNIST.traindata()

# 学習データをニューラルネットワーク用に前処理
train_x = Array{Float64,3}(reshape(train_x, 1, 28*28, :)) # 1 x 784 x 60000 Array{Float64,3}
train_x = permutedims(train_x, [3, 2, 1]) # 60000 x 784 x 1 Array{Float64,3}
train_x = Array{Float64,2}(reshape(train_x, :, 784)) # 60000 x 784 Array{Float64,2}

# 教師データをone-hot-vector形式に変換
train_y = hcat([[i-1 == y ? 1.0 : 0.0 for i in 1:10] for y in train_y]...)' # 60000 x 10 Array{Float64,2}

# 損失関数の履歴
train_loss_list = []

# ハイパーパラメータ
iters_num = 1 # パラメータ更新回数
train_size = size(train_x, 1) # 学習データ枚数: 60,000
batch_size = 100 # ミニバッチサイズ
learning_rate = 0.1 # 学習率

# 2層ニューラルネットワーク
net = TwoLayerNetwork()

# 学習関数
## Juliaの慣習で、副作用のある(実行ごとに結果が変わる)関数には ! をつける
train!(net::Network) = begin
    # ミニバッチ取得: 対象データ群から batch_size 個のデータを抜き出し
    batch_mask = rand(1:train_size, batch_size)
    batch_x = train_x[batch_mask, :]
    batch_t = train_y[batch_mask, :]
    
    # 勾配の計算
    grad = numeric_gradient(net, batch_x, batch_t)
    
    # パラメータの更新
    net.b[1] -= learning_rate * grad["B1"]
    net.b[2] -= learning_rate * grad["B2"]
    net.w[1] -= learning_rate * grad["W1"]
    net.w[2] -= learning_rate * grad["W2"]
    
    # 学習経過の記録
    loss_value = loss(net, batch_x, batch_t)
    push!(train_loss_list, loss_value)
end

# 学習実行
## @timeマクロでプログラムの実行時間計測が可能
@time train!(net)
Output
 82.919140 seconds (149.05 M allocations: 107.032 GiB, 13.88% gc time)

上記のように、たった1回の学習にも80秒程度かかってしまう(Intel© Core i7-7700 3.6 GHz の場合)ため、10,000回学習させようとすると 222時間(≒9日)程度かかる計算になります

この程度の学習に9日間もかかってしまうようでは、とても実用に耐えることは出来ません

そのため、計算の最適化が必要になります

ここからが本稿で最も重要なテーマである、誤差逆伝播法による計算最適化です


誤差逆伝播法

数値微分の難点

ここまでで、ニューラルネットワークの基本的な学習方法である、確率的勾配降下法は実装できました

しかし、勾配の計算を数値微分によって実装すると、計算に時間がかかるという難点があります

そこで、ここでは、重みパラメータの勾配計算を効率良く行う誤差逆伝播法を実装します

計算グラフ

誤差逆伝播法を視覚的に理解する方法として計算グラフ(Computational graph)というものがあります

計算グラフとは計算の過程をグラフによって表したものです

ここで言うグラフとは、データ構造としてのグラフであり、複数のノードエッジ(ノード間を結ぶ直線)によって表現されるものです

計算グラフで計算問題を解く

以下のような簡単な計算問題を、計算グラフを使って解いてみましょう

  • 問: スーパーで1個100円のりんごを2個買ったとき、支払う金額はいくらか?(消費税は10%とする)

計算グラフは以下のように、ノードごとの計算結果が左から右へ伝わるように表現します

computational_graph01.png

上記では、「×2」や「×1.1」を一つの演算として○(ノード)でくくっていますが、「×」のような演算子は単一のノードで表現するのが望ましいです

この場合、「2」と「1.1」は、それぞれ「りんごの個数」と「消費税」という変数として、ノードの外側に表記します

computational_graph02.png

次に、もう少し複雑な計算を計算グラフで解いてみます

  • 問: スーパーでりんご(100円/個)を2個、みかん(150円/個)を3個買ったとき、支払う金額はいくらか?(消費税は10%とする)

この問題を解くための計算グラフは以下のようになります

computational_graph03.png

このように、計算グラフを使って問題を解くには、

  1. 計算式に使われる要素を演算子と変数に分ける
  2. 演算子と変数をノードとしてエッジでつなぐ
  3. 計算グラフ上で計算を左から右へ進める

という流れで作業します

ここで、「計算を左から右へ進める」処理を順伝播(forward propagation)と呼び、ニューラルネットワークの推論処理に対応します

逆に「計算を右から左へ戻る」処理を逆伝播(back propagation)と呼び、ニューラルネットワークの学習処理に対応します

この逆伝播は、この先、微分を計算するにあたって重要な働きをします

局所的計算と逆伝播による微分

計算グラフを使う利点は大きく以下の2点があります

  • 「局所的な計算」を伝播することにより複雑な計算を行うことができる
    • 各ノードは、全体の計算には関与せず「自分に関係する小さな範囲」の計算だけを行う(局所的計算)
    • これにより計算を単純化することができる
    • また、計算途中の結果をすべて、各ノードで保持することも可能
  • 計算グラフを逆伝播することで微分値を効率よく計算できる
    • 順伝播が「局所的な計算」であるのと同様に、逆伝播は「局所的な微分」を表す
    • これにより微分計算を単純化し、計算速度を向上させることができる

例えば上記の問題について、りんごの値段が値上がりした場合、最終的な支払金額にどのように影響するか知りたいとします

これは「りんごの値段に対する支払金額の微分」を求めることに相当します(りんごの値段を $x$, 支払金額を $L$ とした場合、$∂L/∂x$ を計算することに相当します)

計算グラフの逆伝播によって、この問題を解くと以下のようになります

computational_graph04.png

上記のように、計算グラフを右から左へ計算することで「局所的な微分」を伝達することができます

この結果から「りんごの値段に関する支払金額の微分」の値は 2.2 であると導けます

すなわち、りんごが1円値上がりしたら、最終的な支払金額は2.2円増えることを意味します(正確には、りんごの値段がある微小な値だけ増えたら、最終的な支払金額はその微小な値の2.2倍だけ増加することを意味します)

連鎖律

計算グラフの逆伝播

計算グラフの逆伝播を一般化すると以下のように表すことができます

computational_graph05.png

ここで、上記計算グラフは $y = f(x)$ という計算を表現しています

逆伝播の計算手順は、信号 $E$ に対して、ノードの局所的な微分 $\frac{∂y}{∂x}$ を乗算し、次のノードで伝達する、というものになっています

これは、前述したりんごの支払金額計算で考えると分かりやすいです

  • 最初のノード: f(x) = 2x (x: 前のノードの出力値=りんごの値段, 2: りんごの個数)
  • 最後のノード: g(x) = 1.1x (x: 前のノードの出力値, 1.1: 消費税倍率)
  • 最後のノードの逆伝播: 1 * g’(x) = 1 * 1.1 = 1.1
  • 最初のノードの逆伝播: 1.1 * f’(x) = 1.1 * 2 = 2.2

この計算により効率よく微分値を求めることができるのですが、その理由は連鎖律の原理から説明できます

連鎖律

以下のような計算グラフを考えてみます

computational_graph06.png

この計算を微分すると、以下のような逆伝播で表現されます

computational_graph07.png

ここで、合成関数の定理より $\frac{∂z}{∂z} \frac{∂z}{∂y}$ の $∂z$ は “打ち消し合い”、$\frac{∂z}{∂y}$ となります

同様にして $\frac{∂z}{∂z} \frac{∂z}{∂y} \frac{∂y}{∂x} = \frac{∂z}{∂x}$ となります

このような合成関数の微分の性質を連鎖律と呼びます

連鎖律により、計算の一部を “打ち消す” ことができるため、効率よく微分計算ができるという仕組みです


単純な算術ノードの実装

計算グラフの単純なノードを実装していきましょう

ただし、型名は「ノード」ではなく、ニューラルネットワークの「層(レイヤ)」を意味するものとして ***Layer という名前で実装することにします

# すべてのレイヤの基底となる抽象型: 抽象レイヤ
abstract type AbstractLayer end

加算ノード(加算レイヤ)

まず、加算ノード $z = x + y$ について考えます

このノードの逆伝播(偏微分)は以下のようになります

$$ \begin{array}{ll} \frac{∂z}{∂x} = 1 \\ \frac{∂z}{∂y} = 1 \end{array} $$

従って、この計算グラフは以下のようになります

computational_graph_add.png

上図のように、加算ノードの逆伝播は、上流の値がそのまま分岐して流れていきます

これを実装すると以下のようになります

"""
@note:
    forward, backward 
    => forward!, backward! 
"""
# 加算レイヤ
mutable struct AddLayer <: AbstractLayer end

# 加算レイヤ: 順伝播
## x, y: 上流から流れてきる2値 => 下流に流す値
forward!(layer::AddLayer, x::Float64, y::Float64)::Float64 = x + y

# 加算レイヤ: 逆伝播
## dout: 上流から流れてくる微分値 => 下流に流す2値
backward!(layer::AddLayer, dout::Float64)::Tuple{Float64,Float64} = (dout * 1, dout * 1)

乗算ノード(乗算レイヤ)

同様に、乗算ノード $z = x \times y$ について考えると、このノードの逆伝播(偏微分)は以下のようになります

$$ \begin{array}{ll} \frac{∂z}{∂x} = y \\ \frac{∂z}{∂y} = x \end{array} $$

従って、この計算グラフは以下のようになります

computational_graph_mul.png

すなわち、乗算ノードの逆伝播では、上流から流れてきた微分値に対して、順伝播の “ひっくり返した値” を乗算して流す形になります

これを実装すると以下のようになります

# 乗算レイヤ
mutable struct MulLayer <: AbstractLayer
    x::Float64
    y::Float64
end

MulLayer() = MulLayer(0, 0)

# 乗算レイヤ: 順伝播
## x, y: 上流から流れてきる2値 => 下流に流す値
forward!(layer::MulLayer, x::Float64, y::Float64)::Float64 = begin
    # 順伝播時に値を保持しておく
    layer.x = x
    layer.y = y
    x * y
end

# 乗算レイヤ: 逆伝播
## dout: 上流から流れてくる微分値 => 下流に流す2値
backward!(layer::MulLayer, dout::Float64)::Tuple{Float64,Float64} = (dout * layer.y, dout * layer.x)

活性化関数レイヤの実装

今度は、計算グラフの考え方をニューラルネットワークに適用していきましょう

ReLUレイヤ

ReLU活性化関数は以下の式で表されます

$$ y = \begin{cases} x & (x > 0) \\ 0 & (x \leq 0) \end{cases} $$

よって、微分は以下の式で表されます

$$ \frac{∂y}{∂x} = \begin{cases} 1 & (x > 0) \\ 0 & (x \leq 0) \end{cases} $$

従って、ReLU活性化関数を計算グラフで表すと以下のようになります

computational_graph_relu.png

Input: Julia 1.3.0
# ReLUレイヤ
mutable struct ReluLayer <: AbstractLayer
    mask::AbstractArray{Bool}
    ReluLayer() = new()
end

# ReLUレイヤ: 順伝播
forward!(layer::ReluLayer, x::AbstractArray{Float64})::AbstractArray{Float64} = begin
    mask = layer.mask = (x .<= 0) # x <= 0 の要素をマスキング
    out = copy(x) # 入力値をそのまま出力値にコピー
    out[mask] .= zero(Float64) # x <= 0 でマスキングした要素それぞれに 0 代入
    out
end

# ReLUレイヤ: 逆伝播
backward!(layer::ReluLayer, dout::AbstractArray{Float64})::AbstractArray{Float64} = begin
    # 基本的に上流から流れてきた微分値をそのまま下流へ
    dout[layer.mask] .= zero(Float64) # x <= 0 でマスキングした要素それぞれに 0 代入
    dout
end

"""
@test: ReLU
"""
# ReLUレイヤ: 順伝播
## [1.0 -0.5; -2.0 3.0] => [1.0 0.0; 0.0 3.0]
layer = ReluLayer()
forward!(layer, [1.0 -0.5; -2.0 3.0])
Output
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  3.0
Input: Julia 1.3.0
# ReLUレイヤ: マスク状態(x <= 0 のindex)
layer.mask
Output
2×2 BitArray{2}:
 0  1
 1  0
Input: Julia 1.3.0
# ReLUレイヤ: 逆伝播
## [1.0 1.0; 1.0 1.0] => [1.0 0.0; 0.0 1.0]
backward!(layer, [1.0 1.0; 1.0 1.0])
Output
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0

Sigmoidレイヤ

シグモイド関数は以下の式で表されます

$$ y = \frac{1}{1 + \exp(-x)} $$

これを計算グラフで表すと以下のようになります

computational_graph_sigmoid.png

この計算グラフの逆伝播をステップごとに分解して考えると以下のようになっています

  1. 「/」ノードは $y=\frac{1}{x}$ を表し、この微分は解析的に以下のように解くことができます $$ \begin{array}{ll} \frac{∂y}{∂x} &= -\frac{1}{x^2} \\\
    &= -y^2 \end{array} $$
  2. 「+」ノードは上流から流れてきた微分値をそのまま下流に流します(加算ノードの項参照)
  3. 「exp」ノードは $y=e^x$ を表し、この微分は解析的に以下のように解くことができます $$ \frac{∂y}{∂x} = e^x $$
  4. 「×」ノードは、順伝播時の2値を “ひっくり返して” 乗算するため、ここでは -1 を乗算しています

ここで、$\frac{∂L}{∂y}y^2 e^{-x}$ が順伝播の値 $x$, $y$ のみから計算可能であることに注目すると、上記計算グラフは、以下のように「sigmoid」ノードとしてグループ化することができます

computational_graph_sigmoid2.png

このようにグループ化することで途中の計算を省略することができるため、計算効率を向上させることができるという仕組みです

さらに $\frac{∂y}{∂x}y^2 e^{-x}$ は以下のように展開できるため、さらに整理できます

$$ \begin{array}{ll} \frac{∂y}{∂x}y^2 e^{-x} &= \frac{∂y}{∂x}\frac{1}{(1 + e^{-x})^2}e^{-x} & (∵ y = \frac{1}{1 + e^{-x}}) \\ &= \frac{∂y}{∂x}\frac{1}{1 + e^{-x}}\frac{e^{-x}}{1 + e^{-x}} \\ &= \frac{∂y}{∂x}y(1-y) \end{array} $$

従って、最終的に計算グラフは以下のようにまとめることができます

computational_graph_sigmoid3.png

# Sigmoiodレイヤ
mutable struct SigmoidLayer <: AbstractLayer
    out::AbstractArray{Float64} # Sigmoidレイヤが保持しておかなければならないのは順伝播時の出力値 y
end

# Sigmoiodレイヤ: 順伝播
forward!(layer::SigmoidLayer, x::AbstractArray{Float64})::AbstractArray{Float64} =
    layer.out = 1 ./ (1 .+ exp.(-x))

# Sigmoiodレイヤ: 逆伝播
backward!(layer::SigmoidLayer, dout::AbstractArray{Float64})::AbstractArray{Float64} =
    dount .* layer.out .* (1 .- layer.out)

計算グラフを使って途中計算を整理することによって、上記のように逆伝播(微分)を非常に簡単な計算に置き換えられるということは重要です

Affineレイヤ

ニューラルネットワークの順伝播時で行われる計算は以下のようなものでした

$$ 出力行列 = 入力行列 \cdot 重み行列 + バイアス行列 $$

ここで、行列の内積計算は、幾何学分野でアフィン変換と呼ばれます

そのため、ニューラルネットワーク順伝播時の行列計算を行う層をAffineレイヤとして実装することにします

まず、計算グラフを考えましょう

ここでは、バッチサイズを $N$, 入力値の数を 2, 重みを 2×3行列(入力数: 2, 出力数: 3), バイアスを 1×3行列(出力数: 3)として考えると、計算グラフは以下のようになります

computational_graph_affine.png

なお「dot」ノードは、実質的には「×」ノードと同一と考えて良いため、順伝播時の2入力値を “ひっくり返して” 乗算することで逆伝播の流れを作ることができます

ただし、行列の値をひっくり返すということは、転置行列を作ることを意味します

そのため、入力行列 $X$ の側における逆伝播は $\frac{∂L}{∂Y} \cdot W^T$ となり($W^T$ は $W$ の転置行列)、重み行列 $W$ の側における逆伝播は $X^T \cdot \frac{∂L}{∂Y}$ となっています($X^T$ は $X$ の転置行列)

Input: Julia 1.3.0
# バッチ対応Affineレイヤ
mutable struct AffineLayer <: AbstractLayer
    # AbstractArray型にすることでコンストラクタですべてのメンバの値を指定を不要にする
    W::AbstractArray{Float64,2}
    B::AbstractArray{Float64,2}
    X::AbstractArray{Float64,2}
    dW::AbstractArray{Float64,2}
    dB::AbstractArray{Float64,2}
    
    # コンストラクタで重みとバイアスだけ指定できるように定義
    function AffineLayer(W::AbstractArray{Float64,2}, B::AbstractArray{Float64,2})
        layer = new()
        layer.W = W
        layer.B = B
        layer
    end
end

# Affineレイヤ: 順伝播
forward!(layer::AffineLayer, x::Array{Float64,2})::Array{Float64,2} = begin
    layer.X = x
    x * layer.W .+ layer.B
end

# Affineレイヤ: 逆伝播
backward!(layer::AffineLayer, dout::Array{Float64,2})::Array{Float64,2} = begin
    dX = dout * layer.W'
    layer.dW = layer.X' * dout
    layer.dB = sum(dout, dims=1) # バイアスの逆伝播値: ∂L/∂Y の最初の軸に対する和
    dX
end

"""
@test: Affine
"""
# Affineレイヤ
## 重み行列: 1×3(入力数: 1, 出力数: 3)
## バイアス行列: 1×3(出力数: 3)
layer = AffineLayer([5.0 5.0 5.0], [1.0 2.0 3.0])

# 順伝播
## 入力行列: 2×1(バッチサイズ: 2, 特徴量の数: 1)
x = reshape([0.0; 2.0], 2, 1)
forward!(layer, x) # => 2×3 [1 2 3; 11 12 13]
Output
2×3 Array{Float64,2}:
  1.0   2.0   3.0
 11.0  12.0  13.0

Softmaxレイヤ

Softmax活性化関数は以下の数式で表されます

$$ y_i = \frac{e^{x_i}}{\sum_{i=1}^N e^{x_i}} $$

計算グラフで表すと以下のようになります

computational_graph_softmax.png

この計算グラフの①〜⑥の経路について逆伝播を考えてみましょう

  1. $\frac{∂L}{∂y_i}$
  2. 「×」ノードのため、順伝播時の相方を乗算します $$ \frac{∂L}{∂y_i} \cdot e^{x_i} $$
  3. 「/」ノードのため、順伝播時の出力値を2乗し符号反転したものを乗算し、それらを合計します $$ \begin{array}{ll} -\sum_i^N \frac{∂L}{∂y_i} \cdot e^{x_i} \cdot \frac{1}{S^2} \\\
    = -\sum_i^N \frac{∂L}{∂y_i} \cdot \frac{y_i}{S} &(∵ y_i = \frac{e^{x_i}}{S}) \end{array} $$
  4. 「×」ノードのため、順伝播時の相方を乗算します $$ \frac{∂L}{∂y_i} \cdot \frac{1}{S} $$
  5. 「+」ノードのため、上流値をそのまま流します $$ -\sum_{i=1}^{N}\frac{∂L}{∂y_i} \cdot \frac{y_i}{S} $$
  6. 「exp」ノードのため、順伝播時の出力値を乗算します
    • 今回の場合、上流値が④と⑤の2つあるため、これらの合計値に順伝播時の出力値を乗算することになります $$ \begin{array}{ll} (\frac{∂L}{∂y_i} \cdot \frac{1}{S} - \sum_j^N \frac{∂L}{∂y_j} \cdot \frac{y_j}{S}) \cdot e^{x_i} \\\
      = \frac{1}{S}(\frac{∂L}{∂y_i} - \sum_j^N \frac{∂L}{∂y_j}y_j) \cdot e^{x_i} \\\
      = y_i(\frac{∂L}{∂y_i} - \sum_j^N \frac{∂L}{∂y_j}y_j) &(∵ y_i = \frac{e^{x_i}}{S}) \end{array} $$

Softmaxレイヤ+交差エントロピー誤差レイヤ

上記のようにSoftmaxレイヤの逆伝播計算は少々複雑になってしまっています

しかし実は、交差エントロピー誤差レイヤを組み合わせることで非常に簡単な式に展開することができます

ここで改めて、交差エントロピー誤差の式を掲載します

$$ \begin{array}{ll} E = - \sum_{k} t_k \times \ln{y_k} \\ \\ k: データの次元数 \\ y_k: ニューラルネットワークの出力(予測値) \\ t_k: 教師データ(正解値) \ln: 底がe(自然定数)の対数 \end{array} $$

この計算グラフは以下のようになります

computational_graph_entropy.png

このとき、交差エントロピー誤差の入力値 $y_i$ が、Softmaxレイヤの出力値 $y_i = \frac{e^{x_i}}{S}$ であるとします(Softmaxレイヤと交差エントロピー誤差レイヤが連続的に接続されているとします)

この場合、計算としては、Softmaxレイヤの逆伝播の出力値 $y_i(\frac{∂L}{∂y_i} - \sum_j^N \frac{∂L}{∂y_j}y_j)$ に対して、以下のような置換処理を施せば良いです

$$ \frac{∂L}{∂y_i} → -(\frac{t_i}{y_i})\frac{∂E}{∂y_i} $$

したがって、Softmaxレイヤの逆伝播の出力値は以下のように展開できます

$$ \begin{array}{ll} y_i(\frac{∂L}{∂y_i} - \sum_{j=1}^{N}\frac{∂L}{∂y_j}y_j) \\ = y_i(-(\frac{t_i}{y_i})\frac{∂E}{∂y_i} - \sum_{j=1}^{N}-(y_j\frac{t_j}{y_j})\frac{∂E}{∂y_j}) \\ = -t_i\frac{∂E}{∂y_i} + y_i\sum_{j=1}^{N}t_j\frac{∂E}{∂y_j} \\ = -t_i\frac{∂E}{∂y_i} + y_i\frac{∂E}{∂y_i} &(∵ t: one-hot表現のため正解データのみ 1 で、それ以外のデータが 0 となる) \\ = (y_i - t_i)\cdot \frac{∂E}{∂y_i} \end{array} $$

以上より、Softmax+交差エントロピー誤差レイヤは以下のように表現できます

computational_graph_softmax-entropy.png

# Softmax with Cross entropy error Layer
mutable struct SoftmaxWithLossLayer <: AbstractLayer
    loss::Float64 # 最終結果(誤差)
    y::AbstractArray{Float64} # SoftmaxレイヤからCrossEntropyErrorレイヤへのデータ
    t::AbstractArray{Float64} # 教師データ
    SoftmaxWithLossLayer() = new()
end

# 順伝播関数
## y: 予測結果行列, t: 教師データ行列 => loss: 誤差
forward!(layer::SoftmaxWithLossLayer, y::Array{Float64,2}, t::Array{Float64,2})::Float64 = begin
    layer.t = t
    y = layer.y = softmax(y)
    layer.loss = cross_entropy_error(y, t)
end

# 逆伝播関数
## 教師データは 1×N 行列である想定であるため、バッチサイズで除算する必要あり
backward!(layer::SoftmaxWithLossLayer, dout::Float64=1.0)::Array{Float64,2} = (layer.y .- layer.t) .* dout ./ size(layer.y, 1)

誤差逆伝播法によるニューラルネットワーク

これまでに実装したレイヤを組み合わせて、誤差逆伝播法に対応したニューラルネットワークを実装してみましょう

ここでは、以下のような二層ニューラルネットワークを構築します

  • 入力層:
    • 手書き数字の画像データ(サイズ: 28x28)
    • ニューロン数: 28 * 28 = 784
    • 各ニューロンの入力値は 0.0〜1.0 の実数型である必要がある
  • 中間層(隠れ層)1:
    • Affineレイヤ => ReLU活性化レイヤ
    • ニューロン数: 100
  • 出力層:
    • Affineレイヤ => SoftmaxWithLoss誤差伝播レイヤ
    • ニューロン数: 10(0〜9の数字クラスに分類するため)

二層ニューラルネットワークの実装

Input: Julia 1.3.0
# 二層ニューラルネットワーク
mutable struct TwoLayerNet
    affine_layer1::AffineLayer
    relu_layer::ReluLayer
    affine_layer2::AffineLayer
    softmax_layer::SoftmaxWithLossLayer
end

TwoLayerNet(input_size::Int, hidden_size::Int, output_size::Int, weight_init_std::Float64=0.01) = begin
    W1 = weight_init_std .* rand(Float64, input_size, hidden_size)
    b1 = zeros(Float64, 1, hidden_size)
    W2 = weight_init_std .* rand(Float64, hidden_size, output_size)
    b2 = zeros(Float64, 1, output_size)
    a1lyr = AffineLayer(W1, b1)
    relu1lyr = ReluLayer()
    a2lyr = AffineLayer(W2, b2)
    softmaxlyr = SoftmaxWithLossLayer()
    # TwoLayerNet(W1, b1, W2, b2, a1lyr, relu1lyr, a2lyr, softmaxlyr)
    TwoLayerNet(a1lyr, relu1lyr, a2lyr, softmaxlyr)
end

# 二層ニューラルネットワーク: 順伝播(分類まで)
predict!(net::TwoLayerNet, x::Array{Float64,2})::Array{Float64,2} = begin
    a1 = forward!(net.affine_layer1, x)
    z1 = forward!(net.relu_layer, a1)
    a2 = forward!(net.affine_layer2, z1)
    # 本来なら最後にsoftmax関数を噛ませるが
    ## あえて確率を求めなくても a2 の大小で分類できるため計算を省略する
    # softmax(a2)
    a2
end

# 誤差逆伝播: 更新用構造体
struct TwoLayerNetGrads
    W1::Array{Float64,2}
    b1::Array{Float64,2}
    W2::Array{Float64,2}
    b2::Array{Float64,2}
end

# 逆伝播による微分計算
gradient!(net::TwoLayerNet, x::Array{Float64,2}, t::Array{Float64,2}) = begin
    # forward
    loss!(net, x, t)
    # backward
    dout = one(Float64)
    dz2 = backward!(net.softmax_layer, dout)
    da2 = backward!(net.affine_layer2, dz2)
    dz1 = backward!(net.relu_layer, da2)
    da1 = backward!(net.affine_layer1, dz1)
    TwoLayerNetGrads(net.affine_layer1.dW, net.affine_layer1.dB, net.affine_layer2.dW, net.affine_layer2.dB)
end

# パラメータ更新
apply_gradient!(net::TwoLayerNet, grads::TwoLayerNetGrads, learning_rate::Float64) = begin
    net.affine_layer1.W -= learning_rate .* grads.W1
    net.affine_layer1.B -= learning_rate .* grads.b1
    net.affine_layer2.W -= learning_rate .* grads.W2
    net.affine_layer2.B -= learning_rate .* grads.b2
end

"""
@test: 
"""
# 訓練用画像データと教師データをロード
## train_x: 特徴量=<画像データ|28x28 グレースケール画像 60,000枚>{28x28x60000 Array{UInt8, 3}}
## train_y: 目的変数=<数値クラス|[0..9]の数値 60,000個>{60000 Array{Int, 1}}
train_x, train_y = MNIST.traindata()

# 学習データをニューラルネットワーク用に前処理
train_x = Array{Float64,3}(reshape(train_x, 1, 28*28, :)) # 1 x 784 x 60000 Array{Float64,3}
train_x = permutedims(train_x, [3, 2, 1]) # 60000 x 784 x 1 Array{Float64,3}
train_x = Array{Float64,2}(reshape(train_x, :, 784)) # 60000 x 784 Array{Float64,2}

# 教師データをone-hot-vector形式に変換
train_y = Array{Float64,2}(hcat([[i-1 == y ? 1.0 : 0.0 for i in 1:10] for y in train_y]...)') # 60000 x 10 Array{Float64,2}

# 損失関数の履歴
train_loss_list = []

# ハイパーパラメータ
batch_size = 3 # ミニバッチサイズ
learning_rate = 0.1 # 学習率

# 二層ニューラルネットワーク
## 入力サイズ: 28×28 = 784, 隠れ層: 100ニューロン, 出力層: 10ニューロン
net = TwoLayerNet(784, 100, 10)

# 逆伝播による微分
## @timeマクロでプログラム実行時間計測
grads = @time gradient!(net, train_x[1:batch_size, :], train_y[1:batch_size, :])
Output
0.110585 seconds (408.63 k allocations: 20.579 MiB)

前回、数値微分で行った場合は同一の処理に80秒近くかかっていましたが、今回の逆伝播による微分は 1秒未満の時間で計算できています

これが、逆伝播法による計算最適化の効果であり、本稿で最も重要なテーマです

計算最適化が問題なく完了したので、続いて、実際にMNISTの学習をニューラルネットワークに対して行わせてみましょう

Input: Julia 1.3.0
"""
MNIST
"""
# 履歴
train_loss_list = [] # 損失関数の履歴
train_acc_list = [] # 訓練用データでの精度の履歴
test_acc_list = [] # テスト用データでの精度の履歴

# ハイパーパラメータ
iters_num = 10000 # パラメータ更新回数
train_size = size(train_x, 1) # 学習データ枚数: 60,000
batch_size = 100 # ミニバッチサイズ
learning_rate = 0.1 # 学習率

# 2層ニューラルネットワーク
## 入力サイズ: 28×28 = 784, 隠れ層: 50ニューロン, 出力層: 10ニューロン
net = TwoLayerNet(784, 50, 10)

# 学習関数
train!(net::TwoLayerNet) = begin
    # ミニバッチ取得: 対象データ群から batch_size 個のデータを抜き出し
    batch_mask = rand(1:train_size, batch_size)
    batch_x = train_x[batch_mask, :]
    batch_y = train_y[batch_mask, :]
    
    # 誤差逆伝播法によって勾配を求める
    grads = gradient!(net, batch_x, batch_y)
    
    # 更新
    apply_gradient!(net, grads, learning_rate)
    
    # 学習経過の記録
    loss_value = loss!(net, batch_x, batch_y)
    push!(train_loss_list, loss_value)
end

# 予測精度算出関数
accuracy(net::TwoLayerNet, x::Array{Float64,2}, t::Array{Float64,2})::Float64 = begin
    pred = predict!(net, x)
    ind = getindex.(argmax(pred; dims=2), [2]) # 予測値: 行ごとの最大値のindex(1..10)を取得
    t_ind = getindex.(argmax(t; dims=2), [2]) # 教師データ: 行ごとの最大値のindex(1..10)を取得
    length(ind[ind .== t_ind]) / size(t, 1)
end

# 指定回数学習を繰り返す関数
train!(
    net::TwoLayerNet, iters_num::Int, iters_per_epoch::Int,
    train_x::Array{Float64,2}, train_y::Array{Float64,2},
    test_x::Array{Float64,2}, test_y::Array{Float64,2}
) = begin
    for i = 1:iters_num
        train!(net)
        if i % iter_per_epoch == 1
            train_acc = accuracy(net, train_x, train_y)
            test_acc = accuracy(net, test_x, test_y)
            push!(train_acc_list, train_acc)
            push!(test_acc_list, test_acc)
            println("$(i): train_acc=$(train_acc) / test_acc=$(test_acc)")
        end
    end
end

# テスト用画像データと教師データをロード
test_x, test_y = MNIST.testdata()

# データをニューラルネットワーク用に前処理
test_x = Array{Float64,3}(reshape(test_x, 1, 28*28, :))
test_x = permutedims(test_x, [3, 2, 1])
test_x = Array{Float64,2}(reshape(test_x, :, 28*28))
test_y = Array{Float64,2}(hcat([[i-1 == y ? 1.0 : 0.0 for i in 1:10] for y in test_y]...)')

# 1エポックあたりの繰り返し回数
iter_per_epoch = Int(max(train_size / batch_size, 1))

# 初期化
train_loss_list = [] # 損失関数の履歴
train_acc_list = [] # 訓練用データでの精度の履歴
test_acc_list = [] # テスト用データでの精度の履歴
net = TwoLayerNet(784, 50, 10)

# 学習実行: 10,000回実行
@time train!(net, iters_num, iter_per_epoch, train_x, train_y, test_x, test_y)
Output
1: train_acc=0.09871666666666666 / test_acc=0.098
601: train_acc=0.9006333333333333 / test_acc=0.9027
1201: train_acc=0.9204166666666667 / test_acc=0.9221
1801: train_acc=0.9331166666666667 / test_acc=0.9325
2401: train_acc=0.9407666666666666 / test_acc=0.9405
3001: train_acc=0.9474833333333333 / test_acc=0.9458
3601: train_acc=0.9536333333333333 / test_acc=0.9507
4201: train_acc=0.9595 / test_acc=0.9543
4801: train_acc=0.963 / test_acc=0.9595
5401: train_acc=0.9653666666666667 / test_acc=0.9606
6001: train_acc=0.9665666666666667 / test_acc=0.9615
6601: train_acc=0.9697 / test_acc=0.9626
7201: train_acc=0.97225 / test_acc=0.9659
7801: train_acc=0.9737166666666667 / test_acc=0.9675
8401: train_acc=0.97465 / test_acc=0.966
9001: train_acc=0.9769333333333333 / test_acc=0.9661
9601: train_acc=0.9786333333333334 / test_acc=0.9681

 32.471889 seconds (14.40 M allocations: 28.023 GiB, 1.07% gc time)

上記のように 10,000回という大量の学習回数にも関わらず、30秒程度で実行が完了しています

このくらい速ければ、実用にも十分耐えられるでしょう

ニューラルネットワークの学習推移

今回の学習では、ニューラルネットワークの損失関数と予測精度を履歴として保存しました

その推移をグラフ化すると以下のようになります

Input: Julia 1.3.0
"""
@plot: 
"""
# PyPlot利用
using PyPlot

# 損失関数の推移
plot(train_loss_list)

learning-graph.png

今回損失関数として交差エントロピー誤差を採用していますが、これは分類問題における予測結果と教師データ(正解データ)の誤差の大きさを表します

グラフを見ると、誤差が徐々に減少しており、学習が順調に進んでいることが分かります

続いて予測精度のグラフを確認してみましょう

Input: Julia 1.3.0
# 予測精度の推移
plot(train_acc_list; label="train_acc")
plot(test_acc_list; label="test_acc")
legend()

acc-graph.png

今回、訓練用データとテスト用データ両方を使った予測精度を追っていますが、これは過学習(オーバーフィッティング)を防ぐのが目的です

ディープラーニングにおいては、学習回数を増やせば増やすほど、ニューラルネットワークのパラメータは最適化されていきますが、それが訓練用データに対して最適化されすぎると未知のデータに対しての予測精度が低下してしまうことがあります

このような現象を過学習と呼びます

これを防ぐための手段の一つとして、学習回数を適切な回数に調整するという方法があります

すなわち、訓練用データに対する予測精度が上がり続けているのに、テスト用データ(未知のデータ)に対する予測精度が下がり始めた辺りの学習回数を限度とする方法です

今回の場合で言うと、おそらく学習回数 4,000回目くらいから過学習が起こり始めているので、このくらいの回数で学習を切り上げるのが吉です

というわけで、このニューラルネットワークモデルの最終的な予測精度は 95.5%程度 と言えます

Summary

本稿では、DeepLearningがそのような計算によって学習を行うのか、その仕組みを一から実装してみました

だいたい1時間程度で読める記事になっていますが、ここで改めて学習内容を振り返っておきましょう

  1. パーセプトロン
    • パーセプトロンはディープラーニング(ディープニューラルネットワーク)の起源とされるアルゴリズム
    • 0 or 1 の二値からなる入力信号を受け取り、特定の計算処理を行うことで 0 or 1 の二値を出力する
    • 単層パーセプトロンでは線形分離可能な問題しか解けないが、層を重ねることで線形分離不可能な問題も解ける
  2. ニューラルネットワーク
    • パーセプトロンの構造を拡張した計算アルゴリズム
    • 計算処理を活性化関数という関数に置き換えることでより汎用的な計算が可能となっている
    • パーセプトロン同様、層を重ねることで「深い」計算が可能となり、それがディープラーニング・アルゴリズムの骨子である
  3. 勾配降下法
    • ディープラーニングにおける学習の本質はパラメータの最適化であり、その手法として広く用いられているのが勾配降下法である
    • 勾配とは偏微分計算で求められる微小変化量のことであり、関数グラフの傾きのことである
    • 勾配はその点における最適点の方向を向くため、勾配の示す方向にパラメータを調整していくことで最適化ができるという仕組み
  4. 誤差逆伝播法
    • 勾配計算において、数値微分を利用すると計算処理に莫大な時間がかかってしまうため、その対策のために用いられる手法
    • 計算グラフの逆伝播が微分計算になるという性質を利用して微分を行う
    • 連鎖律の原理により重複する計算等を省略できるため、高速に微分計算を行うことができる

その他、ミニバッチ学習や確率的勾配降下法など、学習効率を向上させるための各種テクニックが存在しますが、概念的に特に重要なのは上記の4つです

とは言え、各社が出しているディープラーニング・ライブラリは優秀であるため、直接本稿の知識を必要とする場面は少ないかもしれません

しかし、より深くディープラーニングを学ぼうする際には、このように原始的・本質的な知識が必要になってきますので、覚えておいて損はないと筆者は考えています

Avatar
Ameno Yoya
Webプログラマ

経験はログに残して初めて資産となる

comments powered by Disqus