RでDeep Learningの一種をやってみる

このブログのTips052で、RでDeep Learningをやっているのを紹介してもらったので、自分も試してみました。

「Deep Learningすげぇ！！」という話は良く聞くのですが、亜種がいっぱいあるみたいで、まだあまり調査しきれてません。また時間があれば調査してまとめられると良いのですが。

以下、RでDA(Denoising Autoencoders)を実行するプログラムです。

sigmoid <- function(x){
     return (1 / (1 + exp(-x)))
}

dA <- setRefClass(
        Class="dA",
        fields=list(input="matrix", n_visible="numeric",
                    n_hidden="numeric", W="matrix",
                    W.prime="matrix", hbias="vector",
                    vbias="vector", rng="numeric"),
        methods=list(
          get_corrupted_input=function(input, corruption_level){
            stopifnot(corruption_level < 1)  # アサーション
            return(rbinom(size=1, n=ncol(input) * nrow(input), 
                          prob=1 - corruption_level) * input)  # size=1とすることで２項分布からベルヌーイ分布を得る
          },
         
        # Encode
        get_hidden_values=function(input){
            # print("-------------------get_hidden_values()")
            # print("input=")
            # print(input)
            # print(".self$W=")
            # print(.self$W)
            # print("-------------------")
            return(sigmoid(input %*% .self$W + .self$hbias))
        },
        
        # Decode
        get_reconstructed_input=function(hidden){
            # print("-------------------get_reconstructed_input()")
            # print("hidden=")
            # print(hidden)
            # print(".self$W.prime=")
            # print(.self$W.prime)
            # print("-------------------")
            .self$W.prime <- t(.self$W)
            return(sigmoid(hidden %*% .self$W.prime + .self$vbias))
        },
        
        train=function(lr=0.1, corruption_level=0.3, input=NULL){
            if(is.null(input)==FALSE){
                .self$input <- input
            }
            input <- .self$input
            # print("-------------------train()")
            # print("input=")
            # print(input)
            # print("-------------------")
            tilde_x <- .self$get_corrupted_input(input, corruption_level)
            y <- .self$get_hidden_values(tilde_x)
            z <- .self$get_reconstructed_input(y)
            
            L_h2 <- input - z
            L_h1 <- (L_h2 %*% .self$W) * y * (1 - y)
            
            L_vbias <- L_h2
            L_hbias <- L_h1
            L_W <- (t(tilde_x) %*% L_h1) + (t(L_h2) %*% y)
            
            .self$W <- .self$W + lr * L_W
            .self$hbias <- .self$hbias + lr * colMeans(L_hbias)
            .self$vbias <- .self$vbias + lr * colMeans(L_vbias)
            # print("=============================")
            # print(L_hbias)
            # print(L_vbias)
            # print(colMeans(L_hbias))
            # print(colMeans(L_vbias))
            # print("=============================")
            
        },
        
        negative_log_likelihood=function(corruption_level=0.3){
            tilde_x <- .self$get_corrupted_input(.self$input, corruption_level)
            y <- .self$get_hidden_values(tilde_x)
            z <- .self$get_reconstructed_input(y)
            # print("=============================")
            # print(tilde_x)
            # print(y)
            # print(z)
            # print(log(z))
            # print(log(1-z))
            # print("=============================")
            
            
            cross_entropy = - mean(rowSums(.self$input * log(z) + (1 - .self$input) * log(1 - z)))
            
            return(cross_entropy)
        },
        
        reconstruct=function(x){
            y <- .self$get_hidden_values(x)
            z <- .self$get_reconstructed_input(y)
            return(z)
        }, 
        
        # コンストラクタ
        initialize=function(input=NULL, n_visible=NULL, 
                            n_hidden=NULL, W=NULL, W.prime=NULL, 
                            hbias=NULL, vbias=NULL, rng=NULL){
            if(is.null(input) == TRUE){
              input <<- matrix(NA)
            }
            else{
              input <<- input
            }
            
            if(is.null(n_visible) == TRUE){
              n_visible <- 2
              n_visible <<- n_visible
            }
            else{
              n_visible <<- n_visible
            }
            
            if(is.null(n_hidden) == TRUE){
              n_hidden <- 3
              n_hidden <<- n_hidden
            }
            else{
              n_hidden <<- n_hidden
            }
            
            if(is.null(rng) == TRUE){
              rng <- 1234
              rng <<- rng
            }
            else{
              rng <<- rng
            }
            
            if(is.null(W) == TRUE){
              a <- 1 / n_visible
              set.seed(rng)
              W <-  matrix(runif(n_visible * n_hidden, min=-a, max=a), 
                           n_visible, n_hidden)
              W <<- W
              W.prime <<- t(W)
            }
            else{
              W <<- W
              W.prime <<- t(W)
            }
            
            if(is.null(hbias) == TRUE){
              hbias <<- rep(0, n_hidden)  # initialize h bias 0
            }
            else{
              hbias <<- hbias
            }
            
            if(is.null(vbias)==TRUE){
              vbias <<- rep(0, n_visible)  # initialize v bias 0
            }
            else{
              vbias <<- vbias
            }
        }
    )
)



test_dA <- function(learning_rate=0.1, corruption_level=0.3, 
                    training_epochs=50){
  data <- rbind(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c(1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c(1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c(0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
                c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1),
                c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1),
                c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1),
                c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0))
  rng <- 123  # seed
  
  # construct dA
  da <- dA$new(input=data, n_visible=20, n_hidden=5, rng=rng)
  
  # train
  for(epoch in seq(training_epochs)){
      da$train(lr=learning_rate, corruption_level=corruption_level)
      cost = da$negative_log_likelihood(corruption_level=corruption_level)
      print(paste('Training epoch ', epoch, '  cost is ', cost))
      learning_rate <- learning_rate * 0.95
  }
  
  # test
  x <- rbind(c(1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1),
             c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0))
  
  print(da$reconstruct(x))
}

test_dA()

とりあえず動くところまでは確認（@konsonsanありがとうございます！）。

今回はダミーデータに対してやってるので、またいつかコードレビューをしつつ、サンプルデータとか実データでもやってみたい。