notebooks/02_model/surprise_svd_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": [
    "# Surprise Singular Value Decomposition (SVD)\n",
    "\n",
    "\n",
    "This notebook serves both as an introduction to the [Surprise](http://surpriselib.com/) library, and also introduces the 'SVD' algorithm which is very similar to ALS presented in the ALS deep dive notebook. This algorithm was heavily used during the Netflix Prize competition by the winning BellKor team."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Matrix factorization algorithm\n",
    "\n",
    "The SVD model algorithm is very similar to the ALS algorithm presented in the ALS deep dive notebook. The two differences between the two approaches are:\n",
    "\n",
    "- SVD additionally models the user and item biases (also called baselines in the litterature) from users and items.\n",
    "- The optimization technique in ALS is Alternating Least Squares (hence the name), while SVD uses stochastic gradient descent.\n",
    "\n",
    "### 1.1 The SVD model\n",
    "\n",
    "In ALS, the ratings are modeled as follows:\n",
    "\n",
    "$$\\hat r_{u,i} = q_{i}^{T}p_{u}$$\n",
    "\n",
    "SVD introduces two new scalar variables: the user biases $b_u$ and the item biases $b_i$. The user biases are supposed to capture the tendency of some users to rate items higher (or lower) than the average. The same goes for items: some items are usually rated higher than some others. The model is SVD is then as follows:\n",
    "\n",
    "$$\\hat r_{u,i} = \\mu + b_u + b_i + q_{i}^{T}p_{u}$$\n",
    "\n",
    "Where $\\mu$ is the global average of all the ratings in the dataset. The regularised optimization problem naturally becomes:\n",
    "\n",
    "$$ \\sum(r_{u,i} - (\\mu + b_u + b_i + q_{i}^{T}p_{u}))^2 +     \\lambda(b_i^2 + b_u^2 + ||q_i||^2 + ||p_u||^2)$$\n",
    "\n",
    "where $\\lambda$ is a the regularization parameter.\n",
    "\n",
    "\n",
    "### 1.2 Stochastic Gradient Descent\n",
    "\n",
    "Stochastic Gradient Descent (SGD) is a very common algorithm for optimization where the parameters (here the biases and the factor vectors) are iteratively incremented with the negative gradients w.r.t the optimization function. The algorithm essentially performs the following steps for a given number of iterations:\n",
    "\n",
    "\n",
    "$$b_u \\leftarrow b_u + \\gamma (e_{ui} - \\lambda b_u)$$\n",
    "$$b_i \\leftarrow b_i + \\gamma (e_{ui} - \\lambda b_i)$$  \n",
    "$$p_u \\leftarrow p_u + \\gamma (e_{ui} \\cdot q_i - \\lambda p_u)$$\n",
    "$$q_i \\leftarrow q_i + \\gamma (e_{ui} \\cdot p_u - \\lambda q_i)$$\n",
    "\n",
    "where $\\gamma$ is the learning rate and $e_{ui} =  r_{ui} - \\hat r_{u,i} = r_{u,i} - (\\mu + b_u + b_i + q_{i}^{T}p_{u})$ is the error made by the model for the pair $(u, i)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Surprise implementation of SVD\n",
    "\n",
    "SVD is implemented in the [Surprise](https://surprise.readthedocs.io/en/stable/) library as a recommender module. \n",
    "* Detailed documentations of the SVD module in Surprise can be found [here](https://surprise.readthedocs.io/en/stable/matrix_factorization.html#surprise.prediction_algorithms.matrix_factorization.SVD).\n",
    "* Source codes of the SVD implementation is available on the Surprise Github repository, which can be found [here](https://github.com/NicolasHug/Surprise/blob/master/surprise/prediction_algorithms/matrix_factorization.pyx)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Surprise SVD movie recommender\n",
    "\n",
    "We will use the MovieLens dataset, which is composed of integer ratings from 1 to 5. \n",
    "\n",
    "Surprise supports dataframes as long as they have three colums reprensenting the user ids, item ids, and the ratings (in this order)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "System version: 3.6.8 |Anaconda, Inc.| (default, Feb 11 2019, 15:03:47) [MSC v.1915 64 bit (AMD64)]\n",
      "Surprise version: 1.0.6\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "sys.path.append(\"../../\")\n",
    "import os\n",
    "import surprise\n",
    "import papermill as pm\n",
    "import pandas as pd\n",
    "\n",
    "from reco_utils.common.timer import Timer\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 (rmse, mae, rsquared, exp_var, map_at_k, ndcg_at_k, precision_at_k, \n",
    "                                                     recall_at_k, get_top_k_items)\n",
    "from reco_utils.recommender.surprise.surprise_utils import predict, compute_ranking_predictions\n",
    "\n",
    "print(\"System version: {}\".format(sys.version))\n",
    "print(\"Surprise version: {}\".format(surprise.__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'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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>rating</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>196</td>\n",
       "      <td>242</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>186</td>\n",
       "      <td>302</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>22</td>\n",
       "      <td>377</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>244</td>\n",
       "      <td>51</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\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": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Train the SVD Model\n",
    "\n",
    "Note that Surprise has a lot of built-in support for [cross-validation](https://surprise.readthedocs.io/en/stable/getting_started.html#use-cross-validation-iterators) or also [grid search](https://surprise.readthedocs.io/en/stable/getting_started.html#tune-algorithm-parameters-with-gridsearchcv) inspired scikit-learn, but we will here use the provided tools instead.\n",
    "\n",
    "We start by splitting our data into trainset and testset with the `python_random_split` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "train, test = python_random_split(data, 0.75)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Surprise needs to build an internal model of the data. We here use the `load_from_df` method to build a `Dataset` object, and then indicate that we want to train on all the samples of this dataset by using the `build_full_trainset` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<surprise.trainset.Trainset at 0x1a9cc607860>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 'reader' is being used to get rating scale (for MovieLens, the scale is [1, 5]).\n",
    "# 'rating_scale' parameter can be used instead for the later version of surprise lib:\n",
    "# https://github.com/NicolasHug/Surprise/blob/master/surprise/dataset.py\n",
    "train_set = surprise.Dataset.load_from_df(train, reader=surprise.Reader('ml-100k')).build_full_trainset()\n",
    "train_set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [SVD](https://surprise.readthedocs.io/en/stable/matrix_factorization.html#surprise.prediction_algorithms.matrix_factorization.SVD) has a lot of parameters. The most important ones are:\n",
    "- `n_factors`, which controls the dimension of the latent space (i.e. the size of the vectors $p_u$ and $q_i$). Usually, the quality of the training set predictions grows with as `n_factors` gets higher.\n",
    "- `n_epochs`, which defines the number of iteration of the SGD procedure.\n",
    "\n",
    "Note that both parameter also affect the training time.\n",
    "\n",
    "\n",
    "We will here set `n_factors` to `200` and `n_epochs` to `30`. To train the model, we simply need to call the `fit()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing epoch 0\n",
      "Processing epoch 1\n",
      "Processing epoch 2\n",
      "Processing epoch 3\n",
      "Processing epoch 4\n",
      "Processing epoch 5\n",
      "Processing epoch 6\n",
      "Processing epoch 7\n",
      "Processing epoch 8\n",
      "Processing epoch 9\n",
      "Processing epoch 10\n",
      "Processing epoch 11\n",
      "Processing epoch 12\n",
      "Processing epoch 13\n",
      "Processing epoch 14\n",
      "Processing epoch 15\n",
      "Processing epoch 16\n",
      "Processing epoch 17\n",
      "Processing epoch 18\n",
      "Processing epoch 19\n",
      "Processing epoch 20\n",
      "Processing epoch 21\n",
      "Processing epoch 22\n",
      "Processing epoch 23\n",
      "Processing epoch 24\n",
      "Processing epoch 25\n",
      "Processing epoch 26\n",
      "Processing epoch 27\n",
      "Processing epoch 28\n",
      "Processing epoch 29\n",
      "Took 19.879321813583374 seconds for training.\n"
     ]
    }
   ],
   "source": [
    "svd = surprise.SVD(random_state=0, n_factors=200, n_epochs=30, verbose=True)\n",
    "\n",
    "with Timer() as train_time:\n",
    "    svd.fit(train_set)\n",
    "\n",
    "print(\"Took {} seconds for training.\".format(train_time.interval))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Prediction\n",
    "\n",
    "Now that our model is fitted, we can call `predict` to get some predictions. `predict` returns an internal object `Prediction` which can be easily converted back to a dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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",
       "      <th>0</th>\n",
       "      <td>600.0</td>\n",
       "      <td>651.0</td>\n",
       "      <td>4.119</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>607.0</td>\n",
       "      <td>494.0</td>\n",
       "      <td>3.728</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>875.0</td>\n",
       "      <td>1103.0</td>\n",
       "      <td>4.225</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>648.0</td>\n",
       "      <td>238.0</td>\n",
       "      <td>4.225</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>113.0</td>\n",
       "      <td>273.0</td>\n",
       "      <td>4.043</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   userID  itemID  prediction\n",
       "0   600.0   651.0       4.119\n",
       "1   607.0   494.0       3.728\n",
       "2   875.0  1103.0       4.225\n",
       "3   648.0   238.0       4.225\n",
       "4   113.0   273.0       4.043"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predictions = predict(svd, test, usercol='userID', itemcol='itemID')\n",
    "predictions.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 Remove rated movies in the top k recommendations\n",
    "\n",
    "To compute ranking metrics, we need predictions on all user, item pairs. We remove though the items already watched by the user, since we choose not to recommend them again. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Took 28.51998782157898 seconds for prediction.\n"
     ]
    }
   ],
   "source": [
    "with Timer() as test_time:\n",
    "    all_predictions = compute_ranking_predictions(svd, train, usercol='userID', itemcol='itemID', remove_seen=True)\n",
    "    \n",
    "print(\"Took {} seconds for prediction.\".format(test_time.interval))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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",
       "      <th>75000</th>\n",
       "      <td>496</td>\n",
       "      <td>101</td>\n",
       "      <td>2.981</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75001</th>\n",
       "      <td>496</td>\n",
       "      <td>471</td>\n",
       "      <td>3.196</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75002</th>\n",
       "      <td>496</td>\n",
       "      <td>121</td>\n",
       "      <td>3.282</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75003</th>\n",
       "      <td>496</td>\n",
       "      <td>238</td>\n",
       "      <td>3.577</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75004</th>\n",
       "      <td>496</td>\n",
       "      <td>243</td>\n",
       "      <td>1.930</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       userID  itemID  prediction\n",
       "75000     496     101       2.981\n",
       "75001     496     471       3.196\n",
       "75002     496     121       3.282\n",
       "75003     496     238       3.577\n",
       "75004     496     243       1.930"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_predictions.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.5 Evaluate how well SVD performs \n",
    "\n",
    "The SVD algorithm was specifically designed to predict ratings as close as possible to their actual values. In particular, it is designed to have a very low RMSE (Root Mean Squared Error), computed as:\n",
    "\n",
    "$$\\sqrt{\\frac{1}{N} \\sum(\\hat{r_{ui}} - r_{ui})^2}$$\n",
    "\n",
    "As we can see, the RMSE and MAE (Mean Absolute Error) are pretty low (i.e. good), indicating that on average the error in the predicted ratings is less than 1. The RMSE is of course a bit higher, because high errors are penalized much more.\n",
    "\n",
    "For comparison with other models, we also display Top-k and ranking metrics (MAP, NDCG, etc.). Note however that the SVD algorithm was designed for achieving high accuracy, not for top-rank predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE:\t\t0.957953\n",
      "MAE:\t\t0.754764\n",
      "rsquared:\t0.286992\n",
      "exp var:\t0.287030\n",
      "----\n",
      "MAP:\t0.013018\n",
      "NDCG:\t0.099960\n",
      "Precision@K:\t0.095122\n",
      "Recall@K:\t0.032043\n"
     ]
    }
   ],
   "source": [
    "eval_rmse = rmse(test, predictions)\n",
    "eval_mae = mae(test, predictions)\n",
    "eval_rsquared = rsquared(test, predictions)\n",
    "eval_exp_var = exp_var(test, predictions)\n",
    "\n",
    "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",
    "\n",
    "print(\"RMSE:\\t\\t%f\" % eval_rmse,\n",
    "      \"MAE:\\t\\t%f\" % eval_mae,\n",
    "      \"rsquared:\\t%f\" % eval_rsquared,\n",
    "      \"exp var:\\t%f\" % eval_exp_var, sep='\\n')\n",
    "\n",
    "print('----')\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": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Record results with papermill for tests\n",
    "pm.record(\"rmse\", eval_rmse)\n",
    "pm.record(\"mae\", eval_mae)\n",
    "pm.record(\"rsquared\", eval_rsquared)\n",
    "pm.record(\"exp_var\", eval_exp_var)\n",
    "pm.record(\"map\", eval_map)\n",
    "pm.record(\"ndcg\", eval_ndcg)\n",
    "pm.record(\"precision\", eval_precision)\n",
    "pm.record(\"recall\", eval_recall)\n",
    "pm.record(\"train_time\", train_time.interval)\n",
    "pm.record(\"test_time\", test_time.interval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "1. Ruslan Salakhutdinov and Andriy Mnih. Probabilistic matrix factorization. 2008. URL: http://papers.nips.cc/paper/3208-probabilistic-matrix-factorization.pdf\n",
    "2. Yehuda Koren, Robert Bell, and Chris Volinsky. Matrix factorization techniques for recommender systems. 2009.\n",
    "3. Francesco Ricci, Lior Rokach, Bracha Shapira, and Paul B. Kantor. Recommender Systems Handbook. 1st edition, 2010."
   ]
  }
 ],
 "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}