This module is designed to provide you an introduction to the keras API, deep learning and some of the key components that make DL algorithms run. Throughout this module you will learn about:

Perceptrons
Gradient descent
Activation functions
Learning rate & momentum
Model capacity (width vs. depth)
Learning curves
Batch size

Package Requirements

Loading

Let’s load the keras package along with a couple other packages we’ll use.

library(keras)       # for modeling
library(tidyverse)   # for wrangling & visualization
library(glue)        # for string literals

Installation

Normally, you will be starting out from scratch and need to install and set up keras on your own laptop. First, you need to install keras from CRAN. Once the package is installed, you need to install the Keras and TensorFlow Python packages, which is what the R Keras and TensorFlow packages communicate with. keras simplifies this with install_keras() which allows for:

both GPU & CPU options setups
installation in a virtual or conda environment
setup for Theano & CNTK backends rather than TensorFlow

See https://keras.rstudio.com/ for details.

# install keras if necessary
# install.packages("keras")

# default CPU-based installations of Keras and TensorFlow
# install_keras()

# for GPU installation
# install_keras(tensorflow = "gpu")

For this workshop we will be using a cloud environment to ensure we are all operating in common environment and Keras and TensorFlow have already been installed.

Simple Linear Regression

For our first example, let’s look at a very simple linear problem based on data with:

obs = 1,000
intercept = 30
slope = 5
with a little bit of random noise

n <- 1000   # n observations
b <- 30     # intercept
a <- 5      # slope

set.seed(123)
(df <- tibble(
  x = runif(n, min = -1, max = 1),
  y = b + a*x + rnorm(n)
))

Simple regression with OLS

I’m assuming we’re all familiar with the OLS methodology for linear regression modeling where our objective function (aka loss) is:

\[loss = MSE = \frac{1}{n}\sum^n_{i-1}(Y_i - \hat Y_i)^2\]

If we apply OLS to our data, we get:

estimated intercept = 30.01
estimated slope = 4.97
loss score (MSE) = 1.002

lm_model <- summary(lm(y ~ x, data = df))
lm_model

We can illustrate this model fit to our data:

mse <- lm_model[["sigma"]]

ggplot(df, aes(x, y)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE) +
  ggtitle(glue("MSE = {mse}"))

Simple regression with a perceptron

Now, let’s illustrate performing a similar process but with a basic building block of neural networks – the perceptron.

To model with keras we need our data to be in tensors. We’ll discuss tensors more later but for now just realize that:

1D tensor = vector
2D tensor = matrix

x <- as.matrix(df$x)
y <- df$y

Training a neural network model consists of 3 steps:

Define model architecture
Define how our model is going to learn
Train our model

1. Define model architecture

Defining an architecture includes defining the type of model and the arrangement of layers:

model <- keras_model_sequential() %>%
  layer_dense(units = 1, input_shape = ncol(x))

Define structure of model
- Sequential: linear layers (i.e. MLP, CNN, RNN, LSTM)
- Functional: adds more flexibility (i.e. combining CNN & LSTM to predict probability of cancer)
Define arrangement and shape of layers
- layer_dense: a single layer of nodes
- units = 1: a single perceptron unit in our layer
- input_shape: we need to tell our first layer how many inputs to expect

2. Define how our model is going to learn

The whole goal of training a neural network is to find the optimal set of parameter weights (aka coefficients in OLS-speak).

Our goal is to find weights that minimize the loss score

We define how our model is going to learn with compile():

model %>% compile(
  optimizer = "sgd",
  loss = "mse"
)

Optimizer: neural networks learn via backpropagation. There are various backpropagation algorithms but we’ll start with the most basic… stochastic gradient descent (SGD).
loss: how do we want to measure our model’s error. keras comes with many built-in loss functions and we can even create custom loss functions. Here, we’ll use MSE.

3. Train our model

The last thing we need is how to train our model, which we do with fit():

history <- model %>% fit(x, y, batch_size = 32, epochs = 10)

x: feature tensor
y: target tensor
batch_size: pick observations from our training data, perform forward pass, compute loss score, compute gradient, perform backward pass, update our weight (default = 32).
epoch: 1 epoch = one forward pass and one backward pass of all the training examples. We’re repeating that 10 times (default = 10).

Putting it all together

Let’s put all three steps together and train our model. Here’s a visual depiction:

And here’s the code:

# 1. Define model architecture
model <- keras_model_sequential() %>%
  layer_dense(units = 1, input_shape = ncol(x))

# 2. Define how our model is going to learn
model %>% compile(
  optimizer = "sgd",
  loss = "mse"
)

# 3. Train our model
history <- model %>% fit(x, y, epochs = 10)

Whereas in OLS we call the intercept and slope parameters “coefficients”, in neural networks we call them weights. We can see that after 10 epochs our weights are getting close to the underlying “truth” values (slope = 5, intercept = 30) but also notice that our MSE loss score is still decreasing after 10 epochs.

get_weights(model)

Models are object oriented

Note that modeling with keras/tensorflow in R may feel a bit different than other modeling packages you’ve used in R. Since keras/tensorflow are reticulated from Python, the model is an object oriented object with Python attributes.

Object oriented - our model object changes without assignment:
- model %>% compile() changed our model object by adding the optimizer and loss parameter arguments
- model %>% fit() will continue to build onto our existing model

# let's execute one more epoch. Note how our loss decreases from the last epoch
# above
model %>% fit(x, y, epoch = 1)

# our model weights get updated from this last epoch and continue to get closer
# to the underlying true values
get_weights(model)

Our model is a Python object - there will be things you can not directly access because they are Python objects. However, for most things that you want to access there will be a function to export them:

# the model weights are held in Python numpy arrays
model$weights

# we use helper functions to export these kinds of objects
get_weights(model)

Your Turn (3 min)

Fill in the blanks below and train the model for 25 epochs.
Explore the history object.
What are the final weights for this model? How do they compare to the underlying intercept (30) and slope (5)?

# 1. Define model architecture
model <- keras_model_sequential() %>%
  layer_dense(units = ___, input_shape = ____)

# 2. Define how our model is going to learn
model %>% compile(
  optimizer = "sgd",
  loss = ____
)

# 3. Train our model
history <- model %>% fit(x, y, epochs = ____)

Gradient descent

We can see the progression of the gradient descent process by examining the weights our model produces after each epoch. This is not something you will do often but, rather, helps make the gradient descent process more concrete.

# data frame to dump our results
model_est <- expand_grid(
  epoch = 1:25,
  a_sgd = NA,
  b_sgd = NA
)

# 1. Define model architecture
model <- keras_model_sequential() %>%
  layer_dense(units = 1, input_shape = ncol(x))

# 2. Define how our model is going to learn
model %>% compile(
  optimizer = "sgd",
  loss = "mse"
)

# 3. Train our model for 25 epochs and record the weights after each epoch
for (row in seq_len(nrow(model_est))) {
  
  history <- model %>% fit(x, y, epoch = 1, verbose = FALSE)
  
  current_wts <- get_weights(model)

  model_est[row, "a_sgd"] <- current_wts[[1]]
  model_est[row, "b_sgd"] <- current_wts[[2]]
}

The following table shows the estimate slope (a_sgd) and intercept (b_sgd) produced by SGD after each epoch.

model_est

history

We can visualize our linear model after each epoch with the following. Note how each epoch results in a linear prediction (dotted) that gets closer to the truth (blue line). After 25 epochs our model basically converges to the same results (yellow dotted line).

We can see this with our loss (MSE) that nearly equates the OLS MSE (1.00187).

epoch_pred <- merge(df, model_est, all = TRUE) %>%
  mutate(pred = b_sgd + a_sgd*x)

last_epoch <- filter(epoch_pred, epoch == max(epoch))

ggplot(data = df, aes(x, y)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = "lm", se = FALSE, size = 1.5) +
  geom_line(data = epoch_pred, aes(x, pred, group = epoch), lty = "dotted") +
  geom_line(data = last_epoch, aes(x, pred, group = epoch), 
            lty = "dotted", color = "yellow", size = 2) +
  ggtitle(glue("MSE = {history$metrics$loss}"))

Key Takeaways

A basic single perceptron computes the same transformation as OLS:

\[\hat y = bias + weight_1 \times x_1 + weight_2 \times x_2 + \dots + weight_n \times x_n\]

Neural networks learn via gradient descent - an iterative approach of predicting with a forward pass, measuring the gradient of the error, and performing a backward pass to update the weights based on the gradient.

Binary Classification

Let’s do the same process but now we’ll do so with a binary classification problem (i.e. predicting yes vs. no response).

set.seed(123)
generated <- mlbench::mlbench.simplex(n = 1000, d = 1, sd = .3)
x <- generated$x
y <- ifelse(generated$classes == 1, 0, 1)

(df <- tibble(x = as.vector(x), y = y))

Our generated data has some overlap so there is no linear seperation without having some error. Note than when discussing binary classification problems, we will mainly use the crossentropy (aka log loss) loss function.

glm_model <- glm(y ~ x, family = binomial(link = "logit"), data = df)
crossentropy <- MLmetrics::LogLoss(glm_model$fitted.values, df$y)

ggplot(df, aes(x, y)) +
  geom_point(aes(color = as.factor(y)), size = 2, show.legend = FALSE) +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +
  ggtitle(glue("crossentropy = {crossentropy}")) +
  ylab("probability y = 1")

Sigmoid Activation Function

When predicting a binary response, we typically want to predict a real value between 0-1 representing the probability of the positive binary class. Unfortunately our regular perceptron creates a linear transformation. However, we can apply an activation function to transform this linear transformation to a non-linear transformation.

When predicting a binary response, we use a sigmoid activation to convert our linear transformation to a 0-1 probability of the positive class.

\[sigmoid(y) = \frac{1}{1+e^{-y}}\]

When predicting a binary response, we need to make the following changes to our code:

Add activation = "sigmoid" to the layer that is predicting the output.
Note that since we are predicting the probability from 0-1 for our response, we keep units = 1.
loss - we change loss = "binary_crossentropy" to use the crossentropy / log loss objective function.

model <- keras_model_sequential() %>%
  layer_dense(units = 1, input_shape = ncol(x), activation = "sigmoid")

model %>% compile(
  optimizer = "sgd",
  loss = "binary_crossentropy"
)

(history <- model %>% fit(x, y, epochs = 50, verbose = FALSE))

We see that our loss is quite a bit off from our logistic regression model.

df %>%
  mutate(pred = predict(model, x) %>% as.vector()) %>%
  ggplot(aes(x, y)) +
  geom_point(aes(color = as.factor(y)), size = 2, show.legend = FALSE) +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +
  geom_line(aes(y = pred), lty = "dashed") +
  ggtitle(glue("crossentropy = {min(history$metrics$loss)}")) +
  ylab("probability of y = 1")

However, if we look at our loss scores, we see that they are still improving, its just taking a long time. Plus, it looks like there is more improvement that can be made to our loss.

plot(history)

Learning rate and momentum

An important parameter in gradient descent is the size of the steps which is controlled by the learning rate. If the learning rate is…

too small: the algorithm will take many iterations (steps) to find the minimum
too large: you might jump across the minimum and end up further away than when you started

The default learning rate for SGD is 0.01. Unfortunately with this rate, it will take over 1,000 epochs to reach a loss score comparable to logistic regression. However, we can customize our optimizer with optimizer_sdg():

model <- keras_model_sequential() %>%
  layer_dense(units = 1, input_shape = ncol(x), activation = "sigmoid")

model %>% compile(
  optimizer = optimizer_sgd(lr = 0.1),
  loss = "binary_crossentropy"
)

(history <- model %>% fit(x, y, epochs = 50, verbose = FALSE))

Another common approach to adjust our learning rate is to add momentum. Adding momentum allows our learning rate to adapt. Momentum simply adds a fraction of the previous weight update to the current one.

Let’s add some momentum to our learning rate. We see that our loss improves even more within the same number of epochs.

model <- keras_model_sequential() %>%
  layer_dense(units = 1, input_shape = ncol(x), activation = "sigmoid")

model %>% compile(
  optimizer = optimizer_sgd(lr = 0.1, momentum = 0.5),
  loss = "binary_crossentropy"
)

(history <- model %>% fit(x, y, epochs = 50, verbose = FALSE))

Your Turn (3 min)

Try different combinations of learning rate and momentum. A few rules of 👍:
- Typically, we start assessing learning rates in log values ranges of [1e-1, 1e-7] (i.e. 0.1, 0.01, …, 0.0000001).
- Momentum is typically > 0.5 and often in the 0.9-0.99 range.
Plot the loss learning curve.
How does your final loss compare to logistic regression?

model <- keras_model_sequential() %>%
  layer_dense(units = 1, input_shape = ncol(x), activation = ____)

model %>% compile(
  optimizer = ____,
  loss = "binary_crossentropy"
)

history <- model %>% fit(x, y, epochs = 50, verbose = FALSE)

Key takeaways

Activation functions:
- We use activation functions to transform the perceptron’s linear equation to a non-linear form.
- For binary classification problems we use the “sigmoid” activation to convert our predictions to a 0-1 probability.
Learning rate:
- We can control the rate of learning by increasing & decreasing the learning rate.
- We can make this learning rate adaptive to the curvature of our loss gradient by incorporating momentum.

Non-linear Patterns

As our datasets get larger or include non-linearities, our model needs to become more sophisticated. For this example, we’ll stick with one predictor variable but we’ll add a non-linearity component:

set.seed(123)
df <- tibble(
  x = seq(from = -1, to = 2 * pi, length = n),
  e = rnorm(n, sd = 0.2),
  y = sin(x) + e
)

ggplot(df, aes(x, y)) +
  geom_point(alpha = 0.5) +
  geom_smooth(se = FALSE)

Again, let’s extract our feature and target tensors:

x <- as.matrix(df$x)
y <- df$y

As our underlying model has more complexity, we add hidden layers to capture non-linearities and interactions. We call these neural network models multi-layer perceptrons (MLPs); also referred to as densely connected feed forward networks.

We can add a hidden layers by adding additional layer_dense() functions to our model architecture. For example, the following code would create an MLP with:

3 hidden layers:
- each hidden layer has 16 nodes
- only the first hidden layer requires input_shape
- each hidden layer uses a ReLU activation function (we’ll discuss shortly)
the last layer_dense() is always the output layer
- activation function for output layer is always dependent on the problem
  - regression: NULL
  - binary classification: activation = "signmoid"
  - multi-class classification: activation = "softmax"

model <- keras_model_sequential() %>%
  layer_dense(units = 16, input_shape = ncol(x), activation = "relu") %>%
  layer_dense(units = 16, activation = "relu") %>%
  layer_dense(units = 16, activation = "relu") %>%
  layer_dense(units = 1)

Why ReLU

The rectified linear activation function is very simple; if the linear transformation within the perceptron results in a negative number then the output is 0. If its positive then its that value.

\[ReLU = max(0, z)\]

Benefits (see http://proceedings.mlr.press/v15/glorot11a/glorot11a.pdf):

Simple geometric transformations can produce very complex patterns.
Computational simplicity (easy to compute the gradient)
Representational sparcity (forcing 0s results in sparse outputs)
Linearity (reduces vanishing gradient descent - discussed later)

Main Ingredients

Package Requirements

Loading

Installation

Simple Linear Regression

Simple regression with OLS

Simple regression with a perceptron

1. Define model architecture

2. Define how our model is going to learn

3. Train our model

Putting it all together

Models are object oriented

Your Turn (3 min)

Gradient descent

Key Takeaways

Binary Classification

Sigmoid Activation Function

Learning rate and momentum

Your Turn (3 min)

Key takeaways

Non-linear Patterns

Why ReLU