--- author: "Alessandro Guida, Ph.D." date: "9/16/2017" output: html_document: code_folding: show fig_height: 8 fig_width: 10 highlight: haddock theme: readable toc: no html_notebook: default version: 0.0.1 --- # Simple Keras in R The MNIST dataset is a famous dataset used for computer vision, it consists of handwritten images and has served as the benchmarking reference for training classification algorithms. **It is perfect to showcase neural network techniques.** **The dataset:** The dataset consists of n=42000 images 28x28 px resolution of digits from 0 to 9. It is provided as a flat 42000x785 matrix, where 784 columns are the flatten image (748 = 28X28) + 1 column with the label (digit number shown in the image). Each row component contains a value between one and zero and this describes the intensity of each pixel **The goal:** is to correctly identify the digits from the handwritten images. In this R notebook, we will try 3 differente neural networks (NN): * Model1: a single hidden layer perceptron; * Model2: a multilayer perceptron; * Model3: a convolutional network (CNN); ## Set up In this section we will set up the work environment by * loading the libraries and dependencies; * define the helper functions used in the analysis. ```{r, message=FALSE, warning=FALSE} devtools::install_github("rstudio/tensorflow") library(tensorflow) devtools::install_github("rstudio/keras") library(keras) library(readr) library(knitr) library(magrittr) library(purrr) library(ggplot2) library(plotly) library(dplyr) set.seed(1) # Helper Functions # ============================================================================== #' mutate_cond #' #' My dplyr extention, that allows to perform a mutate condition on a subset #' of the dataframe and return the full dataframe instead of only the modified #' subset #' #' @.data the input dataset #' @condition the condition that needs to be met to filter part of the dataset #' @... the extra argvs #' #' @return the modified dataframe mutate_cond <- function(.data, condition, ..., envir = parent.frame()) { condition <- eval(substitute(condition), .data, envir) .data[condition, ] <- .data[condition, ] %>% mutate(...) .data } #' rotate_90right #' #' Function that takes a matrix as in input and rotates it 90° on the right #' before returning it #' #' @param mymatrix a matrix #' @return a rotated matrix rotate_90right <- function(mymatrix){ # reverse and transpose apply(mymatrix , 2, rev) %>% t() } #' plot_mnist_image #' #' takes a raw from the mnist datset and renders it as an image #' #' @param dataset data.frame with 784 column #' @param row_number row of the dataframe to be printed #' #' @return it retuns an image plot_mnist_image <- function(dataset, row_number, ...){ dataset[row_number,] %>% unlist %>% matrix(ncol = 28, byrow = TRUE) %>% # convert in a 28x28 matrix rotate_90right %>% # rotate the image image(col = grey.colors(255), ...) # plot it } #' plot_mnist_images #' #' Function wrapper of plot_mnist_image() that adds functionally for #' plotting multiple images and arranging them in a grid. ncols #' specifies how many columns in the grid. #' #' @param dataset data.frame with 784 column #' @param image_row_number row of the dataframe to be printed #' @param ncols #' @param image_labels #' #' @return plot a grid of images #' #' @example plot_mnist_images(train_x, 1:5, ncols=3, image_labels = train_y[1:5]) plot_mnist_images <- function(dataset, image_row_number, ncols=1 , image_labels = FALSE){ # some input checks if (!is.vector(image_row_number)) stop("rows should be a vectors") if(!all(purrr::map_lgl(image_row_number, is.numeric))) stop("rows should be a vectors of integers") if(ncols > length(image_row_number)){ stop("you have specified more columns than the number of images to be printed") } if(length(image_labels) != 0){ if((is.vector(image_labels) || is.array(image_labels)) == FALSE) { stop("Image labels shuld be a vector or an array") } if(length(image_labels) != length(image_row_number) && image_labels != FALSE){ stop("image_labels should have the same length as image_row_number") } } # calculate number of raws required in the grid,based on the number of # columns requested. nrows = length(image_row_number) / ncols # add an extra row if there are remainders if(length(image_row_number) %% ncols != 0) nrows = nrows + 1 # set grid configuration par(mfrow=c(nrows,ncols)) for(i_row in 1:length(image_row_number)){ plot_mnist_image(dataset , image_row_number[i_row] , main=ifelse(image_labels, image_labels[i_row], "") ) } } # Implementation of a special for a "not in" '%notin%' <- function(x, table) match(x, table, nomatch = 0L) == 0L ``` ## Data Exploration and pre-processing Before diving into the building a NN, we need to do some data cleaning, and data exploration. Let's begin by loading the data and checking the dimensions of the datasets. The data were downloaded from Kaggle, which provides a traing dataset on which to build the model, and a test dataset, used by kaggle to score the model performance. The test dataset is not labelled (number of columns i 784) to prevent people from cheating. This also means that we will not be able to use it to validate the performance of our model. We will need to split the trainig set in two parts. ```{r, message=FALSE, warning=FALSE} # read input # ============================================================================== train <- readr::read_csv("../input/train.csv") test <- readr:::read_csv("../input/test.csv") # preview input dimensions data.frame( training_set = dim(train) , test_set = dim(test) , row.names = c("rows", "columns") # add rownames ) %>% knitr::kable() ``` The training set as a 1 extra columns, with the outcome variable. We need to do some pre-processing on our training set and split the regressors from the predictors. In fact, we we look at the first six column names we can see that the first column is the "label" (out outcome variable in this case). ```{r} colnames(train) %>% head() ``` Let's split the training dataset in regressors and predictors ```{r} train_y <- train$label # predictor variable train_x <- dplyr::select(train, -label) # regressors ``` Let's now visualize the label balance in our training set. We want to know whether that are some labels that are overrappresented. ```{r, fig.height=4} train_y %>% table() %>% # count how many labels in each cathegory data.frame(label=names(.)) %>% # convert to dataframe ggplot(aes(x= label, y= Freq )) + geom_bar(stat = "identity") ``` and perform a test ```{r} # calculate observed frequency in the labels freq_labels <- train_y %>% table() # calculate expected frequency expected_freq <- sum(freq_labels) / length(freq_labels) # chi squared test data.frame(freq= as.numeric(freq_labels) , expected_freq= expected_freq) %>% chisq.test() ``` According to the test the classes are unbalanced. While, this is something I will keep into consideration in case of problems, we have got so many samples that I do not expect this to be an issue. For now I will just mark it as a potential issue that might need to be explored further. The pixel values are gray scale between 0 and 255. It is good practice to normalize the values when using neural network models. We will run a simple normalization and convert the pixel values to a range from 0 to 1, by dividing each value by 255. ```{r} # normalize gray channel to a 0 to 1 scale train_x = train_x / 255 test = test / 255 ``` Let's now visualize some of the images. Each image is a 28x28 (= 784 px). Each row in the dataset is an image. We can therefore recompose the 784 length vector into a 28 rows x 28 columns matrix of color intensities and plot it as an image. Using the helper function previously defined we can show a few inputs here. ```{r} # Plot some Images (see helper functions for the code) plot_mnist_images(train_x, 1:6, image_labels = train_y[1:6], ncols=3) ``` The first image corresponds to the following data matrix ```{r} train_x[1,] %>% round(digits = 2) %>% # make digits shorter matrix(ncol = 28) %>% # convert to matrix knitr::kable() # generate kable ``` ### Are all the dimentions important? By looking at the pictures, it is clear that the corner and sides are usually absent of signal. Let's first identify those predictors which have no variance and therefore no contribution for the model. ```{r} # Image variancein the training set apply(train_x, 2, var) %>% # calculate variance for each predicto matrix(ncol=28, byrow=TRUE) %>% # recompose in a matrix image(main="Heatmap of the predictor'variance in the train dataset") ``` It would be interesting the explore and to get an idea of how many dimensions in the dataset contain the majority of the signal. In order to perform a quick check, we can run a simple PCA and define how many principal components capture the majority (90%) of the variance. In the MINST datasets, all the predictors are measured on the same scale and have the same unit. We will therefore NOT perform any scaling of the variables, especially because we are interested investigating the variance contribution of each single regressor. ```{r} # Scale: indicates wether the data should be scaled with a standard deviation of 1, before performing PCA. train_pca <- prcomp(train_x, scale = FALSE) ``` Now, we will calculate and plot the proportion of the variance explained by each principal component in a interactive plot (use mouse-over to see details about the plot). ```{r} # getting proportion of variance fro a scree plot pr.var <- train_pca$sdev^2 # calculate the proportion of variance explained by each principal component pve <- pr.var / sum(pr.var) # allow decimal representation up to 10 digits options(scipen = 10) # Generate Plot # ----------------------------------------------------------------------------- p <- pve %>% round(digits = 5) %>% # convert to dataframe data.frame( predictor_n = 1:length(.) , proportion_explained = . , cumulative = cumsum(.) ) %>% # reshape the data for ggplot tidyr::gather(c(proportion_explained, cumulative) , key = "type" , value = "variance") %>% # init ggplot ggplot(aes(y=variance, x=predictor_n, fill = type)) + geom_point(aes(color=type)) + xlab("Principal Components") + ylab(label="variance explained") + theme(legend.position="bottom") # make it interactive ggp <- ggplotly(p, width = 800, height = 400) layout(ggp, legend=list(x=0.65,y=0.5)) ``` 90 % of the variance can be explained with ```r min(which(cumsum(pve) >= 0.9))``` principal components. ### Dealing with multiple-classes Since this is a multi classification problem, when using neural networks, it is good practice to perform **One Hot Encoding (OHE)**. This step allows to make sure that you transform your outcome variable from a vector that contains values for each class value to a matrix with a Boolean for each class value and whether or not a given instance has that class value or not. OHE can be down with the ```to_categorical()``` function from the keras package. ```{r} # One hot encode training target values train_y_mtx <- keras::to_categorical(train_y) # Preview head(train_y_mtx) ``` ### Split training set into training and validation We now want to split our training dataset into two parts. One will be used for training the NN and the other one will be used for validation. ```{r} # global variable for the training set split % TRAININGSET_SPLIT = 0.1 # random sampling from the trainign set val_idx <- sample(1:nrow(train_x) , round(nrow(train_x) * TRAININGSET_SPLIT) , replace = FALSE) # remove validation portion from trainin dataset train_x_sub <- train_x[-val_idx,] train_y_mtx_sub <- train_y_mtx[-val_idx,] # create validation dataset val_x <- train_x[val_idx,] %>% as.matrix() val_y <- train_y[val_idx] ``` Let's make sure we have not introduced some labeling bias with our random sampling by checking the label balance in the training and validation set. ```{r, fig.height=4} train_y %>% data.frame(label=.) %>% # conver to dataframe dplyr::mutate(addindex = 1:length(train_y)) %>% # add index to df dplyr::mutate(dataset = NA) %>% # add dataset variable mutate_cond(addindex %in% val_idx, dataset = "validation") %>% mutate_cond(addindex %notin% val_idx, dataset = "training") %>% dplyr::select(-addindex) %>% # remove col dplyr::group_by(label, dataset) %>% dplyr::summarise(count=n()) %>% # count ggplot(aes(x= label, y= count, fill=dataset)) + # generate plot geom_bar(stat = "identity") + ggtitle("label balance after dataset split") ``` ## MODEL 1: Construct a very simple perceptron model We will build a simple network topology, a single-layer perceptron. This type of neural network is fully connected. We will use the ReLU activation function. The softmax activation function is used in the output layer to make sure that the output values are in the range of 0 and 1 and may be used as predicted probabilities. ```{r} # Initialize a sequential model model <- keras_model_sequential() # Add layers to the model model %>% layer_dense(units = 400, activation = 'relu', input_shape = ncol(train_x)) %>% layer_dense(units = 10, activation = 'softmax') # Preeview the model topology summary(model) ``` ### Compile the model Now that you have set up the architecture of the model, let's compile it with the ```adam optimizer``` and the ```categorical_crossentropy``` loss function. We will monitor the accuracy during the training by passing **'accuracy'** to the metrics argument. ```{r} # Compile the model model %>% compile( loss = 'categorical_crossentropy', optimizer = 'adam', metrics = 'accuracy' ) ``` ### Fit the model In this sections we will define some of the parameters required for fitting the model. We chose to set the number of epochs to 10 (number of times our training set is learned through the network). ```{r} # We need to convert the input in a matrix or Keras will generate errors if we # input a data.frame train_x_sub = as.matrix(train_x_sub) # Fit the model history <- model %>% fit(train_x_sub, train_y_mtx_sub, epochs = 20, batch_size = 50, validation_split=0.2) ``` ### Preview model performance The following plot shows the performance of the model, the top part is the accuracy of the model, the bottom part is the loss function over the number of epochs. ```{r} # Plot the history plot(history) ``` ```{r, eval=FALSE, include=FALSE} # Plot the model loss of the training data plot(history$metrics$loss, main="Model Loss", xlab = "epoch", ylab="loss", col="blue", type="l") #ylim=c(0,1) # Plot the model loss of the test data lines(history$metrics$val_loss, col="green") # Add legend legend("topright", c("train","test"), col=c("blue", "green"), lty=c(1,1)) ``` ```{r , eval=FALSE, include=FALSE} # Plot the accuracy of the training data plot(history$metrics$acc, main="Model Accuracy", xlab = "epoch", ylab="accuracy", col="blue", type="l", ylim=c(0,1)) # Plot the accuracy of the validation data lines(history$metrics$val_acc, col="green") # Add Legend legend("bottomright", c("train","test"), col=c("blue", "green"), lty=c(1,1)) ``` ## Validation Let's now run the model on the validation dataset we prepared previously. This dataset was never see by the NN we have trained so it will represented novel data. We generate a prediction using our model ```{r} # Predict the classes for the test data val_y_pred <- model %>% predict_classes(val_x, batch_size = 128) ``` and plot a contingency table between the expected labels and the predicted labels. ```{r} # Plot contingegy table as heatmap plot table(val_y_pred, val_y) %>% reshape::melt() %>% ggplot(aes(x=val_y, y=val_y_pred)) + geom_tile(aes(fill = value)) + scale_x_continuous("Truth", breaks = 0:9) + scale_y_continuous("Prediction", breaks = 0:9) + geom_text(aes(label = value), color="white") + scale_size(range = c(15,50)) ``` Prediction accuracy: ```r prop.table(table(val_y_pred == val_y))[2]``` ## MODEL2: Construct a multilayer perceptron model In this section we will try to build a more complext model, adding some extra hidden layers. To prevent overfitting, we will introduce the dropout layer. Dropout consists of randomly setting a fraction rate of input units to 0 at each update during training time. It helps preventing over-fitting because it forces the network to explore new paths. * [link to dropout paper](http://www.cs.toronto.edu/~rsalakhu/papers/srivastava14a.pdf) ```{r} # Define topology # ------------------------------------------------------------------------------ # Initialize a sequential model model2 <- keras_model_sequential() # Add layers to the model model2 %>% layer_dense(units = 784, activation = 'relu', input_shape = ncol(train_x)) %>% layer_dropout(0.25) %>% layer_dense(units = 400, activation = 'relu') %>% layer_dropout(0.25) %>% layer_dense(units = 200, activation = 'relu') %>% layer_dropout(0.25) %>% layer_dense(units = 10, activation = 'softmax') # Compile topology # ------------------------------------------------------------------------------ model2 %>% compile( loss = 'categorical_crossentropy', optimizer = 'adam', metrics = 'accuracy' ) # Fit topology # ------------------------------------------------------------------------------ # Fit the model history2 <- model2 %>% fit(train_x_sub, train_y_mtx_sub, epochs = 20, batch_size = 100, validation_split = 0.2) ``` ### Model performance plot ```{r} # Plot the model loss of the training data plot(history2) ``` ## Validation ```{r} # Predict the classes for the test data val_y_pred <- model2 %>% predict_classes(val_x, batch_size = 128) # Plot contingegy table as heatmap plot table(val_y_pred, val_y) %>% reshape::melt() %>% ggplot(aes(x=val_y, y=val_y_pred)) + geom_tile(aes(fill = value)) + scale_x_continuous("Truth", breaks = 0:9) + scale_y_continuous("Prediction", breaks = 0:9) + geom_text(aes(label = value), color="white") + scale_size(range = c(15,50)) ``` Prediction accuracy: ```r prop.table(table(val_y_pred == val_y))[2]``` Adding more layers and increasing the complexity of the network doesn't necessarily improves the performance of the network. More fine tuning is needed to raise the accuracy of the model. ## MODEL3: Construct a convolutional network (CNN) One problem with the application of perceptrons on image analysis, is that by flattening the image, we might be losing some interesting patters. Convolutional networks were specifically design to work with images. Let's now try to build a convolutional neural network. In this section we will create a simple CNN for MNIST that demonstrates how to use all of the aspects of a modern CNN implementation, including convolutional layers, pooling layers and dropout layers. We will now need to reshape the input data so that we can work with a convolutional network. Since the Mnist dataset is just in Gray, we will have only one extra dimension. If it was RGB, it would be 3 dimensions (red, green, blue). ```{r} # we are reshaping the input train_x_sub_cnn <- train_x_sub dim(train_x_sub_cnn) <- c(nrow(train_x_sub_cnn), 28, 28, 1) ``` Define the CNN topology: * The first hidden layer is a convolutional layer called a Convolution2D. We will start with a 32 filters with shape 5x5 and a ReLU activation function. This is the input layer, the input images should have the structure: [pixels][width][height] * Next we define a pooling layer that takes the max called MaxPooling2D. This layers collects the results from the filters in a pool size of 2×2. * Then we have regularization layer using dropout called Dropout. As we mentioned before this layer tries to prevents overfitting. We program it here to exclude 20% of neurons in the layer. * THe 2D matrix data is now converted into a vector by using a "Flatten" layer. It allows the output to be processed by standard fully connected layers. * Following a fully connected layer with 128 neurons and rectifier activation function. * Finally, the output layer has 10 neurons for the 10 classes and a softmax activation function to output probability-like predictions for each class. ```{r} # Define topology # ------------------------------------------------------------------------------ # Initialize a sequential model model3 <- keras_model_sequential() # Add layers to the model model3 %>% layer_conv_2d(filters = 32, kernel_size = c(5,5), activation = 'relu', input_shape = c(28, 28, 1)) %>% #layer_conv_2d(filters = 64, kernel_size = c(3,3), activation = 'relu') %>% layer_max_pooling_2d(pool_size = c(2, 2)) %>% layer_dropout(rate = 0.2) %>% layer_flatten() %>% layer_dense(units = 128, activation = 'relu') %>% #layer_dropout(rate = 0.5) %>% layer_dense(units = 10, activation = 'softmax') # Compile topology # ------------------------------------------------------------------------------ model3 %>% compile( loss = 'categorical_crossentropy', optimizer = 'adam', metrics = 'accuracy' ) # Fit topology # ------------------------------------------------------------------------------ # Fit the model history3 <- model3 %>% fit(train_x_sub_cnn, train_y_mtx_sub, epochs = 20, batch_size = 100, validation_split = 0.2) ``` ### Model performance ```{r} # Plot the model loss of the training data plot(history3) ``` The fact that the loss function of the validation dataset increases a little bit over the number of ephocs is a sign of overfitting behavious and we need to watch it carfully. ## Validation Let's now validate with our validation dataset ```{r} val_x_cnn<- val_x dim(val_x_cnn) <- c(nrow(val_x_cnn), 28, 28, 1) # Predict the classes for the test data val_y_cnn_pred <- model3 %>% predict_classes(val_x_cnn, batch_size = 128) # Plot contingegy table as heatmap plot table(val_y_cnn_pred, val_y) %>% reshape::melt() %>% ggplot(aes(x=val_y, y=val_y_cnn_pred)) + geom_tile(aes(fill = value)) + scale_x_continuous("Truth", breaks = 0:9) + scale_y_continuous("Prediction", breaks = 0:9) + geom_text(aes(label = value), color="white") + scale_size(range = c(15,50)) ``` Prediction accuracy: ```r prop.table(table(val_y_cnn_pred == val_y))[2]``` ## Prediction Let's finally run out CNN model on the TEST dataset provided by kaggle and generate the predictions for the competition. ```{r} test_cnn <- as.matrix(test) dim(test_cnn) <- c(nrow(test_cnn), 28, 28, 1) # Predict the classes for the test data classes <- model3 %>% predict_classes(test_cnn, batch_size = 128) results <- data.frame(ImageId=1:length(classes) , Label=classes) # output results write.table(results, file="prediction.txt" , quote =FALSE, col.names = TRUE, row.names = FALSE, sep=",") ``` The predicted label are exported in the .txt file ```prediction.txt```. Visually check the first 10 predictions and images ```{r} plot_mnist_images(test, 1:10, ncol=4, classes[1:10]) ``` ## Conlusions We have tested 3 different architectures and they all perform well without too much hyper-parameters optimization. MNIST is a simple dataset, so for the task a simple perceptron can do the job almost as well as CNN, although CNN are known from the literature to perform extremely well on more complex image classification problems. # Session info Information to allow reproducibility of the same workenviroment used in the analysis. ```{r} sessionInfo() ```