[Rの]Andrew先生の機械学習4[練習]

3回目後半と4回目をまとめてしまう。
3回目後半はニューラルネットワークのフィードフォワード。
4回目はニューラルネットワークの学習。

対象は前回と同様手書き数字

今回作るニューラスネットワークの構成は以下
・手書き数字の画像20*20ピクセルなので、入力素子400個
・隠れ層が1層あり、隠れ素子素子25個
・1~10の数字を識別するので、出力素子は10個

・学習させるニューラルネットワークのパラメータは、バイアス項を除くと400*25の行列と25*10の行列。
・非線形出力関数はシグモイド関数




R
学習するニューラスネットワークのパラメータはランダムに初期化しておく。 
全てが0だったりすると学習が全く進まなかったり、ローカルミニマムが大量にあるため、こうするという記憶。
randInitializeWeights <- function(L_in, L_out)
{
  epsilon_init = 0.12
　　return(matrix(runif(L_out * (1 + L_in), min = -epsilon_init, max = epsilon_init), L_out, 1 + L_in))
}

コスト関数。誤差逆伝播法で計算する勾配関数も。
nnCostFunction <- function(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda)
{
  m <- length(X[, 1])
  Theta1 <- matrix(nn_params[1:(hidden_layer_size * (input_layer_size + 1))], hidden_layer_size, input_layer_size + 1)
  Theta2 <- matrix(nn_params[(length(Theta1) + 1):length(nn_params)], num_labels, hidden_layer_size + 1)

  a1 <- X
  a1 <- cbind(matrix(1, length(a1[, 1]), 1), a1)
  z2 <- a1 %*% t(Theta1)
  a2 <- sigmoid(z2)
  a2 <- cbind(matrix(1, length(a2[, 1]), 1), a2)
  z3 <- a2 %*% t(Theta2)
  a3 <- sigmoid(z3)

  yvector <- matrix(0, length(y[, 1]), num_labels)
  for (i in 1 : num_labels)
  {
    yvector[, i] <- y[, 1] == i
  }
 
  # ----- cost function -----
  J <- sum(sum(-yvector * log(a3) - (1 - yvector) * log(1 - a3))) / m
 
  # ----- regularized cost function -----
  J = J + lambda / (2 * m) * (sum(sum(Theta1[, 2:length(Theta1[1, ])] ^ 2))+ sum(sum(Theta2[, 2:length(Theta2[1, ])] ^ 2)))

  # ----- back propagation -----
  delta3 = (a3 - yvector)
  delta2 = delta3 %*% Theta2 * sigmoidGradient(cbind(matrix(1, length(z2[, 1]), 1), z2))
  Delta1 = t(delta2[, 2:length(delta2[1, ])]) %*% a1
  Delta2 = t(delta3) %*% a2

  Theta1_grad = Delta1 / m
  Theta2_grad = Delta2 / m;

  # ----- regularized derivatives -----
  Theta1[, 1] <- 0
  Theta2[, 1] <- 0
  Theta1_grad <- Theta1_grad + lambda * Theta1 / m
  Theta2_grad <- Theta2_grad + lambda * Theta2 / m
  grad = cbind(matrix(Theta1_grad, 1, length(Theta1_grad)), matrix(Theta2_grad, 1, length(Theta2_grad)))

  return(cbind(J, grad))
}

nnCostFunction.Value <- function(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda)
{
  return(nnCostFunction(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda)[1])
}

nnCostFunction.Gradient <- function(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda)
{
  result <- nnCostFunction(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda)
  return(result[2:length(result)])
}

ニューラルネットワークのフィードフォワード
nn_predict <- function(Theta1, Theta2, X)
{
    m <- length(X[, 1])
    num_labels <- length(Theta2[, 1])
   
    a1 <- cbind(matrix(1, m, 1), X)
    a2 <- sigmoid(a1 %*% t(Theta1))
    a2 <- cbind(matrix(1, length(a2[, 1]), 1), a2)
    a3 <- sigmoid(a2 %*% t(Theta2))

    p <- matrix(0, m, 1)
    for (i in 1 : m)
    {
        p[i, 1] <- which.max(a3[i, ])
    }
    return(p)
}

メイン関数
ex4 <- function()
{
     require(R.matlab)
    datamat <- readMat('ex4data1.mat')
    X <- datamat{dollar}X
    y <- datamat{dollar}y
    displayData(X)

    input_layer_size <- 400  # 20x20 Input Images of Digits
    hidden_layer_size <- 25 # 25 hidden units
    num_labels <- 10 # 10 labels, from 1 to 10  

    # ----- randomly initialize nn parameters ----- 
    initial_Theta1 <- randInitializeWeights(input_layer_size, hidden_layer_size);
    initial_Theta2 <- randInitializeWeights(hidden_layer_size, num_labels);
    initial_nn_params <- cbind(matrix(initial_Theta1, 1, length(initial_Theta1)), matrix(initial_Theta2, 1, length(initial_Theta2)))

    # ----- train Neural Network -----
    lambda <- 1
    result <- optim(par = initial_nn_params, fn = nnCostFunction.Value, gr = nnCostFunction.Gradient, input_layer_size = input_layer_size, hidden_layer_size = hidden_layer_size, num_labels = num_labels, X = X, y = y, lambda = lambda, method= "CG")

    # ----- extract trained parameters -----
    nn_params <- result$par
    Theta1 <- matrix(nn_params[1:(hidden_layer_size * (input_layer_size + 1))], hidden_layer_size, input_layer_size + 1)
    Theta2 <- matrix(nn_params[(length(Theta1) + 1):length(nn_params)], num_labels, hidden_layer_size + 1)
   
    # ----- feed forward with NN -----
    pred <- nn_predict(Theta1, Theta2, X)
    cat("Training Set Accuracy : ", mean(pred == y) * 100, "\n")
    return(summary(result))
}



結果
初期パラメータがランダムなので結果は毎回変わる。
> ex4()
Training Set Accuracy :  95.14
            Length Class  Mode  
par         10285  -none- numeric
value           1  -none- numeric
counts          2  -none- numeric
convergence     1  -none- numeric
message         0  -none- NULL  
>

学習にそれなりに時間がかかってしまうが、、誤差逆伝播法を使わなかった場合さらに時間がかかる。
これを見るには、上のメインのoptim関数の引数のうちgr = nnCostFunction.Gradientを消すとよい。そうすると、optim関数がパラメータひとつひとつに対して数値微分を行うため、計算量が恐ろしく増大してしまう。