Fastai Collaborative Filtering (from Scratch) with R and Reticulate

We had created a R notebook version of the first portion of movielens python notebook from the Fastai Deep Learning for Coders (Part 1) where high level fastai functions were used to build and fit the model. This notebook tries to create the R version of the second portion of the movielens python notebook where Jeremy creates the collaborative filtering model form scratch.

This content is covered in videos of lecture 5 and lecture 6. It will be helpful to listen to the lectures before going through this notebook since the concepts of the model and approach are discussed in the lecture and this notebook is just a replication attempt of the material from the course using R.

Initial Setup

# import R libraries
library(reticulate)
library(ggplot2)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

Notes about python setup and machine used are covered in this R notebook.

use_python("/home/paperspace/anaconda3/envs/fastai/bin/python", required = TRUE)
use_condaenv("fastai")
py_config()

python:         /home/paperspace/anaconda3/envs/fastai/bin/python
libpython:      /home/paperspace/anaconda3/envs/fastai/lib/libpython3.6m.so
pythonhome:     /home/paperspace/anaconda3/envs/fastai:/home/paperspace/anaconda3/envs/fastai
version:        3.6.4 |Anaconda, Inc.| (default, Dec 21 2017, 21:42:08)  [GCC 7.2.0]
numpy:          /home/paperspace/anaconda3/envs/fastai/lib/python3.6/site-packages/numpy
numpy_version:  1.13.3

NOTE: Python version was forced by use_python function

main = import_main()
bi = import_builtins()

# get relevant python imports
fstai_learner = import_from_path("fastai.learner", "../../fastai")
fstai_coldata = import_from_path("fastai.column_data", "../../fastai")

py_run_string("
from fastai.learner import *
from fastai.column_data import *
              ")

Get Data

The ratings dataset has the ratings for different users and movies. The movies dataset has the movie title information.

datapath = "../../data/ml-latest-small/"
ratings = read.csv(paste0(datapath, "ratings.csv"), stringsAsFactors = FALSE)
head(ratings)

movies = read.csv(paste0(datapath, "movies.csv"), stringsAsFactors = FALSE)
head(movies)

Working with PyTorch Tensors

Amat = matrix(c(1.0, 2.0, 3.0, 4.0), byrow = TRUE, ncol = 2)
Bmat = matrix(c(2.0, 2.0, 10.0, 10.0), byrow = TRUE, ncol = 2)
a = py$T(Amat)
a

 1  2
 3  4
[torch.cuda.FloatTensor of size 2x2 (GPU 0)]

b = py$T(Bmat)
b

  2   2
 10  10
[torch.cuda.FloatTensor of size 2x2 (GPU 0)]

In python when we use a*b, it works fine but won’t work when calling from R since R doesn’t know how to use * operator with Torch tensors. So here I have used the pyTorch mul function for multiplication

a$mul(b)

  2   4
 30  40
[torch.cuda.FloatTensor of size 2x2 (GPU 0)]

a$mul(b)$sum(1L) # note the use of 1L instead of 1 to ensure that we are passing integer

  6
 70
[torch.cuda.FloatTensor of size 2 (GPU 0)]

In PyTorch, we need to define a model as a Python class that inherits from nn.Module. It has a specific method forward which gives the recipe for computing the prediction given inputs. For example, if the prediction r is a dot product of input vectors u and m, then it would be defined in the following class.

py_run_string("
class DotProduct(nn.Module):
    def forward(self, u, m): return (u*m).sum(1)              
              ")

A model can then be defined as the instance of the class and calling the forward method with inputs with give the prediction.

model = main$DotProduct()
model$forward(a, b)

  6
 70
[torch.cuda.FloatTensor of size 2 (GPU 0)]

Data Preparation for Use in PyTorch Model

The ratings dataset has userId and movieId fields. But for passing to pyTorch, we create a sequential index for both.

u_unique = unique(ratings$userId)
user2idx = as.integer(seq(0, length(u_unique) - 1))
names(user2idx) = u_unique
user2idx[1:10]

 1  2  3  4  5  6  7  8  9 10 
 0  1  2  3  4  5  6  7  8  9 

m_unique = unique(ratings$movieId)
movie2idx = as.integer(seq(0, length(m_unique) - 1))
names(movie2idx) = m_unique
movie2idx[1:10]

  31 1029 1061 1129 1172 1263 1287 1293 1339 1343 
   0    1    2    3    4    5    6    7    8    9 

ratings$userIdx = user2idx[as.character(ratings$userId)]
ratings$movieIdx = movie2idx[as.character(ratings$movieId)]
n_users = length(u_unique)
n_movies = length(m_unique)
n_users; n_movies

[1] 671
[1] 9066

Model 1

Each user is \(i\) represent by an embedding vector \(u_i\) consisting of n_factor values. Similarly a movie \(j\) is represented by an embedding vector \(m_j\) consisting of n_factor values. The model of rating \(r_{ij}\) given by user \(i\) to movie \(j\) is: \[ r_{ij} = u_i^Tv_j\] The class below defines the model. The constructor for class passes other inputs and initializes model parameters.

py_run_string("
class EmbeddingDot(nn.Module):
    def __init__(self, n_users, n_movies, n_factors):
        super().__init__()
        self.u = nn.Embedding(n_users, n_factors)
        self.m = nn.Embedding(n_movies, n_factors)
        self.u.weight.data.uniform_(0,0.05)
        self.m.weight.data.uniform_(0,0.05)
        
    def forward(self, cats, conts):
        users,movies = cats[:,0],cats[:,1]
        u,m = self.u(users),self.m(movies)
        return (u*m).sum(1)              
              ")

The x dataframe just includes the sequential user id and movie id. y is a numpy array of ratings with type float32.

x = ratings[, c("userIdx", "movieIdx")]
y = np_array(ratings$rating)$astype(py$np$float32)

The list of validation data rows are selected and n_factors is set to 50.

val_idxs = py$get_cv_idxs(nrow(ratings))
val_idxs = as.integer(val_idxs)
n_factors = 50L

The data loader object is created

data = py$ColumnarModelData$from_data_frame(datapath, val_idxs, r_to_py(x), y, c("userIdx", "movieIdx"), 64L)

Model is defined using the defined class EmbeddingDot. Based on the model, the optimizer opt is defined.

wd=1e-5
model = py$EmbeddingDot(n_users, n_movies, n_factors)$cuda()
opt = py$optim$SGD(model$parameters(), 1e-1, weight_decay=wd, momentum=0.9)

The model details are listed below

model

EmbeddingDot(
  (u): Embedding(671, 50)
  (m): Embedding(9066, 50)
)

Model is fit. In Jupyter notebook, a widget shows a nice output of progress. That output doesn’t render properly within RStudio and the html generated document has too much output. For now, I have turned off output and am just storing the final validation loss metric (MSE loss in this case)

mdlfit = py$fit(model, data, 3L, opt, py$F$mse_loss)

paste0("MSE of validation set = ", round(mdlfit[[1]], 3))

[1] "MSE of validation set = 1.219"

Model run for some more epochs

py$set_lrs(opt, 0.01)
mdlfit = py$fit(model, data, 3L, opt, py$F$mse_loss)

paste0("MSE of validation set = ", round(mdlfit[[1]], 3))

[1] "MSE of validation set = 1.131"

Alternate way to check validation MSE is to explicitly get predictions from model for validation data and compare to validation data

yval_preds = py$predict(model, data$val_dl)
yval=data$val_y[,1]
mse_loss = mean((yval_preds - yval)**2)
mse_loss

[1] 1.131446

Model 2

This is model 1 with added bias term for users and movies. Each user is \(i\) represent by an embedding vector \(u_i\) consisting of n_factor values and a user bias value \(ub_i\). Similarly a movie \(j\) is represented by an embedding vector \(m_j\) consisting of n_factor values and movie bias value \(mb_j\). In addtion, the output is contrained to be between minimum rating r_min and maximum rating r_max using a sigmoid function. The model of rating \(r_{ij}\) given by user \(i\) to movie \(j\) is: \[ r_{ij} = \frac{e^{u_i^Tv_j + b_i + m_j}}{e^{u_i^Tv_j + b_i + m_j} + 1}(r_{max} - r_{min}) + r_{min} \]

py_run_string("
def get_emb(ni,nf):
    e = nn.Embedding(ni, nf)
    e.weight.data.uniform_(-0.01,0.01)
    return e
              
class EmbeddingDotBias(nn.Module):
    def __init__(self, n_users, n_movies, n_factors, min_rating, max_rating):
        super().__init__()
        (self.u, self.m, self.ub, self.mb) = [get_emb(*o) for o in [
              (n_users, n_factors), (n_movies, n_factors), (n_users,1), (n_movies,1)
              ]]
        self.max_rating = max_rating
        self.min_rating = min_rating
              
    def forward(self, cats, conts):
        users,movies = cats[:,0],cats[:,1]
        um = (self.u(users)* self.m(movies)).sum(1)
        res = um + self.ub(users).squeeze() + self.mb(movies).squeeze()
        res = F.sigmoid(res) * (self.max_rating-self.min_rating) + self.min_rating
        return res              
              ")

The model object opt is defined and fit

wd=2e-4
min_rating = min(ratings$rating)
max_rating = max(ratings$rating)
model = py$EmbeddingDotBias(n_users, n_movies, n_factors, min_rating, max_rating)$cuda()
opt = py$optim$SGD(model$parameters(), 1e-1, weight_decay=wd, momentum=0.9)

The model details are listed below

model

EmbeddingDotBias(
  (u): Embedding(671, 50)
  (m): Embedding(9066, 50)
  (ub): Embedding(671, 1)
  (mb): Embedding(9066, 1)
)

mdlfit = py$fit(model, data, 3L, opt, py$F$mse_loss)

paste0("MSE of validation set = ", round(mdlfit[[1]], 3))

[1] "MSE of validation set = 0.807"

Change the learning rate and refit

py$set_lrs(opt, 0.01)
mdlfit = py$fit(model, data, 3L, opt, py$F$mse_loss)

paste0("MSE of validation set = ", round(mdlfit[[1]], 3))

[1] "MSE of validation set = 0.801"

Model 3

This model is a deep learning model with the following layers:

1. Input - Concatenation of user embedding vector and movie embedding vector (size = 2*n_factors)
2. Dropout with dropout rate p1=0.05
3. Linear fully connected layer with output size nh = 10
4. Relu
5. Dropout with dropout rate p2=0.5
6. Linear fully connected layer with output size 1.
7. Apply sigmoid function and scale to be between min and max rating

py_run_string("
class EmbeddingNet(nn.Module):
    def __init__(self, n_users, n_movies, n_factors, min_rating, max_rating, nh=10, p1=0.05, p2=0.5):
        super().__init__()
        (self.u, self.m) = [get_emb(*o) for o in [
            (n_users, n_factors), (n_movies, n_factors)]]
        self.lin1 = nn.Linear(n_factors*2, nh)
        self.lin2 = nn.Linear(nh, 1)
        self.drop1 = nn.Dropout(p1)
        self.drop2 = nn.Dropout(p2)
        self.min_rating = min_rating
        self.max_rating = max_rating
        
    def forward(self, cats, conts):
        users,movies = cats[:,0],cats[:,1]
        x = self.drop1(torch.cat([self.u(users),self.m(movies)], dim=1))
        x = self.drop2(F.relu(self.lin1(x)))
        return F.sigmoid(self.lin2(x)) * (self.max_rating-self.min_rating+1) + self.min_rating-0.5
           ")

Model is defined and fit.

wd=1e-5
model = py$EmbeddingNet(n_users, n_movies, n_factors, min_rating, max_rating)$cuda()
opt = py$optim$Adam(model$parameters(), 1e-3, weight_decay=wd)

The model details are listed below

model

EmbeddingNet(
  (u): Embedding(671, 50)
  (m): Embedding(9066, 50)
  (lin1): Linear(in_features=100, out_features=10)
  (lin2): Linear(in_features=10, out_features=1)
  (drop1): Dropout(p=0.05)
  (drop2): Dropout(p=0.5)
)

mdlfit = py$fit(model, data, 3L, opt, py$F$mse_loss)

paste0("MSE of validation set = ", round(mdlfit[[1]], 3))

[1] "MSE of validation set = 0.784"

Learning rate is changed and model refit

py$set_lrs(opt, 0.01)
mdlfit = py$fit(model, data, 3L, opt, py$F$mse_loss)

paste0("MSE of validation set = ", round(mdlfit[[1]], 3))

[1] "MSE of validation set = 0.828"

Summary

This shows how a model can be developed from scratch. But the model would need to be defined as python class. I am guessing if somebody is developing a model from scratch, it might be better to do it just in python and create wrapper functions which can then by used by R.

sessionInfo()

R version 3.4.3 (2017-11-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.3 LTS

Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /home/paperspace/anaconda3/envs/fastai/lib/libmkl_intel_lp64.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] dplyr_0.7.4    ggplot2_2.2.1  reticulate_1.6

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.14     assertthat_0.2.0 R6_2.2.2         grid_3.4.3       plyr_1.8.4      
 [6] jsonlite_1.5     gtable_0.2.0     magrittr_1.5     scales_0.5.0     pillar_1.1.0    
[11] rlang_0.1.6      lazyeval_0.2.1   bindrcpp_0.2     tools_3.4.3      glue_1.2.0      
[16] munsell_0.4.3    yaml_2.1.18      compiler_3.4.3   pkgconfig_2.0.1  colorspace_1.3-2
[21] bindr_0.1.1      knitr_1.20       tibble_1.4.2