深層学習の基礎

記事の内容


本記事では, 深層学習の初歩的な解説と素朴に実装したプログラムをご紹介します.

概要

基本的なアイデア

機械学習の動機の1つは, データを作り出す関数を特定することです. 例えば線形回帰なら, 線形の関数からデータが生成されると仮定して, その正体を特定します. しかし, 予め「線形」という制約が入ります. 線形性により, 関数の表現力は著しく低下し, 現実のデータでは多少無理が生じます. そこで, 関数の形状を指定することなくデータから学習できれば理想的です. その実現方法の1つが深層学習です. 深層学習では単層のNeural Networkを多数つなげたモデルで学習を行います. 以下, 教師あり学習のみを扱います.

一般的な定式化をしておきます. 手元に, 教師ありデータ\(\{(\boldsymbol{x}_n, \boldsymbol{y}_n)\}_{n=1}^N\)があるとします. このデータから, 入出力の関係 \begin{equation} \boldsymbol{y}_n \simeq \Phi (\boldsymbol{x}_n), \quad \Phi : \mathbb{R}^{D_x}\to\mathbb{R}^{D_y} \end{equation} を学習します. この関数について何か手がかりが欲しいところです. そこで, 入力に対する基本的な変換, 線形変換と非線形変換を複数回施すことにします. すなわち, 以下のような関数を考えます.

\begin{equation} \Phi(\boldsymbol{x}) = \phi^{(L)}\left(W^{(L)}\cdots+ \phi^{(3)}\left( W^{(3)}\phi^{(2)}\left( W^{(2)}\boldsymbol{x}+\boldsymbol{b}^{(2)}\right)+\boldsymbol{b}^{(3)}\right)+\cdots+ \boldsymbol{b}^{(L)} \right) \end{equation}

\(W^{(l)}\)は行列で, \(\boldsymbol{b}^{(l)}\)はベクトルです. それぞれ, 重み, バイアスと呼びます. また, \(\phi^{(l)}\)は非線形関数です. これは深層学習の枠組みでは, 活性化関数と呼ばれます. 表現を簡潔にするため, 上の式を次のように段階的に書きます.

\begin{align} \boldsymbol{a}^{(1)} &= \boldsymbol{x} \\ \boldsymbol{z}^{(2)} &= W^{(2)}\boldsymbol{a}^{(1)} + \boldsymbol{b}^{(2)} \\ \boldsymbol{a}^{(2)} &= \phi\left( \boldsymbol{z}^{(2)}\right) \\ \vdots\\ \boldsymbol{z}^{(L)} &= W^{(L)}\boldsymbol{a}^{(L-1)} + \boldsymbol{b}^{(L)} \\ \boldsymbol{a}^{(L)} &= \phi\left( \boldsymbol{z}^{(L)}\right) \\ \end{align}

これを, \(L\)層の順伝播NeuralNetworkと言います. 深層学習では, この関数がデータの入出力の関係\(\Phi\left( \boldsymbol{x}\right) = \boldsymbol{a}^{(L)}\)を与えるように, パラメータ\(W^{(2)}, \cdots, W^{(L)},\boldsymbol{b}^{(2)},\cdots,\boldsymbol{b}^{L}\)を決めます. 次節では, これらのパラメータを決めるためのアルゴリズムを構成していきます.

アルゴリズム

本節では, 深層学習のアルゴリズムの根幹をなす3つの要素, 順伝播, 逆伝播, 確率的勾配降下法を簡単に解説します.

順伝播

前節の深層学習モデルにより, 重みとバイアスが与えられれば出力\(\Phi\left( \boldsymbol{x}\right)=\boldsymbol{a}^{(L)}\)が計算できます. この計算を順伝播と言います. 深層学習では, この順伝播の計算は何度か行います. 実装は簡単で, 前節のように線形変換と非線形変換を繰り返すだけです. より詳細に定式化しておくと, 以下のようになります.

\begin{align} \boldsymbol{a}^{(1)} &= \boldsymbol{x} \in \mathbb{R}^{D_1} \\ \boldsymbol{z}^{(l)} &= W^{(l)}\boldsymbol{a}^{(l-1)} + \boldsymbol{b}^{(l)},\quad W^{(l)}\in\mathbb{R}^{D_l\times D_{l-1}},\quad \boldsymbol{b}^{(l)}\in\mathbb{R}^{D_l} \quad (l=2,\cdots,L) \\ \boldsymbol{a}^{(l)} &= \phi^{(l)} \left( \boldsymbol{z}^{(l)} \right),\quad \phi^{(l)}:\mathbb{R}^{D_l}\to\mathbb{R}^{D_l} \quad (l=2,\cdots,L) \end{align}

ここで, \(D_x=D_1, D_y=D_L\)とします. 以下, \(\boldsymbol{x}_n\)に対応する\(\boldsymbol{a}, \boldsymbol{z}\)を, \(\boldsymbol{a}_n, \boldsymbol{z}_n\)と表記します.

方針

パラメータ\(W^{(l)}, \boldsymbol{b}^{(l)}\)を定めるための方針を述べます. 基本的には最適化問題として定式化します. \(\Phi\)のデータへの当てはまり具合を評価するために, 次のような指標を用意します.

\begin{equation} E(W, \boldsymbol{b}) = \sum_{n=1}^{N} E_n(W, \boldsymbol{b}) , \quad E_n(W, \boldsymbol{b}) = \frac{1}{2} \| \boldsymbol{y}_n-\Phi\left(\boldsymbol{x}_n\right)\|^2 \end{equation}

ノルムは通常の2ノルムです. この指標をコスト関数と言います. 深層学習では, コスト関数を最小値するようなパラメータを求めます. ということで, 以下の最小化問題を解くことが目標です.

\begin{equation} \underset{W^{(l)}, \boldsymbol{b}^{(l)}}{\min} E(W, \boldsymbol{b}) \end{equation}

最小化問題を解く方法は様々あります. とりあえず, 最もシンプルな最急降下法を用います. 最急降下法は, 目的関数の勾配を利用して反復的に近似解を構成する方法です. 第\(k+1\)反復では, 以下のように更新します. \(\varepsilon\)は学習率と呼ばれます.

\begin{align} W^{(l)}_{k+1} &= W^{(l)}_k - \varepsilon \nabla_{W^{(l)}} E\left( W^{(l)}_k, \boldsymbol{b}^{(l)}_k \right) \\ \boldsymbol{b}^{(l)}_{k+1} &= \boldsymbol{b}^{(l)}_k - \varepsilon \nabla_{b^{(l)}} E\left( W^{(l)}_k, \boldsymbol{b}^{(l)}_k \right) \end{align}

ということで, コスト関数の勾配が求まればよさそうです. しかし, 重みやバイアスのコスト関数への寄与は一般に複雑です. 有限差分近似などの方法もありますが, 精度や計算コストの観点から避けた方が良いでしょう. そこで, 次の誤差逆伝播法を用います.

誤差逆伝播

誤差逆伝播法は, コスト関数の勾配を求める方法です. ここで, 次のようなベクトルを用意します.

\begin{equation} \boldsymbol{\delta}^{(l)}_n = \nabla_{\boldsymbol{z}^{(l)}} E_n(W, \boldsymbol{b}) \end{equation}

連鎖律から, 次のような式が得られます[1].

\begin{align} &\boldsymbol{\delta}^{(L)}_n = {\phi^{(l)}}^{\prime}\left( \boldsymbol{z}^{(L)}_n\right) \odot \left( \boldsymbol{a}^{(L)}_n - \boldsymbol{y}_n\right) \\ &\boldsymbol{\delta}^{(l)}_n = {\phi^{(l)}}^{\prime}\left( \boldsymbol{z}^{(l)}_n\right) \odot \left( {W^{(l+1)}}^{\mathrm{T}}\boldsymbol{\delta}^{(l+1)}_n \right),\quad (l=2,\cdots,L-1)\\ &\nabla_{W^{(l)}} E_n\left( W, \boldsymbol{b} \right) = \boldsymbol{\delta}^{(l)}_n{\boldsymbol{a}^{(l-1)}_n}^{\mathrm{T}}, \quad (l=2,\cdots,L-1)\\ &\nabla_{b^{(l)}} E_n\left( W, \boldsymbol{b} \right) = \boldsymbol{\delta}^{(l)}_n, \quad (l=2,\cdots,L-1) \end{align}

\(\odot\)は要素ごとの積(Hadamard積)です. 勾配を計算するには\(\boldsymbol{\delta}_n\)が必要です. しかし, \(\boldsymbol{\delta}^{(l)}\)は添え字の大きい方から計算する必要があります. これは順伝播とは逆方向の演算なので, 逆伝播と呼ばれます. この計算は本質的に自動微分と同一です. したがって, 実装する際には自動微分のライブラリを呼び出せばOKです(本記事のプログラムでは自動微分を使わずに, 素朴に実装しています).

確率的勾配降下法

これで勾配の計算はできそうです. 深層学習では, 普通の最急降下法ではなく, 確率的勾配降下法を用います. これは, 各反復でデータの中から1サンプル(もしくは一部のサンプル)を用いて最急降下法を適用する方法です. 一部のサンプルを使うと反復ごとに勾配が異なるため, 局所解から抜け出すことができると期待できます. このサンプルは例えば一様分布に従って選びます. 1サンプルを使う場合の確率的勾配降下法の更新式は以下の通りです.

\begin{align} W^{(l)}_{k+1} &= W^{(l)}_k - \varepsilon \nabla_{W^{(l)}} E_n\left( W^{(l)}_k, \boldsymbol{b}^{(l)}_k \right) \\ \boldsymbol{b}^{(l)}_{k+1} &= \boldsymbol{b}^{(l)}_k - \varepsilon \nabla_{b^{(l)}} E_n\left( W^{(l)}_k, \boldsymbol{b}^{(l)}_k \right) \end{align}

勾配は先ほどの誤差逆伝播法から計算できます.

アルゴリズム

深層学習のアルゴリズムの基本的な形は以下の通りです.

【Neural Networkの訓練】

Set the number of units, number of layers \(L\) and activate function.
Initialize \(W^{(l)}, \boldsymbol{b}^{(l)}\) and set the number of iterations and learning rate \(\varepsilon\).
for \(i=1,2, \cdots\)
Choose \(k\) from \(\{1,2,\cdots,N\}\) at uniformly random.
#forward propagete
Set \(\boldsymbol{a}^{(1)}=\boldsymbol{x}_k\).
for \(l=2,3,\cdots,L\)
Compute \(\boldsymbol{z}^{(l)}, \boldsymbol{a}^{(l)}\).
end
#back propagate
Compute \(\boldsymbol{\delta}^{(L)}_k\).
for \(l=L-1,\cdots,3,2\)
Compute \(\boldsymbol{\delta}^{(l)}_k\).
end
#SGD
for \(l=2,3,\cdots,L\)
Update \(W^{(l)}\) and \(\boldsymbol{b}^{(l)}\).
end
end

例1: 分類問題

まず, 分類問題に取り組みます. 下図のように赤い点と青い点の, 2種類が存在するとします. この2つのカテゴリの境界を求めます. 赤い点は\(\mathrm{Beta}(1.1, 3)\), 青い点は\(\mathrm{Beta}(1.9, 1.1)\)から発生させた乱数を座標とします. データサイズは\(N=14\)です.

【コード3の実行結果】

入力層から\(D_1=2, D_2=3, D_3=3, D_4=2\)とします. 活性化関数は, いずれも以下のシグモイド関数です. (厳密には, シグモイド関数を成分とするベクトル値関数です. )

\begin{equation} \sigma(x) = \frac{1}{1+e^{-x}}\ \end{equation}

学習後の関数\(\Phi\)は2次元のベクトル値関数です. 最終的には, 第1成分の方が大きければ赤色の点と同じカテゴリー, 第2成分の方が大きければ青色の点と同じカテゴリーとみなします. 学習率は\(\varepsilon=0.1\), 反復回数は\(100000\)とします. 以下に学習後の分類の境界を示します. 赤い領域は, 赤い点と同じカテゴリに入るとみなされます.

【コード5の実行結果】

分類はうまくいってそうです. 以下に, コスト関数の勾配のノルムの増減の様子を示します.

【コード6の実行結果】

例2: 回帰問題

2つ目の例は回帰問題です. 以下の関数をもとにノイズを加えたデータ点を作成します.

\begin{equation} f(x) = \frac{\sin (2x)}{2x} \end{equation}

ノイズは正規分布から発生させます. また, データサイズは\(N=100\)です.

【コード7の実行結果】

学習は, 入力層から\(D_1=1, D_2=10, D_3=1\)とし, 活性化関数はシグモイド関数です. 学習率は\(\varepsilon=0.1\), 反復回数は\(100000\)とします. 以下に結果を示します.

【コード9の実行結果】

端のあたりの再現に失敗しました... まだ層数やユニット数, 活性化関数, 学習率などは調整中です. 以下にコスト関数の勾配のノルムの増減の様子を示します.

【コード10の実行結果】

コード

【Juliaコード; インポート】
using LinearAlgebra
using ForwardDiff
using Distributions
using Random
using Statistics
using ProgressMeter
using Plots
pyplot()
【Juliaコード2; 関数・複合型の定義】
#各層
mutable struct Layer
    din::Int64 #dimension of input
    dout::Int64 #dimension of output
    zvec::Vector{Float64} #transfomed
    avec::Vector{Float64} #activate
    error::Vector{Float64} #error δ
    weights::Matrix{Float64} #weighs matrix
    bias::Vector{Float64} #bias vector
    act_func::Function #activate function
    act_der::Function #derivative of activate function
end

#Neural Network
mutable struct NeuralNetwork
    nus::Vector{Int64} #number of hidden units in each layer including input and output layer
    layers::Vector{Layer} #number of layers including output layer
end

#initialize Layer
function initLayer(din::Int64, dout::Int64, Winit::Matrix{Float64}, binit::Vector{Float64}, act_func)
    act_der(x) = ForwardDiff.derivative(act_func, x)
    Layer(din, dout, zeros(dout), zeros(dout), zeros(dout), Winit, binit, act_func, act_der)
end

#initialize Neural Network
function initNeuralNetwork(nus::Vector{Int64}, act_funcs)
    #number of layers
    L = length(nus)
    
    #initialize layers except for input layer
    layers = Vector{Layer}(undef, L-1)
    for l in 1:L-1
        din = nus[l]
        dout = nus[l+1]
        layers[l] = initLayer(din, dout, zeros(dout, din), zeros(dout), act_funcs[l])
    end
    NeuralNetwork(nus, layers)
end

#check the dimensions
function check_dims(X::Matrix{Float64}, Y::Matrix{Float64}, nn::NeuralNetwork)
    #size of X and Y
    Dx, N = size(X)
    Dy, M = size(Y)
    input_nus = nn.nus[begin]
    output_nus = nn.nus[end]
    
    #sample size check
    if N != M
        println("Sample size must be same.")
    end
    
    #input and output dimension must be equal to number of units of each layer
    if input_nus != Dx
        println("Input dimension must be same.")
    elseif output_nus != Dy
        println("Output dimension must be same.")
    end
    return Dx, N
end

#initialize the parameters of each layers
function init_params(nn::NeuralNetwork, L::Int64)
    for l in 1:L-1
        layer = nn.layers[l]
        layer.weights = rand(layer.dout, layer.din)
        layer.bias = rand(layer.dout)
    end
end

#activation
function activate(layer::Layer, avec::Vector{Float64})
    layer.zvec = layer.weights * avec + layer.bias
    layer.avec = layer.act_func.(layer.zvec)
end

#forward pass
function forward_pass(x::Vector{Float64}, nn::NeuralNetwork, L::Int64)
    act_prev = x
    
    for l in 1:L-1
        layer = nn.layers[l]
        activate(layer, act_prev)
        act_prev = layer.avec
    end
end

#back propagation
function back_prop_update(layer::Layer, weights::Matrix{Float64}, error::Vector{Float64})
    layer.error = layer.act_der.(layer.zvec) .* (weights' * error)
end

function back_prop(y::Vector{Float64}, nn::NeuralNetwork, L::Int64)
    #第L層の計算
    layer = nn.layers[L-1]
    layer.error = layer.act_func.(layer.zvec) .* (layer.avec-y)
    wgh_next = layer.weights
    err_next = layer.error
    
    #第L層より前の計算
    for l in reverse(1:L-2)
        layer = nn.layers[l]
        back_prop_update(layer, wgh_next, err_next)
        wgh_next = layer.weights
        err_next = layer.error
    end
end

#Stochastic Gradient Descent
function mySGD_update(layer::Layer, ϵ::Float64, avec::Vector{Float64})
    layer.weights -= ϵ * layer.error * avec'
    layer.bias -= ϵ * layer.error
end

function mySGD(x::Vector{Float64}, ϵ::Float64, nn::NeuralNetwork, L::Int64)
    act_prev = x
    for l in 1:L-1
        layer = nn.layers[l]
        mySGD_update(layer, ϵ, act_prev)
        act_prev = layer.avec
    end
end

#calculate the cost
function calc_cost(X::Matrix{Float64}, Y::Matrix{Float64}, nn::NeuralNetwork, N::Int64, L::Int64)
    costvec = zeros(N)
    for n in 1:N
        forward_pass(X[:,n], nn, L)
        costvec[n] = norm(Y[:,n]-nn.layers[end].avec)^2
    end
    return mean(costvec)
end

#train Neural Network
function trainNeuralNetwork(X::Matrix{Float64}, Y::Matrix{Float64}, n_iter::Int64, ϵ::Float64, nn::NeuralNetwork)
    #number of layers except for input layers
    L = length(nn.nus)
    
    #sample size and the dimensions
    Dx,N = check_dims(X, Y, nn)

    #initialize the parameter
    init_params(nn, L)
    
    #initialize an array to save cost
    costs = zeros(n_iter)
    
    @showprogress for i in 1:n_iter
        #choose the sample
        k = rand(1:N)
        x = X[:,k]
        y = Y[:,k]
        
        #forward pass
        forward_pass(x, nn, L)
        
        #back propagation
        back_prop(y, nn, L)
        
        #SGD
        mySGD(x, ϵ, nn, L)
        
        #save the cost
        costs[i] = calc_cost(X, Y, nn, N, L)
    end
    return nn, costs
end

#Neural Network as a function 
function predictNeuralNetwork(x::Vector{Float64}, nn::NeuralNetwork)
    L = length(nn.nus)
    forward_pass(x, nn, L)
    nn.layers[end].avec
end

function predictNeuralNetwork(x::Float64, nn::NeuralNetwork)
    L = length(nn.nus)
    forward_pass([x], nn, L)
    nn.layers[end].avec[1]
end

#Neural Network classifier
inA(y::Vector{Float64}) = y[1]>y[2] ? 1.0 : 0.0
classifyNeuralNetwork(x::Vector{Float64}, nn::NeuralNetwork) = inA(predictNeuralNetwork(x, nn))
【Juliaコード3; データ作成: 分類問題】
#データサイズ
N = 14

#データ生成(xは各座標, yはラベル)
Random.seed!(46)
X = hcat(rand(Beta(1.1,3), 2, div(N,2)), rand(Beta(1.9,1.1), 2, div(N,2)))
Y = vcat(vcat(ones(div(N,2)), zeros(div(N,2)))', vcat(zeros(div(N,2)), ones(div(N,2)))')

fig1 = plot(xticks=0:0.2:1, xlim=[0,1], yticks=0:0.2:1, ylim=[0,1], aspect_ratio=:equal, title="data", legend=false)
for k in 1:N
    if Y[1,k]>Y[2,k]
        plot!([X[1,k]], [X[2,k]], st=:scatter, markershape=:circle, markersize=10, color="red")
    else
        plot!([X[1,k]], [X[2,k]], st=:scatter, markershape=:x, markersize=10, color="blue")
    end
end
#savefig(fig1, "figs-DL/fig1.png")
【Juliaコード4; NNの訓練: 分類問題】
#sigmoid function
σ(x) = 1/(1+exp(-x))

#initialize the Neural Network
nus = [2, 3, 3, 2] #number of units
act_funcs = [σ, σ, σ] #activate functions for each layer
nn = initNeuralNetwork(nus, act_funcs)

#train the Neural Network
n_iter = Int(1e5) #number of iterations
ϵ = 0.1 #learning rate
nn, costs = trainNeuralNetwork(X, Y, n_iter, ϵ, nn)
【Juliaコード5; 学習の結果: 分類問題】
#predict
T = 30
y_pred = zeros(2)
preds = zeros(T,T)
x1s = range(0,1,length=T)
x2s = range(0,1,length=T)
for i in 1:T
    for j in 1:T
        preds[i,j] = classifyNeuralNetwork([x1s[i], x2s[j]], nn)
    end
end

#visualize
fig2 = plot(xticks=0:0.2:1, xlim=[0,1], yticks=0:0.2:1, ylim=[0,1], aspect_ratio=:equal, title="data", legend=false)
for k in 1:N
    if Y[1,k]>Y[2,k]
        plot!([X[1,k]], [X[2,k]], st=:scatter, markershape=:circle, markersize=10, color="red")
    else
        plot!([X[1,k]], [X[2,k]], st=:scatter, markershape=:x, markersize=10, color="blue")
    end
end
for i in 1:T
    for j in 1:T
        x1 = x1s[i]
        x2 = x2s[j]
        if preds[i,j] == 1.0
            plot!([x1], [x2], st=:scatter, color="red", markerstrokewidth=0, alpha=0.4)
        else
            plot!([x1], [x2], st=:scatter, color="blue", markerstrokewidth=0, alpha=0.4)
        end
    end
end

#savefig(fig2, "figs-DL/fig2.png")
【Juliaコード6; コスト関数の様子: 分類問題】
#cost function
fig3 = plot(1:n_iter, costs, xscale=:log10, yscale=:log10, leg=false, xlabel="iter", ylabel="cost", title="cost function")
#savefig(fig3, "figs-DL/fig3.png")
【Juliaコード7; データの作成: 回帰問題】
#推定したい関数
true_func(x) = sin(2*x)/x/2

#データの作成
Random.seed!(42)
N = 100
X = 4*(rand(N) .- 0.5)
Y = true_func.(X) + 0.2*(rand(N) .-0.5)

#プロットで確認
xs = -2:0.01:2
fig4 = plot(xlabel="x", ylabel="y")
plot!(xs, x->true_func(x), label="true func", color=:blue)
plot!(X, Y, label="data", st=:scatter, alpha=0.5, markerstrokewidth=0, color=:blue)
#savefig(fig4, "figs-DL/fig4.png")
【Juliaコード8; NNの訓練: 回帰問題】
#データのサイズを変換
Xt = reshape(X, 1, length(X))
Yt = reshape(Y, 1, length(Y))

#initialize the Neural Network
Random.seed!(42)
nus = [1, 10, 1] #number of units
act_funcs = [σ, σ] #activate functions for each layer
nn = initNeuralNetwork(nus, act_funcs)

#train the Neural Network
n_iter = Int(1e5) #number of iterations
ϵ = 0.1 #learning rate
nn, costs = trainNeuralNetwork(Xt, Yt, n_iter, ϵ, nn)
【Juliaコード9; 学習の結果: 回帰問題】
fig5 = plot(xs, x->predictNeuralNetwork(x,nn), label="estimated", legend=:topright, xlabel="x", ylabel="y", color=:red)
plot!(xs, x->true_func(x), label="true func", color=:blue)
plot!(X, Y, label="data", st=:scatter, alpha=0.5, markerstrokewidth=0, color=:blue)
#savefig(fig5, "figs-DL/fig5.png")
【Juliaコード10; コスト関数の様子】
#cost function
fig6 = plot(1:n_iter, costs, xscale=:log10, yscale=:log10, leg=false, xlabel="iter", ylabel="cost", title="cost function")
#savefig(fig6, "figs-DL/fig6.png")
参考文献

      [1]C.F.Higham, D.J.Higham, Deep Learning: An Introduction for Applied Mathematicians, arXiv, 2018
      [2]須山敦志, ベイズ深層学習, 講談社, 2020
      [3]岡谷貫之, 深層学習, 講談社, 2019