深層学習の基礎（修正版）

最終更新：2022/2/24

記事の内容

以前書いた記事の改訂版です. 内容を修正し, Fluxでの実装例も追加しました. 主に[1]を参考にしました.

概要

問題設定と方針

機械学習, 特に教師あり学習は, データを説明する適当な関数を定める問題として, 以下のように定式化されます.

【目標】

手元にサイズ\(N\)のデータ\(X=\left\{ \boldsymbol{x}_n\right\}_{n=1}^{N}, Y=\left\{ \boldsymbol{y}_n\right\}_{n=1}^{N}\)があるとき, 適当なパラメータ\(\boldsymbol{w}\)を持つ写像\(\Phi_{w}\colon \mathbb{R}^{D_x}\to\mathbb{R}^{D_y}\)を, そこそこの予測性能を保ちつつ, \(\boldsymbol{y}_n\simeq\Phi_{w}(\boldsymbol{x}_n)\)が達成されるように定めたい.

"適当な関数"として, 予め"線形関数"や"周期関数"のように限定して議論しても良いのですが, 当然, 写像の表現能力は制限されます. そこで, 複雑な表現を可能にするために, 以下のような形を考えてみましょう：

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

これはいわゆるニューラルネットワークです. この場合, 写像のパラメータ\(\boldsymbol{w}\)は以下のように表されます：

\begin{equation} \boldsymbol{w} = \left\{ W^{(2)},\dots,W^{(L)},\boldsymbol{b}^{(2)},\dots,\boldsymbol{b}^{(L)}\right\}. \end{equation}

ここで, \(W\)達は重みと呼ばれ, \(\boldsymbol{b}\)達はバイアスと呼ばれます. また, \(\phi\)達は活性化関数と呼ばれ, 通常, tanhやシグモイド関数などの非線形関数が選択されます. 要するに, ニューラルネットワークとは, 線形変換と非線形変換を交互に繰り返した写像です：

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

ニューラルネットワークは深層学習の基本的なツールであり, 深層学習の方針は以下の一文に集約されます.

【深層学習の方針】

データを用いて, ニューラルネットワーク\(\Phi_{w}(\boldsymbol{x})=\boldsymbol{a}^{(L)}\)がいい感じになるように\(\boldsymbol{w}\)を調整する.

"いい感じ"の基準には様々考えられますが, データへの当てはまりの良さを要求するのであれば, 以下のような関数で測るのが自然でしょう：

\begin{equation} E(\boldsymbol{w}) = \sum_{n=1}^NE_n(\boldsymbol{w}),\quad E_n(\boldsymbol{w}) = \frac{1}{2}\|\boldsymbol{y}_n-\Phi_{w}(\boldsymbol{x}_n)\|^2. \end{equation}

データによく当てはまっているほど, コスト関数\(E\)の値は小さくなりますから, 以下の最適化問題を解くのが目標です：

\begin{equation} \underset{\boldsymbol{w}}{\mathrm{argmin}}\ E(\boldsymbol{w}). \end{equation}

この最適化問題はニューラルネットワークの訓練と呼ばれ, 解析的に解くのは難しいので, 最急降下法で反復的にパラメータを更新していきます：

\begin{equation} \boldsymbol{w}^{(k+1)} = \boldsymbol{w}^{(k)} - \alpha_k \nabla_{\theta}E\left( \boldsymbol{w}^{(k)}\right),\quad (k=0,1,2,\dots). \end{equation} 問題となるのは勾配の計算です. 目的関数への寄与が複雑な大量のパラメータの勾配を求めるのは一般に計算コストが高く, 工夫が必要です. そこで登場するのが誤差逆伝播法です. まあ, 実際にプログラムを書く際には自動微分で良いのですが... 次節では, 以上の方針を実現するための訓練アルゴリズムの全容をご紹介します.

アルゴリズム

本節では, ニューラルネットワークの訓練アルゴリズムを構成する順伝播, 誤差逆伝播法, 確率的勾配降下法について説明します.

順伝播

コスト関数の勾配計算には, コスト関数の現在の値が必要であり, そのコスト関数はデータとニューラルネットワークの二乗誤差の和で定義したので, 予測値\(\boldsymbol{a}^{(L)}\)を求める必要が生じます. データから予測値を計算する過程を順伝播と言いますが, 要するにニューラルネットワークの定義式にデータを放り込んで計算するだけです.

誤差逆伝播

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

\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( \boldsymbol{\theta} \right) = \boldsymbol{\delta}^{(l)}_n{\boldsymbol{a}^{(l-1)}_n}^{\mathrm{T}}, \quad (l=2,\cdots,L)\\ &\nabla_{b^{(l)}} E_n\left( \boldsymbol{\theta} \right) = \boldsymbol{\delta}^{(l)}_n, \quad (l=2,\cdots,L). \end{align}

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

確率的勾配降下法

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

\begin{align} W^{(l)}_{k+1} &= W^{(l)}_k - \alpha_k \nabla_{W^{(l)}} E_n\left(\boldsymbol{w}^{(k)} \right) \\ \boldsymbol{b}^{(l)}_{k+1} &= \boldsymbol{b}^{(l)}_k - \alpha_k \nabla_{b^{(l)}} E_n\left(\boldsymbol{w}^{(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,0]^{\mathrm{T}}\), 青バツをベクトル\([0,1]^{\mathrm{T}}\)で表しておきます.

【コード6の実行結果】

隠れ層が2層のニューラルネットワークで分類を試みます. 隠れ層のユニット数は10としておきます. ニューラルネットワークの出力を2次元にしておいて, 第1成分の方が大きければ赤丸, 第2成分の方が大きければ青バツと予測することにします. 活性化関数は全てシグモイド関数, コスト関数はこれまで通り二乗和誤差とします. ニューラルネットワークを式で表示しておくと以下のようになります：

\begin{equation} \sigma\left( W^{(4)}\sigma\left( W^{(3)}\sigma(W^{(2)}+\boldsymbol{b}^{(2)})+\boldsymbol{b}^{(3)}\right) + \boldsymbol{b}^{(4)}\right). \end{equation}

実験結果

以上の設定のもと, 実際に訓練をしてみましょう. コードは本記事の末尾の方に載せています. 訓練回数は\(10^5\)回とし, ステップ幅を一律に0.1としました. そこそこうまく分類できていますね. 分類の境界はやや粗く, コスト関数も思ったほど減少していません. 素朴な実装としてはまずまずでしょうか...

【コード7の実行結果】

Flux

Fluxとは

以上のようにゼロからコードを書いても良いのですが, やはり問題ごとにいちいち書き換えるのも面倒だし, 速度的にも限界があります. Juliaの場合, Flux.jlが便利です. 以下に簡単な実装例を示します. もちろん, アレンジすれば単純な深層学習モデル以外の発展的なモデル簡単に実装できます. 中の人は, 変分推論法や変分オートエンコーダなどもFluxで書いています.

基本的な使い方

Flux.jlを使ってみましょう. まず, ニューラルネットワークを定義してみます. Chainの中に, 各層を表すDenseを入れます. Denseにはインプットの次元, アウトプットの次元, 活性化関数を指定します.

【ニューラルネットワークの作成・初期化】
#define and initialize neural network
nn = Chain(
        Dense(2,10,sigmoid),
        Dense(10,10,sigmoid),
        Dense(10,10,sigmoid),
        Dense(10,2,sigmoid)
        )

以下のような結果が返ってきます. これで初期化完了です.

【ニューラルネットワークの作成・初期化：結果】
Chain(
  Dense(2, 10, sigmoid),                # 30 parameters
  Dense(10, 10, sigmoid),               # 110 parameters
  Dense(10, 10, sigmoid),               # 110 parameters
  Dense(10, 2, sigmoid),                # 22 parameters
)                   # Total: 8 arrays, 272 parameters, 1.562 KiB.

次に, コスト関数を定義します.

【コスト関数の定義】
#define cost function
function cost(x,y,nn)
    norm(y-nn(x))^2
end

とりあえず現時点でのコスト関数の値を計算してみます.

【訓練前のコスト計算】
#cost before training
sum([cost(X[:,n],Y[:,n],nn) for n in 1:N])
【訓練前のコスト計算：結果】
8.114218200862476

この状態から訓練して, コスト関数を減少させてみましょう. 訓練前に更新方法を指定します. 今回はADAMと呼ばれる方法を使ってみます. ステップ幅を適応的に決めてくれるやつです.

【optimizerの呼び出し】
#define optimizer
opt = ADAM()
【optimzerの呼び出し：結果】
ADAM(0.001, (0.9, 0.999), 1.0e-8, IdDict{Any, Any}())

では, 確率的勾配降下方をイメージして, 10番目のデータだけを用いて, パラメータを1回だけ更新してみます. train関数には, コスト関数, モデルのパラメータ, 使うデータ, 更新方法を指定します. ニューラルネットワークのパラメータはparamsで取り出すことができ, パラメータだけを別で管理できます.

【訓練】
#training neural network
n = 10
Flux.train!(cost,params(nn),[(X[:,n],Y[:,n])],opt)

パラメータを更新できたので, この時点でのコスト関数の値を計算してみます. うまくいっていれば値は先ほどより減少しているはずです.

【訓練後のコスト計算】
#cost after training
sum([cost(X[:,n],Y[:,n],nn) for n in 1:N])
【訓練後のコスト計算：結果】
8.123143853089669

データを1個分しか使っていないので, 今回は失敗しました. コスト関数の値が増加しています. 以上をまとめて次のように書けば, ニューラルネットワークの訓練が実行できます.

【ニューラルネットワークの訓練】
#set the random seed
Random.seed!(42)

#initialize the model
nn = Chain(
        Dense(2,10,sigmoid),
        Dense(10,10,sigmoid),
        Dense(10,10,sigmoid),
        Dense(10,2,sigmoid)
        )

#cost function 
function cost(x,y,nn)
    norm(y-nn(x))^2
end

#training
function my_train(data,nn,n_train)
    @unpack X,Y,N = data
    costs = zeros(n_train)
    opt = ADAM(0.01)
    @showprogress for k in 1:n_train
        n = rand(1:N)
        Flux.train!(cost,params(nn),[(X[:,n],Y[:,n])],opt)
        costs[k] = mean([cost(X[:,n],Y[:,n]) for n in 1:N])
    end
    return costs
end

分類してみる （Fluxで実装）

先ほどの分類問題を, 今書いたコードで解いてみましょう.

実験結果

先程の素朴な実装とは多少条件が異なります. 層は1層増え, optimizerにはADAMを指定しています. 実行結果は以下の通りです. コスト関数は先ほどに比べ, 大きく減少しています. また, 分類の境界も滑らかですね.

【コード8の実行結果】

コード

【Juliaコード1; 初期化】
#mathematics
using LinearAlgebra
using Flux

#statistics
using Random
using Statistics

#visualize
using Plots
pyplot()

#macros
using UnPack
using ProgressMeter
【Juliaコード2; ニューラルネットワークの定義】
#create abstract layer, Neural Network
abstract type AbstractLayer end
abstract type AbstractNeuralNetwork end

#create layer
mutable struct Layer <: AbstractLayer
    d_in::Int32
    d_out::Int32
    zvec::Vector{Float32}
    avec::Vector{Float32}
    error::Vector{Float32}
    weights::Matrix{Float32}
    bias::Vector{Float32}
end

#create neural network
mutable struct NeuralNetwork <: AbstractNeuralNetwork
    n_hiddens::Vector{Int32}
    layers::Vector{Layer}
end

#initialize the layers
function init_layer(d_in, d_out, W_init, b_init)
    return Layer(d_in, d_out, zeros(d_out), zeros(d_out), zeros(d_out), W_init, b_init)
end
init_layer(d_in, d_out) = init_layer(d_in, d_out, randn(d_out,d_in), zeros(d_out))

#initialize the neural network
function init_nn(n_hiddens)
    L = length(n_hiddens)
    layers = Vector{Layer}(undef,L-1)
    for l in 1:L-1
        d_in = n_hiddens[l]
        d_out = n_hiddens[l+1]
        layers[l] = init_layer(d_in, d_out)
    end
    return NeuralNetwork(n_hiddens,layers)
end
【Juliaコード3; 順伝播と逆伝播】
#activation function 
function sigmoid(x)
   1 ./ (1 .+ exp.(-x)) 
end

#derivative of sigmoid
function sigmoid_prime(x)
    sigmoid(x) .* (1 .- sigmoid(x))
end

#forward pass
function forward!(x,nn,L)
    avec_prev = x
    for l in 1:L-1
        nn.layers[l].zvec = nn.layers[l].weights * avec_prev + nn.layers[l].bias
        nn.layers[l].avec = sigmoid(nn.layers[l].zvec)
        avec_prev = nn.layers[l].avec
    end
end

#back propagation
function back_prop!(y,nn,L)
    nn.layers[L-1].error = sigmoid_prime(nn.layers[L-1].zvec) .* (nn.layers[L-1].avec-y)
    for l in reverse(1:L-2)
        nn.layers[l].error = sigmoid_prime(nn.layers[l].zvec) .* (nn.layers[l+1].weights' * nn.layers[l+1].error)
    end
end
【Juliaコード4; 訓練用の関数】
#loss function 
function cost(X,Y,nn,L)
    N = size(X,2)
    costs = zeros(N)
    for n in 1:N
        forward!(X[:,n],nn,L)
        costs[n] = norm(Y[:,n]-nn.layers[L-1].avec)^2
    end
    mean(costs)
end

#stochastic gradient descent
function SGD_update!(nn,ε,xvec,L)
    avec_prev = xvec
    for l in 1:L-1
        nn.layers[l].weights -= ε * nn.layers[l].error * avec_prev'
        nn.layers[l].bias -= ε * nn.layers[l].error
        avec_prev = nn.layers[l].avec
    end
end

#training
function my_train(data, nn, n_train, ε)
    @unpack X,Y,N = data
    L = length(nn.n_hiddens)
    costs = zeros(n_train)
    @showprogress for k in 1:n_train
        n = rand(1:N)
        x = X[:,n]
        y = Y[:,n]
        forward!(x,nn,L)
        back_prop!(y,nn,L)
        SGD_update!(nn,ε,x,L)
        costs[k] = cost(X,Y,nn,L)
    end
    return costs
end
【Juliaコード5; 実験用の関数】
#plot the data and return the figure
function plot_data(X, Y)
    _,N = size(X)
    fig = 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[k]==1
            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
    return fig
end

function plot_data(fig, X, Y)
    _,N = size(X)
    fig = plot!(xticks=0:0.2:1, xlim=[0,1], yticks=0:0.2:1, ylim=[0,1], aspect_ratio=:equal, legend=false)
    for k in 1:N
        if Y[k]==1
            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
    return fig
end

#prediction
function pred(x₁,x₂,nn)
    L = length(nn.n_hiddens)
    forward!([x₁,x₂],nn,L)
    return nn.layers[L-1].avec[1]>nn.layers[L-1].avec[2] ? 0 : 1
end
【Juliaコード6; データの作成】
#create the data
N = 16
X = [
    0.11 0.30 0.35 0.32 0.50 0.71 0.78 0.92  0.10 0.15 0.40 0.50 0.60 0.70 0.80 0.90
    0.10 0.93 0.54 0.33 0.19 0.78 0.51 0.32  0.80 0.40 0.75 0.95 0.50 0.30 0.10 0.80
]
labels = vcat(zeros(div(N,2)), ones(div(N,2)))
Y = zeros(2,N)
for n in 1:N
    if labels[n] == 0.0
        Y[1,n] = 1
    else
        Y[2,n] = 1
    end
end
data = (X=X,Y=Y,N=N)

#plot the data
fig1 = plot_data(X,labels)
savefig(fig1, "figs-DL/fig1.png")
【Juliaコード7; 素朴な実装で実験】
#set the random seed
Random.seed!(42)

#initialize the neural network
n_hiddens = [2,10,10,10,2]
nn = init_nn(n_hiddens)

#training neural entwork
ε = 0.1
n_train = Int(1e5)
@time costs = my_train(data, nn, n_train, ε)

#prediction
p1 = plot(0:0.1:1, 0:0.1:1, (x₁,x₂)->pred(x₁,x₂,nn), st=:heatmap, c=cgrad(:coolwarm), alpha=0.6, clim=(0,1))
p1 = plot_data(p1,X,labels)
p2 = plot(costs,xlabel="iter",ylabel="costs",title="costs",xscale=:log10,yscale=:log10,label=false)
fig2 = plot(p1,p2,size=(1000,400))
savefig(fig2, "figs-DL/fig2.png")
【Juliaコード8; Flux.jlで実験】
#set the random seed
Random.seed!(42)

#initialize the model
nn = Chain(
        Dense(2,10,sigmoid),
        Dense(10,10,sigmoid),
        Dense(10,10,sigmoid),
        Dense(10,2,sigmoid)
        )

#cost function 
function cost(x,y,nn)
    norm(y-nn(x))^2
end

#training
function my_train(data,nn,n_train)
    @unpack X,Y,N = data
    costs = zeros(n_train)
    opt = ADAM(0.01)
    @showprogress for k in 1:n_train
        n = rand(1:N)
        Flux.train!(cost,params(nn),[(X[:,n],Y[:,n])],opt)
        costs[k] = mean([cost(X[:,n],Y[:,n]) for n in 1:N])
    end
    return costs
end

#train the model
n_train = Int(1e5)
@time costs = my_train(data,nn,n_train)

#prediction 
function pred(x₁,x₂,nn)
    tmp = nn([x₁,x₂])
    return tmp[1]>tmp[2] ? 0 : 1
end

#visualize
p1 = plot(0:0.02:1, 0:0.02:1, (x₁,x₂)->pred(x₁,x₂,nn), st=:heatmap, c=cgrad(:coolwarm), alpha=0.6, clim=(0,1))
p1 = plot_data(p1,X,labels)
p2 = plot(costs,xlabel="iter",ylabel="costs",title="costs",xscale=:log10,yscale=:log10,label=false)
fig3 = plot(p1,p2,size=(1000,400))
savefig(fig3, "figs-DL/fig3.png")
参考文献
      [1]C.F.Higham, D.J.Higham, Deep Learning: An Introduction for Applied Mathematicians, arXiv, 2018
      [2]岡谷貫之, 深層学習, 講談社, 2019