ディープラーニングをJuila実装しながら理解する
まえがき
ディープラーニングってどうやってプログラミングするんだ？ってのを解決したかったというのと，Julia言語に慣れるためにやってみた．
コードはここにまとめた．
https://github.com/Soyukke/julia-ayame-machine-learning
内容
掛け算，足し算，Mean Squared Errorを実装してIrisデータセットを使って学習するというところまで．
計算のおおまかな流れ
- 掛け算や足し算を組み合わせてモデルの構築
- lossを計算する
- backward (誤差逆伝搬)する
- パラメータを更新する
モデルを構築するとき，グラフが構築されていくことをイメージした実装をする．
たとえば足し算なら x1 + x2 = yとなるので，yの下にx1とx2がぶら下がっているイメージ．
型
演算子に入力される型をTensorとする．
演算fに対して，Tensorを入力してTensorが出力される．Tensorは繋がりを持っているので，structで定義する．
# ベースとなる型 abstract type Tensor end # 計算グラフの末端を意味する型 abstract type LeafTensor <: Tensor end # 演算子から出力されるTensorの型 abstract type FunctionalTensor <: Tensor end # 具体 supertypeがTensorなタイプ # VariableTensorはニューラルネットにおいてパラメータとなり，最適化されるTensor mutable struct VariableTensor <: LeafTensor x # Tensorが保持する値, ベクトルやマトリックス grad # 誤差逆伝搬で計算した勾配を蓄積するための変数 end # ConstantTensor mutable struct ConstantTensor <: LeafTensor x # Tensorが保持する値, ベクトルやマトリックス grad # 誤差逆伝搬で計算した勾配を蓄積するための変数 end # 以降は演算して出力されるTensorのstructを定義していく # y = x1 + x2 mutable struct AddTensor <: FunctionalTensor x # y.x tensors::Array{Tensor, 1} # x1, x2 end # y = x1 * x2 mutable struct MulTensor <: FunctionalTensor x # y.x tensors::Array{Tensor, 1} # x1, x2 end
forwardの実装
実装の試案
演算子が行われるとき，演算子f::Functionにいくつかの入力がある．
具体例として2項演算子を考える．
fの引数は二つのTensor型のオブジェクトなのでt1::Tensor, t2::Tensorが引数である．
出力はTensorだが，具体的な演算子に関する型である．例えば足し算であればスーパータイプがTensorであるAddTensor型が出力される．以下のような関数で演算を実装できる．
function add(t1::Tensor, t2::Tensor) x3 = t1.x + t2.x return AddTensor(x3, [t1, t2]) end
掛け算など他の演算も同様の実装方法で可能である．
また，FunctionalTensorがスーパータイプである型，AddTensor, MulTensorなどは微分を計算することができる．backward計算時に微分が計算され，下に繋がるTensorに微分が渡されていく．
しかし，渡す必要がない場合があるのである．それはConstantTensor同士の演算である．
function add(t1::ConstantTensor, t2::ConstantTensor) x3 = t1.x + t2.x return ConstantTensor(x3) end
足し算のみに限らず，他の演算子でもこれはおそらく成り立つ．演算子を定義するたびにConstantTensorを入力したときの関数を定義するのは面倒だ．そこで，関数の抽象化を考えてみる．
演算する関数の抽象化
引数として演算子の種類としてstrutctの名前であるEnzansiTensor，forwar計算するときの関数，および入力されるTensorのリストtensorsを与えることで計算できる．
入力されるtensorsがすべてConstantTensor型である場合のみ`ConstatnTensorがreturnされるようにすればよい．
# 抽象化された関数 # 演算子を計算するときにConstant同士の演算の出力はConstantであるようにする # Array{T, 1} where T <: Tensor という書き方でTensor以下のタイプのTのArrayを網羅できる function abstract_operator(EnzansiTensor::DataType, f::Function, tensors::Array{T, 1} where T <: Tensor) # Tensor.xを展開してfに代入する y = f(map(tensor -> tensor.x, tensors)...) return EnzansiTensor(y, tensors) end # ConstantTensor同士の演算 function abstract_operator(EnzansiTensor::DataType, f::Function, tensors::Array{ConstantTensor, 1}) y = f(map(tensor -> tensor.x, tensors)...) return ConstantTensor(y) end
Array{T, 1} where T <: Tensorという構文でTensorのサブタイプのArrayをすべて網羅することができる．
仮にArray{Tensor, 1}であると，ReLUTensor <: Tensor = trueであるが，Array{ReLUTensor, 1} <: Array{Tensor, 1} = falseであるためabstract_operatorですべてのTensorのArrayをカバーできない．
具体的な関数を実装するときは以下のようにして実装できる．
function ReLU(t::Tensor) return abstract_operator(ReLUTensor, x -> x .* (x .> 0), [t]) end
backward の実装
backward演算では，計算グラフの頂点から下へ下へ微分値を引数として渡していくようなイメージ
t2 = f(t1)
t3 = g(t2)
という演算があったとき，t3を最下層にあるt1のパラメータに関する微分を計算してみる．ただし，$t3$はスカラであることを想定する．実用上も最終的には損失関数を通してスカラになるのでこれでよい．
backwardを計算するときの計算の流れを追ってみる. sizeをtensor.xとtensor.gradを同じにするために毎回転置している．表記が面倒なので$\frac{\part y}{\part x}$をdydxという風に表記する．
dt3dt3 = 1dt3d2 = dt3dt3 * dt3dt2dt3dt1 = dt3dt2 * dt2dt1
よって，一つ上のTensorの微分がわかれば最下層の微分の値が計算できる
演算t2 = f(t1)を実装するときにその微分であるdt2dt1も同時に実装すればよい
例えば足し算ではt3 = t1 + t2で出てくるdt3dt1とdt3dt2の二通りの微分がある．n項演算子はn個の偏微分が存在する．これをgrad_fnという関数で計算してそれぞれの偏微分をリストで返す．このとき，一つ上の偏微分を引数として受け取ることにする．
t4 = f(t3)という関数があったときにdt4dt1, dt4dt2の二通りが存在してそれぞれdt4dt1 = dt4dt3 * dt3dt1, dt4dt2 = dt4dt3 * dt3dt2である．なので，引数は一つ上のTensorの微分であるdt4dt3．
# t4 = f(t3) # t3 = t1 + t2 # 足し算なのでAddTensor型である function grad_fn(t3::AddTensor, dt4dt3) # 足し算されたTensorを取り出す t1, t2 = t3.tensors # それぞれの偏微分 # dt2/dt1 * dt1/dt3 = dt4dt3 * I = dt4dt3 # dt2/dt1 * dt1/dt4 = dt4dt3 * I = dt4dt3 return [dt4dt3, dt4dt3] end
このようにgrad_fnを各演算に対して実装してあげる．
最下層へのbackward
さて，偏微分の実装はできたが値が返されるだけで何もしていない．
必要なのはパラメータに関する偏微分の値をVariableTensorのgradというメンバ変数(メンバ変数という呼び方であってるのかわからない)に足しこんでいくということが重要だ．あとはその勾配を使って最急降下法などの最適化アルゴリズムを使ってパラメータの更新ができる．
実は簡単でこんな感じで実装できる
# 最初はbackward(tensor)って呼び出す． function backward(t::Tensor) backward(t, 1) end # 再帰 function backward(t::Tensor, grad_up) # 微分はtが持つTensorsの数だけあるから，リストで返すことにする． vgrad = grad_fn(t, grad_up) # 偏微分と対応するTensorのペアをつくる # [(tensor1, dt3dt1), (tensor2, dt3dt2)] といった感じ pair = zip(t.tensors, vgrad) # 掘り下げ, 各ペアに対してbackwardを呼び出す map(x -> backward(x[1], x[2]), pair) end # ConstantTensorとかは偏微分を受け取っても何もすることがない function backward(t::LeafTensor, grad_up) end # VariableTensorは勾配を蓄積していく function backward(t::VariableTensor, grad_up) # 偏微分は分母にベクトル/マトリックスが来ているのでt.xとの対応は転置の関係にある t.grad += grad_up' end
VariableTensorの変数を書き換える
パラメータを更新するとき, VariableTensorだけxを更新するので，関数fをVariableTensorにだけ適用する関数を用意する．
function tensor_map(f::Function, t::Tensor) # 任意の関数をVariableTensorに作用させる # 掘り下げていく map(t -> tensor_map(f, t), t.tensors) end function tensor_map(f::Function, t::LeafTensor) # 末端だけどVariableTensor以外のtypeはパラメータを変えない nothing end function tensor_map(f::Function, t::VariableTensor) # VariableTensorにだけ作用させる f(t) end
VariableTensorの勾配を0にする
function zero_grad(tensor::Tensor) f = t -> t.grad = zero(t.grad) tensor_map(f, tensor) end
最適化関数
シンプルに勾配に学習率lrだけかけた値でパラメータを更新する最急降下法
引数に学習率を受け取り，最適化関数を返す．
function SGD(lr=1e-3) f = t::Tensor -> t.x = t.x - lr * t.grad optimizer = t::Tensor -> tensor_map(f, t) return optimizer end
あとは具体的な演算を実装するだけ
モジュール化した．
https://github.com/Soyukke/julia-ayame-machine-learning/blob/master/DL.jl
Irisデータセットを使ってみた
データの読み込み
Irisデータセットを使用する．
fileはここからダウンロードしたものを使用した．
https://github.com/pandas-dev/pandas/blob/master/pandas/tests/data/iris.csv
CSVというパッケージを使用して読み込む．(openで読み込んだほうが速いと思う)
using CSV filename = "iris.csv" data = CSV.read(filename) n_data, n_type = size(data) # n_train = 150 * 0.8 = 120 n_train::Int64 = n_data * 0.8 # n_test = 150 - 120 = 30 n_test::Int64 = n_data - n_train
Irisデータセットでは，花の特徴量が4つと花の名前というデータになっている．
機械学習モデルに花の特徴量を入力して，花の名前を予測することを考える．花の種類は
unique(data[end])を実行して
["Iris-setosa", "Iris-versicolor", "Iris-virginica"]の三種類であることがわかる．
なのでモデルの出力ははonehot形式にする．
Iris-setosaを予測する場合，出力が[x1, x2, x3]となり，Iris-setosaに対応するx1が最大値を取ればよいとする．Iris-versicolorであればx2が最大，Iris-virginicaならx3が最大になればいい．
# 花の名前 -> Onehot # irisデータセットでは3種類の花の名前がある 花の名前s = unique(data[end]) # n_hana = 3 n_hana = length(花の名前s) vid2onehot = function(x::Int64) onehot = zeros(n_hana, 1) onehot[x, 1] = 1.0 return onehot end # 花の名前をonehot形式に変換する関数 花toOnehot = x::String -> vid2onehot(findfirst(花の名前s .== x)) # argmax(onehot)を花の名前に変える関数 id2name = x -> 花の名前s[x] # 入力の次元 in_features = n_type - 1 # 出力の次元 out_features = n_hana
dataをArray形式にする．
data = map(i -> [data[i, :]...], 1:n_data)
これを学習用とテスト用のデータに分割するが，まずはdataをシャッフルする．Arrayをシャッフルする関数を用意する．StatsBaseというパッケージを使用した．
using StatsBase # Arrayをシャッフルする関数 function shuffle(vx) n = length(vx) return sample(vx, n, replace=false) end
学習データとテストデータに分割して，さらに入力用のvtrain_xと教師データ用のvtrain_yに分割する．
data = shuffle(data) # 学習データ vtrain_x = map(i -> Float64.(data[i][1:in_features]), 1:n_train) vtrain_y = map(i -> data[i][n_type], 1:n_train) vtrain_y = map(花toOnehot, vtrain_y) traindataset = collect(zip(vtrain_x, vtrain_y)) # テストデータ vtest_x = map(i -> Float64.(data[i][1:in_features]), n_train+1:n_data) vtest_y = map(i -> data[i][n_type], n_train+1:n_data) vtest_y = map(花toOnehot, vtest_y) testdataset = collect(zip(vtest_x, vtest_y))
モデルの構築
機械学習モデルを構築していく
という単純なモデルを構築する．
コードではこのようにしてmodelを実装できる．
# 最適化するパラメータ W1 = VariableTensor(randn(32, in_features)/sqrt(32*in_features)) b1 = VariableTensor(zeros(32, 1)) W2 = VariableTensor(randn(out_features, 32)/sqrt(out_features*32)) b2 = VariableTensor(zeros(out_features, 1)) # 層 l1 = x -> W1 * x + b1 l2 = x -> W2 * x + b2 model = x::Tensor -> l2(ReLU(l1(x)))
もう少しおしゃれな書き方をしてみよう．
要は関数の合成をしているだけなので関数を合成する関数を実装する．
∘(f::Function, g::Function) = x::Tensor-> f(g(x::Tensor))
これを使うと...
model = l2 ∘ ReLU ∘ l1
ふつくしい．．．（∘をコピーペーストして入力しているところには目を瞑る）
損失関数
これだけでよい
loss_func = (predict::Tensor, data::Tensor) -> MSE(predict, data)
最適化関数
らくらく．学習率は1e-4にした．
optimizer = SGD(1e-4)
学習
trainのコード．詳細はコードのコメントに書いた．
# train_lossの値とテストデータの正解率のログを取る用． f = open("loss.dat", "w") # 1000 epochの学習 for epoch in 1:1000 # train @show epoch # train_lossを初期化 train_loss::Float64 = 0 # traindatasetをシャッフルする． traindataset = shuffle(traindataset) for (batch_idx, (x, y)) in enumerate(traindataset) # 入力する特徴量 x = ConstantTensor(x) # 教師用データ. onehot形式 y = ConstantTensor(y) # forward計算 predict = model(x) # lossを計算 loss = loss_func(predict, y) train_loss += loss.x # W1.grad, W2.grad, b1.grad, b2.gradに微分を足しこんでいく backward(loss) # .gradを使ってW1.x, W2.x, b1.x, b2.xを更新する # batch_size = 30で学習する. つまり，30回間隔でパラメータを更新する． if batch_idx % 30 == 0 || batch_idx == n_train # パラメータの更新 optimizer(loss) # .gradを0にする zero_grad(loss) end end @show train_loss # test 正解数::Int64 = 0 for (x, y) in testdataset # 入力 x = ConstantTensor(x) # 教師データ y = ConstantTensor(y) p = model(x) pname = id2name(argmax(p.x)) name = id2name(argmax(y.x)) 正解数 += pname == name end @printf(f, "%d, % .3e, %.3f\n", epoch, train_loss, 正解数 / n_test) end close(f) end
実行結果
まずはMSEの学習データのロスから．それっぽい見た目をしている．
続いて正解率．今回はテストデータとして30 / 150個を使っている．分類問題なのでMSEではなくCrossEntropyを使わないとだめかな？と思ったが割と問題なさそうだった．
