notebooks/02_model/cornac_bpr_deep_dive.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<i>Copyright (c) Microsoft Corporation. All rights reserved.</i>\n",
    "\n",
    "<i>Licensed under the MIT License.</i>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Personalized Ranking (BPR)\n",
    "\n",
    "This notebook serves as an introduction to Bayesian Personalized Ranking (BPR) model for implicit feedback.  In this tutorial, we focus on learning the BPR model using matrix factorization approach, hence, the model is sometimes also named BPRMF.\n",
    "\n",
    "The implementation of the model is from [Cornac](https://github.com/PreferredAI/cornac), which is a framework for recommender systems with a focus on models leveraging auxiliary data (e.g., item descriptive text and image, social network, etc)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0 Global Settings and Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "System version: 3.6.8 |Anaconda, Inc.| (default, Feb 21 2019, 18:30:04) [MSC v.1916 64 bit (AMD64)]\n",
      "Cornac version: 1.1.2\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "sys.path.append(\"../../\")\n",
    "import os\n",
    "import cornac\n",
    "import papermill as pm\n",
    "import pandas as pd\n",
    "from reco_utils.dataset import movielens\n",
    "from reco_utils.dataset.python_splitters import python_random_split\n",
    "from reco_utils.evaluation.python_evaluation import map_at_k, ndcg_at_k, precision_at_k, recall_at_k\n",
    "from reco_utils.recommender.cornac.cornac_utils import predict_ranking\n",
    "from reco_utils.common.timer import Timer\n",
    "from reco_utils.common.constants import SEED\n",
    "\n",
    "print(\"System version: {}\".format(sys.version))\n",
    "print(\"Cornac version: {}\".format(cornac.__version__))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "# Select MovieLens data size: 100k, 1m, 10m, or 20m\n",
    "MOVIELENS_DATA_SIZE = '100k'\n",
    "\n",
    "# top k items to recommend\n",
    "TOP_K = 10\n",
    "\n",
    "# Model parameters\n",
    "NUM_FACTORS = 200\n",
    "NUM_EPOCHS = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 BPR Algorithm\n",
    "\n",
    "### 1.1 Personalized Ranking from Implicit Feedback\n",
    "\n",
    "The task of personalized ranking aims at providing each user a ranked list of items (recommendations).  This is very common in scenarios where recommender systems are based on implicit user behavior (e.g. purchases, clicks).  The available observations are only positive feedback where the non-observed ones are a mixture of real negative feedback and missing values.\n",
    "\n",
    "One usual approach for item recommendation is directly predicting a preference score $\\hat{x}_{u,i}$ given to item $i$ by user $u$.  BPR uses a different approach by using item pairs $(i, j)$ and optimizing for the correct ranking given preference of user $u$, thus, there are notions of *positive* and *negative* items.  The training data $D_S : U \\times I \\times I$ is defined as:\n",
    "\n",
    "$$D_S = \\{(u, i, j) \\mid i \\in I^{+}_{u} \\wedge j \\in I \\setminus I^{+}_{u}\\}$$\n",
    "\n",
    "where user $u$ is assumed to prefer $i$ over $j$ (i.e. $i$ is a *positive item* and $j$ is a *negative item*).\n",
    "\n",
    "\n",
    "### 1.2 Objective Function\n",
    "\n",
    "From the Bayesian perspective, BPR maximizes the posterior probability over the model parameters $\\Theta$ by optimizing the likelihood function $p(i >_{u} j | \\Theta)$ and the prior probability $p(\\Theta)$.\n",
    "\n",
    "$$p(\\Theta \\mid >_{u}) \\propto p(i >_{u} j \\mid \\Theta) \\times p(\\Theta)$$\n",
    "\n",
    "The joint probability of the likelihood over all users $u \\in U$ can be simplified to:\n",
    "\n",
    "$$ \\prod_{u \\in U} p(>_{u} \\mid \\Theta) = \\prod_{(u, i, j) \\in D_S} p(i >_{u} j \\mid \\Theta) $$\n",
    "\n",
    "The individual probability that a user $u$ prefers item $i$ to item $j$ can be defined as:\n",
    "\n",
    "$$ p(i >_{u} j \\mid \\Theta) = \\sigma (\\hat{x}_{uij}(\\Theta)) $$\n",
    "\n",
    "where $\\sigma$ is the logistic sigmoid:\n",
    "\n",
    "$$ \\sigma(x) = \\frac{1}{1 + e^{-x}} $$\n",
    "\n",
    "The preference scoring function $\\hat{x}_{uij}(\\Theta)$ could be an arbitrary real-valued function of the model parameter $\\Theta$.  Thus, it makes BPR a general framework for modeling the relationship between triplets $(u, i, j)$ where different model classes like matrix factorization could be used for estimating $\\hat{x}_{uij}(\\Theta)$.\n",
    "\n",
    "For the prior, one of the common pratices is to choose $p(\\Theta)$ following a normal distribution, which results in a nice form of L2 regularization in the final log-form of the objective function.\n",
    "\n",
    "$$ p(\\Theta) \\sim N(0, \\Sigma_{\\Theta}) $$\n",
    "\n",
    "To reduce the complexity of the model, all parameters $\\Theta$ are assumed to be independent and having the same variance, which gives a simpler form of the co-variance matrix $\\Sigma_{\\Theta} = \\lambda_{\\Theta}I$.  Thus, there are less number of hyperparameters to be determined.\n",
    "\n",
    "The final objective of the maximum posterior estimator:\n",
    "\n",
    "$$ J = \\sum_{(u, i, j) \\in D_S} \\text{ln } \\sigma(\\hat{x}_{uij}) - \\lambda_{\\Theta} ||\\Theta||^2 $$\n",
    "\n",
    "where $\\lambda_\\Theta$ are the model specific regularization paramerters.\n",
    "\n",
    "\n",
    "### 1.3 Learning with Matrix Factorization\n",
    "\n",
    "#### Stochastic Gradient Descent\n",
    "\n",
    "As the defined objective function is differentible, gradient descent based method for optimization is naturally adopted.  The gradient of the objective $J$ with respect to the model parameters:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{\\partial J}{\\partial \\Theta} & = \\sum_{(u, i, j) \\in D_S} \\frac{\\partial}{\\partial \\Theta} \\text{ln} \\ \\sigma(\\hat{x}_{uij}) - \\lambda_{\\Theta} \\frac{\\partial}{\\partial \\Theta} ||\\Theta||^2 \\\\\n",
    "& \\propto \\sum_{(u, i, j) \\in D_S} \\frac{-e^{-\\hat{x}_{uij}}}{1 + e^{-\\hat{x}_{uij}}} \\cdot  \\frac{\\partial}{\\partial \\Theta} \\hat{x}_{uij} - \\lambda_{\\Theta} \\Theta\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Due to slow convergence of full gradient descent, we prefer using stochastic gradient descent to optimize the BPR model.  For each triplet $(u, i, j) \\in D_S$, the update rule for the parameters:\n",
    "\n",
    "$$ \\Theta \\leftarrow \\Theta + \\alpha \\Big( \\frac{e^{-\\hat{x}_{uij}}}{1 + e^{-\\hat{x}_{uij}}} \\cdot \\frac{\\partial}{\\partial \\Theta} \\hat{x}_{uij} + \\lambda_\\Theta \\Theta \\Big) $$\n",
    "\n",
    "#### Matrix Factorization for Preference Approximation\n",
    "\n",
    "As mentioned earlier, the preference scoring function $\\hat{x}_{uij}(\\Theta)$ could be approximated by any real-valued function.  First, the estimator $\\hat{x}_{uij}$ is decomposed into:\n",
    "\n",
    "$$ \\hat{x}_{uij} = \\hat{x}_{ui} - \\hat{x}_{uj} $$\n",
    "\n",
    "The problem of estimating $\\hat{x}_{ui}$ is a standard collaborative filtering formulation, where matrix factorization approach has shown to be very effective.  The prediction formula can written as dot product between user feature vector $w_u$ and item feature vector $h_i$:\n",
    "\n",
    "$$ \\hat{x}_{ui} = \\langle w_u , h_i \\rangle = \\sum_{f=1}^{k} w_{uf} \\cdot h_{if} $$\n",
    "\n",
    "The  derivatives of matrix factorization with respect to the model parameters are:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial}{\\partial \\theta} \\hat{x}_{uij} = \n",
    "\\begin{cases}\n",
    "    (h_{if} - h_{jf})  & \\text{if } \\theta = w_{uf} \\\\\n",
    "    w_{uf}             & \\text{if } \\theta = h_{if} \\\\\n",
    "    -w_{uf}            & \\text{if } \\theta = h_{jf} \\\\\n",
    "    0                  & \\text{else}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "In theory, any kernel can be used to estimate $\\hat{x}_{ui}$ besides the dot product $ \\langle \\cdot , \\cdot \\rangle $.  For example, k-Nearest-Neighbor (kNN) has also been shown to achieve good performance.\n",
    "\n",
    "#### Analogies to AUC optimization\n",
    "\n",
    "By optimizing the objective function of BPR model, we effectively maximizing [AUC](https://towardsdatascience.com/understanding-auc-roc-curve-68b2303cc9c5) measurement.  To keep the notebook focused, please refer to the [paper](https://arxiv.org/ftp/arxiv/papers/1205/1205.2618.pdf) for details of the analysis (Section 4.1.1)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Cornac implementation of BPR\n",
    "\n",
    "BPR is implemented in the [Cornac](https://cornac.readthedocs.io/en/latest/index.html) framework as part of the model collections.\n",
    "* Detailed documentations of the BPR model in Cornac can be found [here](https://cornac.readthedocs.io/en/latest/models.html#bayesian-personalized-ranking-bpr).\n",
    "* Source codes of the BPR implementation is available on the Cornac Github repository, which can be found [here](https://github.com/PreferredAI/cornac/blob/master/cornac/models/bpr/recom_bpr.pyx).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Cornac BPR movie recommender\n",
    "\n",
    "\n",
    "### 3.1 Load and split data\n",
    "\n",
    "To evaluate the performance of item recommendation, we adopted the provided `python_random_split` tool for the consistency.  Data is randomly split into training and test sets with the ratio of 75/25.\n",
    "\n",
    "\n",
    "Note that Cornac also cover different [built-in schemes](https://cornac.readthedocs.io/en/latest/eval_methods.html) for model evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████████████████████████████████████████████████████████| 4.81k/4.81k [00:08<00:00, 590KB/s]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>userID</th>\n",
       "      <th>itemID</th>\n",
       "      <th>rating</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>196</td>\n",
       "      <td>242</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>186</td>\n",
       "      <td>302</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>22</td>\n",
       "      <td>377</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>244</td>\n",
       "      <td>51</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>166</td>\n",
       "      <td>346</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   userID  itemID  rating\n",
       "0     196     242     3.0\n",
       "1     186     302     3.0\n",
       "2      22     377     1.0\n",
       "3     244      51     2.0\n",
       "4     166     346     1.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = movielens.load_pandas_df(\n",
    "    size=MOVIELENS_DATA_SIZE,\n",
    "    header=[\"userID\", \"itemID\", \"rating\"]\n",
    ")\n",
    "\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "train, test = python_random_split(data, 0.75)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Cornac Dataset\n",
    "\n",
    "To work with models implemented in Cornac, we need to construct an object from [Dataset](https://cornac.readthedocs.io/en/latest/data.html#module-cornac.data.dataset) class.\n",
    "\n",
    "Dataset Class in Cornac serves as the main object that the models will interact with.  In addition to data transformations, Dataset provides a bunch of useful iterators for looping through the data, as well as supporting different negative sampling techniques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of users: 943\n",
      "Number of items: 1642\n"
     ]
    }
   ],
   "source": [
    "train_set = cornac.data.Dataset.from_uir(train.itertuples(index=False), seed=SEED)\n",
    "\n",
    "print('Number of users: {}'.format(train_set.num_users))\n",
    "print('Number of items: {}'.format(train_set.num_items))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Train the BPR model\n",
    "\n",
    "The BPR has a few important parameters that we need to consider:\n",
    "\n",
    "- `k`: controls the dimension of the latent space (i.e. the size of the vectors  $w_u$  and  $h_i$ ).\n",
    "- `max_iter`: defines the number of iterations of the SGD procedure.\n",
    "- `learning_rate`: controls the step size $\\alpha$ in the gradient update rules.\n",
    "- `lambda_reg`: controls the L2-Regularization $\\lambda$ in the objective function.\n",
    "\n",
    "Note that different values of `k` and `max_iter` will affect the training time.\n",
    "\n",
    "We will here set `k` to 200, `max_iter` to 100, `learning_rate` to 0.01, and `lambda_reg` to 0.001. To train the model, we simply need to call the `fit()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "bpr = cornac.models.BPR(\n",
    "    k=NUM_FACTORS,\n",
    "    max_iter=NUM_EPOCHS,\n",
    "    learning_rate=0.01,\n",
    "    lambda_reg=0.001,\n",
    "    verbose=True,\n",
    "    seed=SEED\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████| 100/100 [00:07<00:00, 13.27it/s, correct=92.19%, skipped=9.38%]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization finished!\n",
      "Took 7.6953 seconds for training.\n"
     ]
    }
   ],
   "source": [
    "with Timer() as t:\n",
    "    bpr.fit(train_set)\n",
    "print(\"Took {} seconds for training.\".format(t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 Prediction and Evaluation\n",
    "\n",
    "Now that our model is trained, we can produce the ranked lists for recommendation.  Every recommender models in Cornac provide `rate()` and `rank()` methods for predicting item rated value as well as item ranked list for a given user.  To make use of the current evaluation schemes, we will through `predict_rating()` and `predict_ranking()` functions inside `cornac_utils` to produce the predictions.\n",
    "\n",
    "Note that BPR model is effectively designed for item ranking.  Hence, we only measure the performance using ranking metrics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Took 1.7393803596496582 seconds for prediction.\n"
     ]
    }
   ],
   "source": [
    "with Timer() as t:\n",
    "    all_predictions = predict_ranking(bpr, train, usercol='userID', itemcol='itemID', remove_seen=True)\n",
    "print(\"Took {} seconds for prediction.\".format(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>userID</th>\n",
       "      <th>itemID</th>\n",
       "      <th>prediction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>75000</td>\n",
       "      <td>811</td>\n",
       "      <td>755</td>\n",
       "      <td>0.117239</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>75001</td>\n",
       "      <td>811</td>\n",
       "      <td>287</td>\n",
       "      <td>2.579992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>75002</td>\n",
       "      <td>811</td>\n",
       "      <td>181</td>\n",
       "      <td>3.743980</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>75003</td>\n",
       "      <td>811</td>\n",
       "      <td>96</td>\n",
       "      <td>1.959841</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>75004</td>\n",
       "      <td>811</td>\n",
       "      <td>83</td>\n",
       "      <td>1.122369</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       userID  itemID  prediction\n",
       "75000     811     755    0.117239\n",
       "75001     811     287    2.579992\n",
       "75002     811     181    3.743980\n",
       "75003     811      96    1.959841\n",
       "75004     811      83    1.122369"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_predictions.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MAP:\t0.109077\n",
      "NDCG:\t0.403395\n",
      "Precision@K:\t0.354989\n",
      "Recall@K:\t0.180183\n"
     ]
    }
   ],
   "source": [
    "k = 10\n",
    "eval_map = map_at_k(test, all_predictions, col_prediction='prediction', k=k)\n",
    "eval_ndcg = ndcg_at_k(test, all_predictions, col_prediction='prediction', k=k)\n",
    "eval_precision = precision_at_k(test, all_predictions, col_prediction='prediction', k=k)\n",
    "eval_recall = recall_at_k(test, all_predictions, col_prediction='prediction', k=k)\n",
    "\n",
    "print(\"MAP:\\t%f\" % eval_map,\n",
    "      \"NDCG:\\t%f\" % eval_ndcg,\n",
    "      \"Precision@K:\\t%f\" % eval_precision,\n",
    "      \"Recall@K:\\t%f\" % eval_recall, sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Record results with papermill for tests\n",
    "pm.record(\"map\", eval_map)\n",
    "pm.record(\"ndcg\", eval_ndcg)\n",
    "pm.record(\"precision\", eval_precision)\n",
    "pm.record(\"recall\", eval_recall)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "1. Rendle, S., Freudenthaler, C., Gantner, Z., & Schmidt-Thieme, L. (2009, June). BPR: Bayesian personalized ranking from implicit feedback. https://arxiv.org/ftp/arxiv/papers/1205/1205.2618.pdf\n",
    "2. Pan, R., Zhou, Y., Cao, B., Liu, N. N., Lukose, R., Scholz, M., & Yang, Q. (2008, December). One-class collaborative filtering. https://cseweb.ucsd.edu/classes/fa17/cse291-b/reading/04781145.pdf\n",
    "3. **Cornac** - A Comparative Framework for Multimodal Recommender Systems. https://cornac.preferred.ai/"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}