【Rで機械学習】ニューラルネットワークを自作してみる。

やっぱりニューラルネットワークの基礎くらいはきちんとおさえておいた方が良いだろうなぁ、ということで、取り急ぎこちらの書籍を読んでみました。

ニューラルネットワーク自作入門

ニューラルネットワーク自作入門

基礎から丁寧に解説してくれている本です。前提となる知識はほぼ不要です。行列の積と微分がちょっぴり出てきますが、それも必要な解説はすべて含まれています。これを読めばニューラルネットワークの基礎の基礎はひとまずおさえられるはず、という作りになっています。
この本の良いところは、実際にPyhtonでニューラルネットワークを作ってみよう、というところで、中に実際のコードが記述されてあります。ニューラルネットワークの基礎はそれほど難しくありませんから、Pythonのコードも非常にシンプルです。Pythonの入り口に立つのにも良いのではないかと思います。

さて、じゃあ自分でもやってみよう、と思ったのですが、さすがにここでPythonのコードを丸パクリするのはどう考えてもまずいので、本書に記載されているPythonのコードをRに置き換えてみよう、ということをやってみました。

前提

まず、上で紹介した本をすでに読んでいることが前提です。本書で解説されていることをここで改めて解説することはしません。
使用言語はRです。ver.3.4.0を使っています。特別なパッケージは不要です。

前提はこれだけです。
では、早速はじめましょう。

ニューラルネットワークを自作

大まかな方針

まず、今回は本書にならって入力層1、隠れ層1、出力層1のネットワークを作成します。後々、隠れ層の数を増やせるようにしても良いですが、今回はまず、この層構成でやることを前提とします。
入力層、隠れ層のノード数はインプットデータとそのラベルの数に左右されます。隠れ層のノード数は一般に入力層のノード数以下、出力層のノード数以上が良いとされています。インプット、アウトプットに合わせてノード数が決まるようにコードを組んでも良いですが、ひとまずは自分で値を設定するようにしました。

ニューラルネットワークにおける計算の流れは大まかに

  • 入力から出力に至る順方向の計算
  • 出力データの誤差を逆伝播させる逆方向の計算

の二つのフェーズに分けられます。
順方向の計算では重み行列とインプットデータベクトルの行列積の算出を行い、逆方向の計算ではノードの重み行列を利用した誤差の分配と誤差を使った重み更新を行います。ノードから出力される値はシグモイド関数を通します。
このように考えると、大まかに次のように分解できます。

  • 行列積を算出する部分
  • 入力された値をシグモイド関数に通す部分
  • 逆伝播させる誤差を算出する部分
  • 更新させる重みの差を算出する部分

基本的なニューラルネットワークはそこまで複雑なものではありませんから、上記のパーツを組み合わせた、比較的単純なコードで実現することが可能です。

Rにおける行列演算

ニューラルネットワークでは行列積が重要な役割を担っています。
ここでRにおける行列演算について、簡単に振り返ってみましょう。

mat1 <- matrix(1:4,ncol=2,byrow = T)
mat2 <- matrix(1:4,ncol=2,byrow = F)

> mat1
     [,1] [,2]
[1,]    1    2
[2,]    3    4
> mat2
     [,1] [,2]
[1,]    1    3
[2,]    2    4

Rでは行列をmatrix()で作ることができます。
行列積は%*%演算子で計算できます。*だと行列の各要素同士の積になるので注意。

> mat1 %*% mat2
     [,1] [,2]
[1,]    5   11
[2,]   11   25
> mat1 * mat2
     [,1] [,2]
[1,]    1    6
[2,]    6   16

ただ、今回作成するニューラルネットワークではこのような行列積は出てきません。
代わりに、行列とベクトルの積が出てきます。といってもやることは同じなんですけれど。

mat1 <- matrix(1:4,ncol=2,byrow = T)
vec <- c(1,2)

> mat1
     [,1] [,2]
[1,]    1    2
[2,]    3    4
> vec
[1] 1 2
> mat1 %*% vec
     [,1]
[1,]    5
[2,]   11

Rの場合はこれでOKです。
Pythonではベクトルを転置していますが、Rの場合、そこは空気を読んでくれているようです。

ちなみに、Rでの行列の転置はt()で行います。

> mat1
     [,1] [,2]
[1,]    1    2
[2,]    3    4
> t(mat1)
     [,1] [,2]
[1,]    1    3
[2,]    2    4

このあたりは問題ないですよね?
では、実際のコーディングに移りましょう。

初期設定

まずは初期設定です。
全体のコードの中では計算実行の少し手前で行いますが、どんな変数を使うのかを示すため、最初に記載します。

# set node number
# in this example, only 1 hidden layer is used
# in the near future, this code will be expanded to multi-hidden layers
numInput <- 3
numHidden <- 3
numOutput <- 3

# set learning rate
lr <- 0.01

# set learning epoch
epoch <- 5

最初に各層のノード数を設定しています。入力層と出力層のノード数はデータの形式に依存するため、いろんな値になりえますが、ひとまず例示としてすべてのノード数を3としています。
学習率は0.01とゆっくりにしています。大きすぎても小さすぎてもダメです。
epochは同一の入力データを使って学習を繰り返す回数です。例示では5としています。この値によっても、学習の精度が変わってきます。
このあたりの詳細は本書を読んでみてください。

次に、重みの初期値です。
本書を読んでいない方には???な話かもしれません。
ニューラルネットワークでは、入力層から隠れ層、隠れ層から出力層にいたるノードのつながりに重みをもたせています。各ノードで行っている計算は同じで、学習しているのはこの「重み」の部分です。そして、この「重み」の初期値がいい加減だと学習が正しく行われないことがわかっています。

# set initial weight-matrix
# by uniform distribution
# limited in +-1/sqrt(number of input layer)
Whi <- matrix(runif(n=numInput*numHidden,max=(1/sqrt(numInput))*2),ncol=numInput,byrow=T) - (1/sqrt(numInput))
Woh <- matrix(runif(n=numHidden*numOutput,max=(1/sqrt(numHidden))*2),ncol=numHidden,byrow=T) - (1/sqrt(numHidden))

Whiは入力層と隠れ層の間の重み、Wohは隠れ層と出力層の間の重みです。本書では一様乱数を使う方法と正規分布を使う方法の二つが紹介されていますが、今回は一様乱数を使用する方式でコーディングしました。
また、本書にある通り、重みの初期値は経験則として上限下限を設けることが推奨されています。そのため、上記のような重みの初期値設定をしています。

必要な関数定義を行う

初期設定ができたので次に必要な関数を作っていきます。
といっても、そんな面倒なものはありません。
ひとつひとつ見ていきましょう。

# sigmoid function
# in this example, logistic function is used
sigm <- function(x){
  1 / (1 + exp(-x))
}

まずはシグモイド関数です。本書にならい、シグモイド関数にはロジスティック関数を使用します。

# forward
funcForward <- function(inputVector, weightMatrix){
  outVec <- weightMatrix %*% inputVector
  sapply(outVec,sigm)
}

次に順方向の流れです。前層の各ノードに入力されたデータを各重みに応じて次層に分配し、次層のノードでシグモイド関数に通し、その値をもって次層の出力とします。この計算は入力層→隠れ層でも隠れ層→出力層でも同じです。

誤差の計算は出力層とそれより前の層で計算方法が少し違うので、関数にはしません。

# backward delta weight
funcBackward <- function(Error,output1,output2,weightMatrix){
  deltaWeight <- lr * (Error * output1 * (1 - output1)) %*% t(output2)
  out <- deltaWeight + weightMatrix
  return(out)
}

最後に誤差を利用した重みの更新です。分配された誤差と各層の出力値の積に学習率を乗じることで、重みの変化量を算出しています。この変化量を現在の重みに加えることで重みを更新していきます。

これで必要な関数定義ができたので、全体の計算の流れを作っていきましょう。

計算の全体像

作成した関数を組み合わせて「入力→出力→重み更新」のサイクルを作ります。

# define input-output cycle
funcCycle <- function(inputData, targets, Woh, Whi){
  # calc input-to-hidden and hidden-to-output
  Ohi <- funcForward(inputData, Whi)
  Ooh <- funcForward(Ohi, Woh)
  # calc error
  Eo <- targets - Ooh
  Eh <- t(Woh) %*% Eo
  # create update weight matrix
  Woh <- funcBackward(Eo,Ooh,Ohi,Woh)
  Whi <- funcBackward(Eh,Ohi,inputData,Whi)
  # return updated weight matrix
  return(list(Woh,Whi))
}

前層から次の層に渡す値の計算、誤差の計算、重みの変化量の計算まで行い、重みの変化量を戻り値としてとります。これが1サイクルで、インプットされたデータ数×epoch数だけ繰り返すようにすれば良いわけです。

このような感じに。

for(i in 1:epoch){
  for(n in 1:nrow(inputData)){
    outWeight <- funcCycle(as.numeric(inputData[n,1:3]), targets[n,], Woh, Whi)
    Woh <- outWeight[[1]]
    Whi <- outWeight[[2]]
  }
}

はい、これで計算部分は出来上がりました。

実験してみる

インプットを作る際の注意点

本書にも注意として書かれていることなのですが、インプットデータは0より大きい1までの値となっている必要があります。
というののも、隠れ層、出力層から出力される値はシグモイド関数を通しているので、少なくとも0以上1以下の値をとります。それに対して入力層の値だけその範囲を超えた値になってしまうと、入力層と隠れ層の間の重み更新の変化量がおかしなことになってしまいます。
また、ノードから出力される値が0をとってしまうと、重みの更新ができなくなってしまいます。重みの変化量に出力の値を使うため。試しに計算してみるとわかるでしょう。

以上のような事情を踏まえて、簡単なトイタスクを作ってみましょう。

※本書ではMINSTのデータを使っていますが、ここではもっとずっと単純な例を使います。

インプットデータの作成

今回はインプットの量もアウトプットの量も少ない、簡単な例で試してみます。
インプットデータ生成のコードを先に示します。

inputData <- lapply(seq(1,100000),function(i){
  x <- floor(runif(1,min=1,max=4))
  outVec <- c(rnorm(n=3,x*10,sd=3),x)
  return(outVec)
})
inputData <- do.call(rbind,inputData)
for(i in 1:3){
  inputData[,i] <- inputData[,i] - min(inputData[,i])
  inputData[,i] <- (inputData[,i]/(max(inputData[,i])) * 0.99) + 0.01
}

まず、1,2,3のいずれかの整数を一様乱数で発生させます。それをxとします。
次に、平均10x、標準偏差3の正規分布に従う乱数を3つ発生させます。それを100,000回繰り返します。
発生させた乱数はそれぞれ3次元空間上のx,y,z座標に相当し、はじめに発生させた一様乱数xはその三次元空間上の点がどのグループに属しているかを示します。
そして、得られた座標を注意点に従って0より大きく1までの間の値に変換します。
これを3次元空間上にグループごとに色分けしてプロットするとこのような感じになります。

f:id:wanko_sato:20171001133439p:plain

三色団子ができました。
で、この座標をインプットとして、その点が黒、赤、緑どのグループに属するかをニューラルネットワークに学習させます。

したがって、入力層のノード数は3、出力層のノード数も3になります。隠れ層のノード数は入力層以下、出力層以上が望ましいので、3となります。
最初に設定したノード数はこのためのものなのでした。

学習用のラベルの作成

本書にもある通り、学習用のラベルにも少し工夫が必要です。

targets <- lapply(inputData[,4],function(x){
  out <- rep(0.01,numOutput)
  out[as.numeric(x)] <- 0.99
  return(out)
})
targets <- do.call(rbind,targets)

何をやっているかというと、出力層の各ノードに対応する3つの値を作り、ノード1をグループ1である確率、ノード2をグループ2である確率、ノード3をグループ3である確率とみなし、入力データがグループ1であれば(0.99,0.01,0.01)というベクトルを作る、というようなことです。

> head(targets)
     [,1] [,2] [,3]
[1,] 0.01 0.01 0.99
[2,] 0.01 0.99 0.01
[3,] 0.01 0.01 0.99
[4,] 0.99 0.01 0.01
[5,] 0.01 0.01 0.99
[6,] 0.01 0.01 0.99

このような感じのデータになっています。ニューラルネットワークは出力がこの値に近づいていくよう、重みを調整していく作業を中で繰り返しているわけです。
ここで値の上限を0.99、下限を0.01にしているのにはわけがあります。各ノードから出力される値はシグモイド関数を通ってきているのでした。シグモイド関数は限りなく0に近づくことができますが0になるのは無限遠です。1についても同様です。0または1を目標値としてしまうと出力値が無限にならなければならず、そのような値を出力するように重みを調整してしまうとおかしなことが起こってしまいます。したがって、許容できる範囲として、上限0.99、下限0.01としているわけです。
このあたりの詳しいことは本書を参考になさってください。

評価用のデータの作成

と、ここまでインプットデータのことばかり書いてきましたが、ニューラルネットワークが適切に学習できているかを評価するためのデータが必要です。

testData <- lapply(seq(1,300),function(i){
  x <- floor(runif(1,min=1,max=4))
  outVec <- c(rnorm(n=3,x*10,sd=3),x)
  return(outVec)
})
testData <- do.call(rbind,testData)
for(i in 1:3){
  testData[,i] <- testData[,i] - min(testData[,i])
  testData[,i] <- (testData[,i]/(max(testData[,i])) * 0.99) + 0.01
}

インプットデータと同じ方法で300個の評価用データを作りました。学習精度はこのデータを使って評価します。

学習させてみよう

以上で準備が整ったので、実際に学習させてみましょう。
関数定義、データの準備は出来ているので、初期設定と学習の実行部分を書きます。

#####################################################################################
# initialize

# set node number
# in this example, only 1 hidden layer is used
# in the near future, this code will be expanded to multi-hidden layers
numInput <- 3
numHidden <- 3
numOutput <- 3

# set learning rate
lr <- 0.01

# set learning epoch
epoch <- 5

# set initial weight-matrix
# by uniform distribution
# limited in +-1/sqrt(number of input layer)
Whi <- matrix(runif(n=numInput*numHidden,max=(1/sqrt(numInput))*2),ncol=numInput,byrow=T) - (1/sqrt(numInput))
Woh <- matrix(runif(n=numHidden*numOutput,max=(1/sqrt(numHidden))*2),ncol=numHidden,byrow=T) - (1/sqrt(numHidden))

#####################################################################################
# execute cycle
for(i in 1:epoch){
  for(n in 1:nrow(inputData)){
    outWeight <- funcCycle(as.numeric(inputData[n,1:3]), targets[n,], Woh, Whi)
    Woh <- outWeight[[1]]
    Whi <- outWeight[[2]]
  }
}

ちょっぴり時間がかかりますので少し待ちましょう。
学習が完了したら、次のコードで評価します。

#####################################################################################
# evaluate error rate
testOut <- lapply(seq(1,nrow(testData)),function(x){
  funcForward(funcForward(testData[x,1:3],Whi),Woh)
})
testOut <- do.call(rbind,testOut)
testTag <- sapply(seq(1,nrow(testOut)),function(x){
  which.max(testOut[x,1:3])
})

テストデータを入力として学習済みのニューラルネットワークに通し、どんな出力が得られるかを確認します。

> head(testOut)
            [,1]      [,2]        [,3]
[1,] 0.008558283 0.4132223 0.875065271
[2,] 0.014413098 0.4018476 0.800517026
[3,] 0.085244262 0.3700818 0.406289701
[4,] 0.067191789 0.3720620 0.457506149
[5,] 0.007589988 0.4153517 0.887068359
[6,] 0.897187592 0.2914357 0.007308831

こような感じになりました。すでに説明したように、各列が各グループへの適合確率です。一行目ですと、グループ1である確率が0.8%、グループ2である確率が41%、グループ3である確率が87%であることを示しています。1行目の場合は3列目が一番確率が高いので答えはグループ3、とします。それを全テストデータに対して行ったものがtestTagになります。
testTagを実際のグループと比較し、正解率を算出します。

# calc error rate
eRate <- sapply(seq(1:300),function(x){
  if(testData[x,4]==testTag[x]){
    out <- 1
  }else{
    out <- 0
  }
  return(out)
})
eRate <- sum(eRate)/300

とすると、

> eRate
[1] 0.9433333

正答率94.3%とでました。非常に単純なコードで、かつ単純なデータですが、まぁほどほどなのではないでしょうか。

実は紆余曲折あった

インプットデータの量の問題

と、ここまで順風満帆にきたように書いていますが、ちょっとだけつまずいたところがありました。「インプットデータの量」です。
本書中では、MINSTのデータを用い、100のインプットでもそこそこの精度がでる、というような記述がありました。学習精度に影響を及ぼすのは学習率であったりepochであったりする、と。そのなかでは「インプットデータの量」については特に明確な記載はなかったのです。

ですが、今回の簡単なデータで試してみて、インプットデータの量も学習精度に大きく影響を及ぼしうるものだということを再認識しました。

f:id:wanko_sato:20171001140941p:plain

インプットデータの量と学習精度の関連を示したグラフです。
一番最初はインプットデータ3000個でテストを行いました。その際の学習精度は70%程度でほとんど使い物にならないレベルです。

f:id:wanko_sato:20171001141131p:plain

団子の真ん中の、本来赤になるべきプロットがすべて緑になっており、グループ2として認識されていません。

f:id:wanko_sato:20171001141303p:plain

データ数を100倍の300,000個に増やすとようやくグループ2の半分くらいがグループ2として認識されるようになり、

f:id:wanko_sato:20171001141436p:plain

1,000,000個でグループ2の3/4くらいがちゃんと認識されるようになりました。
これでようやく学習精度96%程度です。正直かなり微妙な線だな・・・と。

インプットデータの質の問題

果たしてこれって量だけの問題なのかな?とも思うのです。
というのも、インプットされるデータの種類がこの場合はx,y,zの3種類なわけで、そこから得られる特徴って、特にグループの境界あたりでは微妙だったりするんじゃなかろうかな、と。
逆に、MINSTみたいにインプットが784ピクセルの画像だったりすると入力されるデータが非常に豊富なわけで、特徴が抽出しやすく、少ないトレーニングデータでもそこそこの精度が達成できているんじゃなかろうかな、と。

そのあたりの実態ははっきりとはわからないですが、ニューラルネットワークに向いているデータと向いていないデータがあるのじゃないのかな、と考えてみたりするわけです。

コード全体

というわけで、最後に今回作成したコード全体を掲載して終わりにしたいと思います。

#####################################################################################
# define function

# sigmoid function
# in this example, logistic function is used
sigm <- function(x){
  1 / (1 + exp(-x))
}

# forward
funcForward <- function(inputVector, weightMatrix){
  outVec <- weightMatrix %*% inputVector
  sapply(outVec,sigm)
}

# backward delta weight
funcBackward <- function(Error,output1,output2,weightMatrix){
  deltaWeight <- lr * (Error * output1 * (1 - output1)) %*% t(output2)
  out <- deltaWeight + weightMatrix
  return(out)
}

# define input-output cycle
funcCycle <- function(inputData, targets, Woh, Whi){
  # calc input-to-hidden and hidden-to-output
  Ohi <- funcForward(inputData, Whi)
  Ooh <- funcForward(Ohi, Woh)
  # calc error
  Eo <- targets - Ooh
  Eh <- t(Woh) %*% Eo
  # create update weight matrix
  Woh <- funcBackward(Eo,Ooh,Ohi,Woh)
  Whi <- funcBackward(Eh,Ohi,inputData,Whi)
  # return updated weight matrix
  return(list(Woh,Whi))
}

#####################################################################################
# create inputData
inputData <- lapply(seq(1,100000),function(i){
  x <- floor(runif(1,min=1,max=4))
  outVec <- c(rnorm(n=3,x*10,sd=3),x)
  return(outVec)
})
inputData <- do.call(rbind,inputData)
for(i in 1:3){
  inputData[,i] <- inputData[,i] - min(inputData[,i])
  inputData[,i] <- (inputData[,i]/(max(inputData[,i])) * 0.99) + 0.01
}

targets <- lapply(inputData[,4],function(x){
  out <- rep(0.01,numOutput)
  out[as.numeric(x)] <- 0.99
  return(out)
})
targets <- do.call(rbind,targets)

#####################################################################################
# create inputData
testData <- lapply(seq(1,300),function(i){
  x <- floor(runif(1,min=1,max=4))
  outVec <- c(rnorm(n=3,x*10,sd=3),x)
  return(outVec)
})
testData <- do.call(rbind,testData)
for(i in 1:3){
  testData[,i] <- testData[,i] - min(testData[,i])
  testData[,i] <- (testData[,i]/(max(testData[,i])) * 0.99) + 0.01
}

#####################################################################################
# initialize

# set node number
# in this example, only 1 hidden layer is used
# in the near future, this code will be expanded to multi-hidden layers
numInput <- 3
numHidden <- 3
numOutput <- 3

# set learning rate
lr <- 0.01

# set learning epoch
epoch <- 5

# set initial weight-matrix
# by uniform distribution
# limited in +-1/sqrt(number of input layer)
Whi <- matrix(runif(n=numInput*numHidden,max=(1/sqrt(numInput))*2),ncol=numInput,byrow=T) - (1/sqrt(numInput))
Woh <- matrix(runif(n=numHidden*numOutput,max=(1/sqrt(numHidden))*2),ncol=numHidden,byrow=T) - (1/sqrt(numHidden))

#####################################################################################
# execute cycle
for(i in 1:epoch){
  for(n in 1:nrow(inputData)){
    outWeight <- funcCycle(as.numeric(inputData[n,1:3]), targets[n,], Woh, Whi)
    Woh <- outWeight[[1]]
    Whi <- outWeight[[2]]
  }
}

#####################################################################################
# evaluate error rate
testOut <- lapply(seq(1,nrow(testData)),function(x){
  funcForward(funcForward(testData[x,1:3],Whi),Woh)
})
testOut <- do.call(rbind,testOut)
testTag <- sapply(seq(1,nrow(testOut)),function(x){
  which.max(testOut[x,1:3])
})

# calc error rate
eRate <- sapply(seq(1:300),function(x){
  if(testData[x,4]==testTag[x]){
    out <- 1
  }else{
    out <- 0
  }
  return(out)
})
eRate <- sum(eRate)/300

いやぁ、こうして改めてみるとずいぶんコンパクトなコードでニューラルネットワークできるんですねぇ。