Building A Neural Net from Scratch Using R – Part 1

Set Seed

Before we start, let’s set a seed value to ensure reproducibility of the results.

set.seed(69)

Architecture Definition

To understand the matrix multiplications better and keep the numbers digestible, we will describe a very simple 3-layer neural net i.e. a neural net with a single hidden layer. The \(1^{st}\) layer will take in the inputs and the \(3^{rd}\) layer will spit out an output.

The input layer will have two (input) neurons, the hidden layer four (hidden) neurons, and the output layer one (output) neuron.

Our input layer has two neurons because we’ll be passing two features (columns of a dataframe) as the input. A single output neuron because we’re performing binary classification. This means two output classes - 0 and 1. Our output will actually be a probability (a number that lies between 0 and 1). We’ll define a threshold for rounding off this probability to 0 or 1. For instance, this threshold can be 0.5.

In a deep neural net, multiple hidden layers are stacked together (hence the name “deep”). Each hidden layer can contain any number of neurons you want.

In this series, we’re implementing a single-layer neural net which, as the name suggests, contains a single hidden layer.

• n_x: the size of the input layer (set this to 2).
• n_h: the size of the hidden layer (set this to 4).
• n_y: the size of the output layer (set this to 1).

Neural networks flow from left to right, i.e. input to output. In the above example, we have two features (two columns from the input dataframe) that arrive at the input neurons from the first-row of the input dataframe. These two numbers are then multiplied by a set of weights (randomly initialized at first and later optimized).

An activation function is then applied on the result of this multiplication. This new set of numbers becomes the neurons in our hidden layer. These neurons are again multiplied by another set of weights (randomly initialized) with an activation function applied to this result. The final result we obtain is a single number. This is the prediction of our neural-net. It’s a number that lies between 0 and 1.

Once we have a prediction, we then compare it to the true output. To optimize the weights in order to make our predictions more accurate (because right now our input is being multiplied by random weights to give a random prediction), we need to first calculate how far off is our prediction from the actual value. Once we have this loss, we calculate the gradients with respect to each weight.

The gradients tell us the amount by which we need to increase or decrease each weight parameter in order to minimize the loss. All the weights in the network are updated as we repeat the entire process with the second input sample (second row).

After all the input samples have been used to optimize weights, we say that one epoch has passed. We repeat this process for multiple number of epochs till our loss stops decreasing.

At this point, you might be wondering what an activation function is. An activation function adds non-linearity to our network and enables it to learn complex features. If you look closely, a neural network consists of a bunch multiplications and additions. It’s linear and we know that a linear classification model will not be able to learn complex features in high dimensions.

Here are a few popular activation functions -

We will use tanh() and sigmoid() activation functions in our neural net. Because tanh() is already available in base-R, we will implement the sigmoid() function ourselves later on.

Dry Run

For now, let’s see how the numbers flow through the above described neural-net by writing out the equations for a single sample (one input row).

For one input sample \(x^{(i)}\) where \(i\) is the row-number:

First, we calculate the output \(Z\) from the input \(x\). We will tune the parameters \(W\) and \(b\). Here, the superscript in square brackets tell us the layer number and the one in parenthesis tell use the neuron number. For instance \(z^{[1] (i)}\) is the output from the \(i^{{th}}\) neuron of the \(1^{{st}}\) layer.

\[z^{[1] (i)} = W^{[1]} x^{(i)} + b^{[1] (i)}\tag{1}\]

Then we’ll pass this value through the tanh() activation function to get \(a\).

\[a^{[1] (i)} = \tanh(z^{[1] (i)})\tag{2}\]

After that, we’ll calculate the value for the final output layer using the hidden layer values.

\[z^{[2] (i)} = W^{[2]} a^{[1] (i)} + b^{[2] (i)}\tag{3}\]

Finally, we’ll pass this value through the sigmoid() activation function and obtain our output probability. \[\hat{y}^{(i)} = a^{[2] (i)} = \sigma(z^{ [2] (i)})\tag{4}\]

To obtain our prediction class from output probabilities, we round off the values as follows. \[y^{(i)}_{prediction} = \begin{cases} 1 & \mbox{if } a^{[2](i)} > 0.5 \\ 0 & \mbox{otherwise } \end{cases}\tag{5}\]

Once, we have the prediction probabilities, we’ll compute the loss in order to tune our parameters (\(w\) and \(b\) can be adjusted using gradient-descent).

Given the predictions on all the examples, we will compute the cost \(J\) the cross-entropy loss as follows:
\[J = - \frac{1}{m} \sum\limits_{i = 0}^{m} \large\left(\small y^{(i)}\log\left(\hat{y}^{(i)}\right) + (1-y^{(i)})\log\left(1- \hat{y}^{ (i)}\right) \large \right) \small \tag{6}\]

Once we have our loss, we need to calculate the gradients. I’ve calculated them for you so you don’t differentiate anything. We’ll directly use these values -

• \(dZ^{[2]} = A^{[2]} - Y\)
• \(dW^{[2]} = \frac{1}{m} dZ^{[2]}A^{[1]^T}\)
• \(db^{[2]} = \frac{1}{m}\sum dZ^{[2]}\)
• \(dZ^{[1]} = W^{[2]^T} * g^{[1]'} Z^{[1]}\) where \(g\) is the activation function.
• \(dW^{[1]} = \frac{1}{m}dZ^{[1]}X^{T}\)
• \(db^{[1]} = \frac{1}{m}\sum dZ^{[1]}\)

Now that we have the gradients, we will update the weights. We’ll multiply these gradients with a number known as the learning rate. The learning rate is represented by \(\alpha\).

• \(W^{[2]} = W^{[2]} - \alpha * dW^{[2]}\)
• \(b^{[2]} = b^{[2]} - \alpha * db^{[2]}\)
• \(W^{[1]} = W^{[1]} - \alpha * dW^{[1]}\)
• \(b^{[1]} = b^{[1]} - \alpha * db^{[1]}\)

This process is repeated multiple times until our model converges i.e. we have learned a good set of weights that fit our data well.

Load and Visualize the Data

Since, the goal of the series is to understand how neural-networks work behind the scene, we’ll use a small dataset so that our focus is on building our neural net.

We’ll use a planar dataset that looks like a flower. The output classes cannot be separated accurately using a straight line.

df <- read.csv(file = "planar_flower.csv")

df <- df[sample(nrow(df)), ]

head(df)

##           x1        x2 y
## 209  1.53856  3.242555 0
## 347 -0.05617 -0.808464 0
## 386 -3.85811  1.423514 1
## 112  0.82630  0.044276 1
## 104  0.31350  0.004274 1
## 111  2.28420  0.352476 1

train_test_split_index <- 0.8 * nrow(df)

train <- df[1:train_test_split_index,]
head(train)

##           x1        x2 y
## 209  1.53856  3.242555 0
## 347 -0.05617 -0.808464 0
## 386 -3.85811  1.423514 1
## 112  0.82630  0.044276 1
## 104  0.31350  0.004274 1
## 111  2.28420  0.352476 1

test <- df[(train_test_split_index+1): nrow(df),]
head(test)

##          x1       x2 y
## 210 -0.0352 -0.03489 0
## 348  2.7257 -0.54170 0
## 19  -2.2235  0.42137 1
## 362  2.3366 -0.40412 0
## 143 -1.4984  3.55267 0
## 4   -3.2264 -0.81648 0

Preprocess

Neural networks work best when the input values are standardized. So, we’ll scale all the values to to have their mean=0 and standard-deviation=1.

Standardizing input values speeds up the training and ensures faster convergence.

To standardize the input values, we’ll use the scale() function in R. Note that we’re standardizing the input values (X) only and not the output values (y).

X_train <- scale(train[, c(1:2)])

y_train <- train$y
dim(y_train) <- c(length(y_train), 1) # add extra dimension to vector

X_test <- scale(test[, c(1:2)])

y_test <- test$y
dim(y_test) <- c(length(y_test), 1) # add extra dimension to vector

The output below tells us the shape and size of our input data.

## Shape of X_train (row, column): 
##  320 2 
## Shape of y_train (row, column) : 
##  320 1 
## Number of training samples: 
##  320 
## Shape of X_test (row, column): 
##  80 2 
## Shape of y_test (row, column) : 
##  80 1 
## Number of testing samples: 
##  80

Because neural nets are made up of a bunch matrix multiplications, let’s convert our input and output to matrices from dataframes. While dataframes are a good way to represent data in a tabular form, we choose to convert to a matrix type because matrices are smaller than an equivalent dataframe and often speed up the computations.

We will also change the shape of X and y by taking its transpose. This will make the matrix calculations slightly more intuitive as we’ll see in the second part. There’s really no difference though. Some of you might find this way better, while others might prefer the non-transposed way. I feel this this makes more sense.

We’re going to use the as.matrix() method to construct out matrix. We’ll fill out matrix row-by-row.

X_train <- as.matrix(X_train, byrow=TRUE)
X_train <- t(X_train)
y_train <- as.matrix(y_train, byrow=TRUE)
y_train <- t(y_train)

X_test <- as.matrix(X_test, byrow=TRUE)
X_test <- t(X_test)
y_test <- as.matrix(y_test, byrow=TRUE)
y_test <- t(y_test)

Here are the shapes of our matrices after taking the transpose.

## Shape of X_train: 
##  2 320 
## Shape of y_train: 
##  1 320 
## Shape of X_test: 
##  2 80 
## Shape of y_test: 
##  1 80

Build a neural-net

Now that we’re done processing our data, let’s move on to building our neural net. As discussed above, we will broadly follow the steps outlined below.

1. Define the neural net architecture.
2. Initialize the model’s parameters from a random-uniform distribution.
3. Loop:
   • Implement forward propagation.
   • Compute loss.
   • Implement backward propagation to get the gradients.
   • Update parameters.

getLayerSize <- function(X, y, hidden_neurons, train=TRUE) {
  n_x <- dim(X)[1]
  n_h <- hidden_neurons
  n_y <- dim(y)[1]   
  
  size <- list("n_x" = n_x,
               "n_h" = n_h,
               "n_y" = n_y)
  
  return(size)
}

layer_size <- getLayerSize(X_train, y_train, hidden_neurons = 4)
layer_size

## $n_x
## [1] 2
## 
## $n_h
## [1] 4
## 
## $n_y
## [1] 1

initializeParameters <- function(X, list_layer_size){

    m <- dim(data.matrix(X))[2]
    
    n_x <- list_layer_size$n_x
    n_h <- list_layer_size$n_h
    n_y <- list_layer_size$n_y
        
    W1 <- matrix(runif(n_h * n_x), nrow = n_h, ncol = n_x, byrow = TRUE) * 0.01
    b1 <- matrix(rep(0, n_h), nrow = n_h)
    W2 <- matrix(runif(n_y * n_h), nrow = n_y, ncol = n_h, byrow = TRUE) * 0.01
    b2 <- matrix(rep(0, n_y), nrow = n_y)
    
    params <- list("W1" = W1,
                   "b1" = b1, 
                   "W2" = W2,
                   "b2" = b2)
    
    return (params)
}

init_params <- initializeParameters(X_train, layer_size)
lapply(init_params, function(x) dim(x))

## $W1
## [1] 4 2
## 
## $b1
## [1] 4 1
## 
## $W2
## [1] 1 4
## 
## $b2
## [1] 1 1

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

forwardPropagation <- function(X, params, list_layer_size){
    
    m <- dim(X)[2]
    n_h <- list_layer_size$n_h
    n_y <- list_layer_size$n_y
    
    W1 <- params$W1
    b1 <- params$b1
    W2 <- params$W2
    b2 <- params$b2
    
    b1_new <- matrix(rep(b1, m), nrow = n_h)
    b2_new <- matrix(rep(b2, m), nrow = n_y)
    
    Z1 <- W1 %*% X + b1_new
    A1 <- sigmoid(Z1)
    Z2 <- W2 %*% A1 + b2_new
    A2 <- sigmoid(Z2)
    
    cache <- list("Z1" = Z1,
                  "A1" = A1, 
                  "Z2" = Z2,
                  "A2" = A2)

    return (cache)
}

fwd_prop <- forwardPropagation(X_train, init_params, layer_size)
lapply(fwd_prop, function(x) dim(x))

## $Z1
## [1]   4 320
## 
## $A1
## [1]   4 320
## 
## $Z2
## [1]   1 320
## 
## $A2
## [1]   1 320

End of Part 1

We have reached the end of Part1 here. In the next and final part, we will implement backpropagation and evaluate our model. Stay tuned!