o_test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# O. SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This task is classical recommender system. It is well known since Netflix competition in 2006-2009. One\n",
    "of the competitors Simon Funk has [really nice description](http://sifter.org/simon/journal/20061211.html) of his method that uses SGD to find matrix factorization. It is good, because we don't need to deal with huge sparse matrices.\n",
    "\n",
    "Another useful thing was [suprise](http://surpriselib.com/) library. It does exactly what is required in our task and \n",
    "it has convenient methods to get sample data for testing purposes. We will try to implement our own algorithm\n",
    "to find matrix factorization and compare the results with those received using this library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from surprise import SVD, Dataset, accuracy\n",
    "from surprise.prediction_algorithms.baseline_only import BaselineOnly\n",
    "from surprise.model_selection import cross_validate, train_test_split"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can load smaller 100K dataset or larger 1M dataset (it has comparable number of ratings/users/movies as in our task).\n",
    "\n",
    "**100K** contains 100,000 ratings from 1000 users on 1700 movies.\n",
    "\n",
    "**1M** contains 1,000,000 ratings from 6000 users on 4000 movies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data = Dataset.load_builtin('ml-100k')\n",
    "data = Dataset.load_builtin('ml-1m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "trainset, testset = train_test_split(data, test_size=.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's try baseline algorithm, that doesn't take into account user-movie interaction but just uses global \n",
    "average rating, per movie average and per user average rating:\n",
    "$\\hat{r}_{ui} = \\mu + b_u + b_i$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimating biases using als...\n",
      "RMSE: 0.9069\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.90690923722444172"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = BaselineOnly()\n",
    "model.fit(trainset)\n",
    "predictions = model.test(testset)\n",
    "accuracy.rmse(predictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's try SVD approach that should produce better results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE: 0.8712\n",
      "CPU times: user 12.6 s, sys: 345 ms, total: 13 s\n",
      "Wall time: 13.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "model = SVD(n_factors=20, n_epochs=10, lr_all=0.01)\n",
    "model.fit(trainset)\n",
    "predictions = model.test(testset)\n",
    "accuracy.rmse(predictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "0.87 vs 0.90 which is a nice improvement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's implement our own SGD to find matrix factorization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_factors = 20\n",
    "n_epochs = 10\n",
    "init_mean = 0\n",
    "init_std = .1\n",
    "lr = .01\n",
    "reg = .02\n",
    "n_users = trainset.n_users\n",
    "n_items = trainset.n_items\n",
    "global_mean = trainset.global_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorithm is pretty straightforward. The idea is that we have 4 unknowns (2 matrices and 2 vectors):\n",
    "\n",
    "$b_u$ - per user average rating minus global average, vector of size `n_users`<br>\n",
    "$b_i$ - per movie(item) average rating minus global average, vector of size `n_items`<br>\n",
    "$P$ - matrix with rows representing users of size `n_users x n_factors`<br>\n",
    "$Q$ - matrix with rows representing movies of size `n_items x n_factors`<br>\n",
    "\n",
    "We initialize them with some random values and then iterate over each known user-movie-raiting tuple and compute \n",
    "error. Then we update just a little bit all the weights in matrices to minimize the error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 38s, sys: 1 s, total: 1min 39s\n",
      "Wall time: 1min 38s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "bu = np.zeros(n_users, np.double)\n",
    "bi = np.zeros(n_items, np.double)\n",
    "pu = np.random.normal(init_mean, init_std, (n_users, n_factors))\n",
    "qi = np.random.normal(init_mean, init_std, (n_items, n_factors))\n",
    "\n",
    "for current_epoch in range(n_epochs):\n",
    "    for u, i, r in trainset.all_ratings():\n",
    "        err = r - (global_mean + bu[u] + bi[i] + qi[i] @ pu[u])\n",
    "        \n",
    "        bu[u] += lr * (err - reg * bu[u])\n",
    "        bi[i] += lr * (err - reg * bi[i])\n",
    "\n",
    "        pu[u] += lr * (err * qi[i] - reg * pu[u])\n",
    "        qi[i] += lr * (err * pu[u] - reg * qi[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It takes a bit longer than surprise library. But let's hope the RMSE is in on par."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.87152340556689434"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rmse = 0\n",
    "predictions = []\n",
    "for u, i, r in testset:\n",
    "    est = global_mean\n",
    "    if u in trainset._raw2inner_id_users:\n",
    "        est += bu[trainset.to_inner_uid(u)]\n",
    "    if i in trainset._raw2inner_id_items:\n",
    "        est += bi[trainset.to_inner_iid(i)]\n",
    "    if u in trainset._raw2inner_id_users and i in trainset._raw2inner_id_items:\n",
    "        est += qi[trainset.to_inner_iid(i)] @ pu[trainset.to_inner_uid(u)]\n",
    "    lower_bound, higher_bound = trainset.rating_scale\n",
    "    est = min(higher_bound, est)\n",
    "    est = max(lower_bound, est)\n",
    "    predictions.append(dict(u=u, i=i, est=est))\n",
    "    rmse += (r - est) ** 2\n",
    "rmse = (rmse / len(testset)) ** 0.5\n",
    "rmse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've got pretty much the same RMSE 🙌 Small discrepancy might be due random initialization."
   ]
  }
 ],
 "metadata": {
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}