Practice 07 - MF FunkSVD with Cython.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recommender Systems 2021/22\n",
    "\n",
    "### Practice - FunkSVD and SVD++ implemented with Python and Cython\n",
    "\n",
    "FunkSVD is one of the simplest and most known matrix factorization models for rating prediction. It was proposed by Simon Funk in a now famous post on his website (with a cow grazing around Auckland): https://sifter.org/~simon/journal/20061211.html\n",
    "\n",
    "SVD++ is an extension of FunkSVD that learns also the biases: global, user and item."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Movielens10M: Verifying data consistency...\n",
      "Movielens10M: Verifying data consistency... Passed!\n",
      "DataReader: current dataset is: Movielens10M\n",
      "\tNumber of items: 10681\n",
      "\tNumber of users: 69878\n",
      "\tNumber of interactions in URM_all: 10000054\n",
      "\tValue range in URM_all: 0.50-5.00\n",
      "\tInteraction density: 1.34E-02\n",
      "\tInteractions per user:\n",
      "\t\t Min: 2.00E+01\n",
      "\t\t Avg: 1.43E+02\n",
      "\t\t Max: 7.36E+03\n",
      "\tInteractions per item:\n",
      "\t\t Min: 0.00E+00\n",
      "\t\t Avg: 9.36E+02\n",
      "\t\t Max: 3.49E+04\n",
      "\tGini Index: 0.57\n",
      "\n",
      "\tICM name: ICM_all, Value range: 1.00 / 69.00, Num features: 10126, feature occurrences: 128384, density 1.19E-03\n",
      "\tICM name: ICM_genres, Value range: 1.00 / 1.00, Num features: 20, feature occurrences: 21564, density 1.01E-01\n",
      "\tICM name: ICM_tags, Value range: 1.00 / 69.00, Num features: 10106, feature occurrences: 106820, density 9.90E-04\n",
      "\tICM name: ICM_year, Value range: 6.00E+00 / 2.01E+03, Num features: 1, feature occurrences: 10681, density 1.00E+00\n",
      "\n",
      "\n",
      "Warning: 87 (0.12 %) of 69878 users have no sampled items\n"
     ]
    }
   ],
   "source": [
    "from Data_manager.split_functions.split_train_validation_random_holdout import split_train_in_two_percentage_global_sample\n",
    "from Data_manager.Movielens.Movielens10MReader import Movielens10MReader\n",
    "\n",
    "data_reader = Movielens10MReader()\n",
    "data_loaded = data_reader.load_data()\n",
    "\n",
    "URM_all = data_loaded.get_URM_all()\n",
    "\n",
    "URM_train, URM_test = split_train_in_two_percentage_global_sample(URM_all, train_percentage = 0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<69878x10681 sparse matrix of type '<class 'numpy.float64'>'\n",
       "\twith 8000043 stored elements in Compressed Sparse Row format>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "URM_train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What do we need for FunkSVD and SVD++?\n",
    "\n",
    "* Loss function\n",
    "* User factor and Item factor matrices\n",
    "* Computing prediction\n",
    "* Update rule\n",
    "* Training loop and some patience\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_users, n_items = URM_train.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The two methods are based on two latent factor matrices $ W \\in R^{U \\times E} , V \\in R^{E \\times I}$ with E the embedding size, and biases $ \\mu \\in R , b_u \\in R^{U}, b_i \\in R^{I}$\n",
    "\n",
    "#### How to compute the predictions\n",
    "FunkSVD: $ \\hat{r}_{ui} = \\sum_{j=0}^{E} W_{uj}H_{ji}$\n",
    "\n",
    "SVD++: $ \\hat{r}_{ui} = \\mu + b_u + b_i + \\sum_{j=0}^{E} W_{uj}H_{ji}$\n",
    "\n",
    "\n",
    "#### The loss function we are interested in minimizing are\n",
    "\n",
    "$L_{FunkSVD} = |R - WH|_F + \\alpha|W|_2 + \\beta|H|_2$\n",
    "\n",
    "$L_{SVD++} = |R - (\\mu + b_u + b_i + WH)|_F + \\alpha|W|_2 + \\beta|H|_2 + \\gamma (\\mu + b_u + b_i) $\n",
    "\n",
    "Notice that in this case the loss is the Frobenius norm, not the 2 norm, hence we want to minimize the prediction error only for the existing ratings not the missing ones. In practice this means one samples only among the observed interactions.\n",
    "\n",
    "While this approach works well for rating prediction, it does not for ranking. A model must be trained also on negative data if it expected to learn how to distinguish between positive and negative data. A good strategy is to randomly sample unobserved interactions during training assigning them a rating of 0.\n",
    "\n",
    "#### Gradients\n",
    "\n",
    "$\\frac{\\partial}{\\partial W} L = -2(R - WH)H + 2\\alpha W $\n",
    "\n",
    "$\\frac{\\partial}{\\partial H} L = -2(R - WH)W + 2\\alpha H $\n",
    "\n",
    "$\\frac{\\partial}{\\partial \\mu} L = -2(R - WH) + 2\\gamma \\mu $\n",
    "\n",
    "$\\frac{\\partial}{\\partial b_u} L = -2(R - WH) + 2\\gamma b_u $\n",
    "\n",
    "$\\frac{\\partial}{\\partial b_i} L = -2(R - WH) + 2\\gamma b_i $\n",
    "\n",
    "#### The update is going to be (we can remove the coefficients)\n",
    "$ W = W - \\frac{\\partial}{\\partial W}$, or \n",
    "\n",
    "$ W = W + l((R - WH)H - \\alpha W)$, with $l$ the learning rate\n",
    "\n",
    "... and similarly for the other parameters: $H$, $\\mu$, $b_u$ and $b_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: We create the dense latent factor matrices\n",
    "### In a MF model you have two matrices, one with a row per user and the other with a column per item. The other dimension, columns for the first one and rows for the second one is called latent factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_factors = 10\n",
    "\n",
    "user_factors = np.random.random((n_users, num_factors))\n",
    "item_factors = np.random.random((n_items, num_factors))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.81782221, 0.32623483, 0.14184231, ..., 0.49757594, 0.36608581,\n",
       "        0.56447353],\n",
       "       [0.97263807, 0.12753114, 0.10881674, ..., 0.49452373, 0.31971186,\n",
       "        0.45143693],\n",
       "       [0.29415558, 0.18812572, 0.8754103 , ..., 0.09303427, 0.79470217,\n",
       "        0.25591769],\n",
       "       ...,\n",
       "       [0.19603707, 0.19376484, 0.8821932 , ..., 0.43127781, 0.67879198,\n",
       "        0.0986332 ],\n",
       "       [0.50865935, 0.64839804, 0.7321411 , ..., 0.63986981, 0.81789353,\n",
       "        0.24421367],\n",
       "       [0.18517082, 0.66859662, 0.98310463, ..., 0.77754424, 0.7527969 ,\n",
       "        0.69585431]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "user_factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.34346435, 0.2274097 , 0.41063912, ..., 0.59227481, 0.10454079,\n",
       "        0.00179158],\n",
       "       [0.71472291, 0.68239324, 0.25680354, ..., 0.79687973, 0.86179236,\n",
       "        0.76447388],\n",
       "       [0.98953008, 0.78386193, 0.15352571, ..., 0.12199865, 0.94654505,\n",
       "        0.86581823],\n",
       "       ...,\n",
       "       [0.2179097 , 0.71748106, 0.01822124, ..., 0.88652553, 0.75715021,\n",
       "        0.82133859],\n",
       "       [0.6848683 , 0.20372009, 0.03798786, ..., 0.02428757, 0.52423793,\n",
       "        0.25527629],\n",
       "       [0.686449  , 0.64856278, 0.05689972, ..., 0.22115793, 0.43541106,\n",
       "        0.56553645]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_factors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: We sample an interaction and compute the prediction of the current FunkSVD model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5006967"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "URM_train_coo = URM_train.tocoo()\n",
    "\n",
    "sample_index = np.random.randint(URM_train_coo.nnz)\n",
    "sample_index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(43649, 2479, 1.0)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "user_id = URM_train_coo.row[sample_index]\n",
    "item_id = URM_train_coo.col[sample_index]\n",
    "rating = URM_train_coo.data[sample_index]\n",
    "\n",
    "(user_id, item_id, rating)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.229260877924115"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predicted_rating = np.dot(user_factors[user_id,:], item_factors[item_id,:])\n",
    "predicted_rating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The first predicted rating is a random prediction, essentially\n",
    "\n",
    "### Step 3: We compute the prediction error and update the latent factor matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.229260877924115"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prediction_error = rating - predicted_rating\n",
    "prediction_error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The error is positive, so we need to increase the prediction our model computes. Meaning, we have to increase the values latent factor matrices\n",
    "\n",
    "### Which latent factors we modify? All the factors of the item and user we used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copy original value to avoid messing up the updates\n",
    "H_i = item_factors[item_id,:]\n",
    "W_u = user_factors[user_id,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.18565804, 0.93119643, 0.99308814, 0.86275186, 0.95212794,\n",
       "       0.1236942 , 0.52131264, 0.14175339, 0.40714719, 0.75519078])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.66744419, 0.44321036, 0.00922959, 0.51790834, 0.76529915,\n",
       "       0.6326098 , 0.75277024, 0.94348194, 0.77057488, 0.78101555])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W_u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Apply the update rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "learning_rate = 1e-4\n",
    "regularization = 1e-5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.41388688, -2.07588421, -2.21385262, -1.92330414, -2.12254922,\n",
       "       -0.27575297, -1.16214941, -0.31601473, -0.90764501, -1.68352507])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "user_factors_update = prediction_error * H_i - regularization * W_u\n",
    "user_factors_update"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.48790907, -0.98804082, -0.02058509, -1.15456144, -1.70606098,\n",
       "       -1.41025352, -1.67812645, -2.10326879, -1.71781651, -1.74109496])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_factors_update = prediction_error * W_u - regularization * H_i\n",
    "item_factors_update"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "user_factors[user_id,:] += learning_rate * user_factors_update \n",
    "item_factors[item_id,:] += learning_rate * item_factors_update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's check what the new prediction for the same user-item interaction would be"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.2272237878516385"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predicted_rating = np.dot(user_factors[user_id,:], item_factors[item_id,:])\n",
    "predicted_rating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The value is higher than before, we are moving in the right direction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### And now? Sample another interaction and repeat... a lot of times"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### WARNING: Initialization must be done with random non-zero values ... otherwise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Prediction is 0.00\n",
      "Prediction error is 1.00\n"
     ]
    }
   ],
   "source": [
    "user_factors = np.zeros((n_users, num_factors))\n",
    "item_factors = np.zeros((n_items, num_factors))\n",
    "\n",
    "predicted_rating = np.dot(user_factors[user_id,:], item_factors[item_id,:])\n",
    "\n",
    "print(\"Prediction is {:.2f}\".format(predicted_rating))\n",
    "\n",
    "prediction_error = rating - predicted_rating\n",
    "\n",
    "print(\"Prediction error is {:.2f}\".format(prediction_error))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_i = item_factors[item_id,:]\n",
    "W_u = user_factors[user_id,:]\n",
    "\n",
    "user_factors[user_id,:] += learning_rate * (prediction_error * H_i - regularization * W_u)\n",
    "item_factors[item_id,:] += learning_rate * (prediction_error * W_u - regularization * H_i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Prediction after the update is 0.00\n",
      "Prediction error is 1.00\n"
     ]
    }
   ],
   "source": [
    "predicted_rating = np.dot(user_factors[user_id,:], item_factors[item_id,:])\n",
    "\n",
    "print(\"Prediction after the update is {:.2f}\".format(predicted_rating))\n",
    "print(\"Prediction error is {:.2f}\".format(rating - predicted_rating))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Since the matrices are multiplied, if we initialize one of them as zero, the updates will always be zero and the model will not be able to learn."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's put all together in a training loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 100000 in 3.39 seconds, loss is 2.64. Samples per second 29534.75\n",
      "Iteration 200000 in 6.81 seconds, loss is 2.62. Samples per second 29347.44\n",
      "Iteration 300000 in 9.86 seconds, loss is 2.59. Samples per second 30417.08\n",
      "Iteration 400000 in 13.17 seconds, loss is 2.56. Samples per second 30380.99\n",
      "Iteration 500000 in 16.15 seconds, loss is 2.53. Samples per second 30967.11\n",
      "Iteration 600000 in 19.51 seconds, loss is 2.51. Samples per second 30750.35\n",
      "Iteration 700000 in 22.75 seconds, loss is 2.48. Samples per second 30767.58\n",
      "Iteration 800000 in 26.11 seconds, loss is 2.46. Samples per second 30644.79\n",
      "Iteration 900000 in 29.29 seconds, loss is 2.44. Samples per second 30724.14\n",
      "Iteration 1000000 in 32.55 seconds, loss is 2.41. Samples per second 30725.01\n"
     ]
    }
   ],
   "source": [
    "URM_train_coo = URM_train.tocoo()\n",
    "\n",
    "num_factors = 10\n",
    "learning_rate = 1e-4\n",
    "regularization = 1e-5\n",
    "\n",
    "user_factors = np.random.random((n_users, num_factors))\n",
    "item_factors = np.random.random((n_items, num_factors))\n",
    "\n",
    "loss = 0.0\n",
    "start_time = time.time()\n",
    "\n",
    "for sample_num in range(1000000):\n",
    "    \n",
    "    # Randomly pick sample\n",
    "    sample_index = np.random.randint(URM_train_coo.nnz)\n",
    "\n",
    "    user_id = URM_train_coo.row[sample_index]\n",
    "    item_id = URM_train_coo.col[sample_index]\n",
    "    rating = URM_train_coo.data[sample_index]\n",
    "\n",
    "    # Compute prediction\n",
    "    predicted_rating = np.dot(user_factors[user_id,:], item_factors[item_id,:])\n",
    "        \n",
    "    # Compute prediction error, or gradient\n",
    "    prediction_error = rating - predicted_rating\n",
    "    loss += prediction_error**2\n",
    "    \n",
    "    # Copy original value to avoid messing up the updates\n",
    "    H_i = item_factors[item_id,:]\n",
    "    W_u = user_factors[user_id,:]  \n",
    "    \n",
    "    user_factors_update = prediction_error * H_i - regularization * W_u\n",
    "    item_factors_update = prediction_error * W_u - regularization * H_i\n",
    "    \n",
    "    user_factors[user_id,:] += learning_rate * user_factors_update \n",
    "    item_factors[item_id,:] += learning_rate * item_factors_update    \n",
    "    \n",
    "    # Print some stats\n",
    "    if (sample_num +1)% 100000 == 0:\n",
    "        elapsed_time = time.time() - start_time\n",
    "        samples_per_second = sample_num/elapsed_time\n",
    "        print(\"Iteration {} in {:.2f} seconds, loss is {:.2f}. Samples per second {:.2f}\".format(sample_num+1, elapsed_time, loss/sample_num, samples_per_second))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What do we see? The loss generally goes down but may oscillate a bit.\n",
    "### How long do we train such a model?\n",
    "\n",
    "* An epoch: a complete loop over all the train data\n",
    "* Usually you train for multiple epochs. Depending on the algorithm and data 10s or 100s of epochs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated time with the previous training speed is 2603.74 seconds, or 43.40 minutes\n"
     ]
    }
   ],
   "source": [
    "estimated_seconds = 8e6 * 10 / samples_per_second\n",
    "print(\"Estimated time with the previous training speed is {:.2f} seconds, or {:.2f} minutes\".format(estimated_seconds, estimated_seconds/60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This model is relatively quick\n",
    "\n",
    "### Let's see what we can do with Cython\n",
    "### First step, just compile it. We do not have the data at compile time, so we put the loop in a function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "def do_some_training(URM_train):\n",
    "\n",
    "    URM_train_coo = URM_train.tocoo()\n",
    "    n_users, n_items = URM_train_coo.shape\n",
    "\n",
    "    num_factors = 10\n",
    "    learning_rate = 1e-4\n",
    "    regularization = 1e-5\n",
    "\n",
    "    user_factors = np.random.random((n_users, num_factors))\n",
    "    item_factors = np.random.random((n_items, num_factors))\n",
    "\n",
    "    loss = 0.0\n",
    "    start_time = time.time()\n",
    "\n",
    "    for sample_num in range(1000000):\n",
    "\n",
    "        # Randomly pick sample\n",
    "        sample_index = np.random.randint(URM_train_coo.nnz)\n",
    "\n",
    "        user_id = URM_train_coo.row[sample_index]\n",
    "        item_id = URM_train_coo.col[sample_index]\n",
    "        rating = URM_train_coo.data[sample_index]\n",
    "\n",
    "        # Compute prediction\n",
    "        predicted_rating = np.dot(user_factors[user_id,:], item_factors[item_id,:])\n",
    "\n",
    "        # Compute prediction error, or gradient\n",
    "        prediction_error = rating - predicted_rating\n",
    "        loss += prediction_error**2\n",
    "\n",
    "        # Copy original value to avoid messing up the updates\n",
    "        H_i = item_factors[item_id,:]\n",
    "        W_u = user_factors[user_id,:]  \n",
    "\n",
    "        user_factors_update = prediction_error * H_i - regularization * W_u\n",
    "        item_factors_update = prediction_error * W_u - regularization * H_i\n",
    "\n",
    "        user_factors[user_id,:] += learning_rate * user_factors_update \n",
    "        item_factors[item_id,:] += learning_rate * item_factors_update    \n",
    "\n",
    "        # Print some stats\n",
    "        if (sample_num +1)% 100000 == 0:\n",
    "            elapsed_time = time.time() - start_time\n",
    "            samples_per_second = sample_num/elapsed_time\n",
    "            print(\"Iteration {} in {:.2f} seconds, loss is {:.2f}. Samples per second {:.2f}\".format(sample_num+1, elapsed_time, loss/sample_num, samples_per_second))\n",
    "\n",
    "    return loss, samples_per_second"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 100000 in 2.80 seconds, loss is 2.62. Samples per second 35729.33\n",
      "Iteration 200000 in 5.65 seconds, loss is 2.58. Samples per second 35412.93\n",
      "Iteration 300000 in 8.46 seconds, loss is 2.55. Samples per second 35450.74\n",
      "Iteration 400000 in 11.42 seconds, loss is 2.53. Samples per second 35031.09\n",
      "Iteration 500000 in 14.21 seconds, loss is 2.50. Samples per second 35182.58\n",
      "Iteration 600000 in 16.98 seconds, loss is 2.48. Samples per second 35336.41\n",
      "Iteration 700000 in 19.87 seconds, loss is 2.45. Samples per second 35222.50\n",
      "Iteration 800000 in 22.69 seconds, loss is 2.43. Samples per second 35263.04\n",
      "Iteration 900000 in 25.48 seconds, loss is 2.41. Samples per second 35318.26\n",
      "Iteration 1000000 in 28.43 seconds, loss is 2.39. Samples per second 35168.39\n"
     ]
    }
   ],
   "source": [
    "loss, samples_per_second = do_some_training(URM_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated time with the previous training speed is 2274.77 seconds, or 37.91 minutes\n"
     ]
    }
   ],
   "source": [
    "estimated_seconds = 8e6 * 10 / samples_per_second\n",
    "print(\"Estimated time with the previous training speed is {:.2f} seconds, or {:.2f} minutes\".format(estimated_seconds, estimated_seconds/60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The compiler is just porting in C all operations that the python interpreter would have to perform, dynamic tiping included. Have a look at the html reports in the Cython_examples folder\n",
    "\n",
    "### Now try to add some types: If you use a variable only as a C object, use primitive tipes\n",
    "\n",
    "* cdef int namevar\n",
    "* cdef double namevar\n",
    "* cdef float namevar\n",
    "* cdef double[:] singledimensionarray\n",
    "* cdef double[:,:] bidimensionalmatrix\n",
    "\n",
    "### Some operations are still done with sparse matrices, those cannot be correctly optimized because the compiler does not know how what is the type of the data.\n",
    "\n",
    "### To address this, we create typed arrays in which we put the URM_train data\n",
    "For example, this operation: user_id = URM_train_coo.row[sample_index]\n",
    "\n",
    "Becomes:\n",
    "cdef int user_id\n",
    "cdef int[:] URM_train_coo_row = URM_train_coo.row\n",
    "user_id = URM_train_coo_row[sample_index]\n",
    "\n",
    "### We can also skip the creation of the items_in_user_profile array and replace the np.random call with the faster native C function rand()\n",
    "\n",
    "\n",
    "### We now use types for all main variables\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "from libc.stdlib cimport rand, srand, RAND_MAX\n",
    "\n",
    "def do_some_training(URM_train):\n",
    "\n",
    "    URM_train_coo = URM_train.tocoo()\n",
    "    n_users, n_items = URM_train_coo.shape\n",
    "    cdef int n_interactions = URM_train.nnz\n",
    "    \n",
    "    cdef int sample_num, sample_index, user_id, item_id, factor_index\n",
    "    cdef double rating, predicted_rating, prediction_error\n",
    "\n",
    "    cdef int num_factors = 10\n",
    "    cdef double learning_rate = 1e-4\n",
    "    cdef double regularization = 1e-5\n",
    "    \n",
    "    cdef int[:] URM_train_coo_row = URM_train_coo.row\n",
    "    cdef int[:] URM_train_coo_col = URM_train_coo.col\n",
    "    cdef double[:] URM_train_coo_data = URM_train_coo.data\n",
    "    \n",
    "    cdef double[:,:] user_factors = np.random.random((n_users, num_factors))\n",
    "    cdef double[:,:] item_factors = np.random.random((n_items, num_factors))\n",
    "    cdef double H_i, W_u\n",
    "    cdef double item_factors_update, user_factors_update\n",
    "                \n",
    "    cdef double loss = 0.0\n",
    "    cdef long start_time = time.time()\n",
    "\n",
    "    for sample_num in range(URM_train.nnz):\n",
    "\n",
    "        # Randomly pick sample\n",
    "        sample_index = rand() % n_interactions\n",
    "\n",
    "        user_id = URM_train_coo_row[sample_index]\n",
    "        item_id = URM_train_coo_col[sample_index]\n",
    "        rating = URM_train_coo_data[sample_index]\n",
    "\n",
    "        # Compute prediction\n",
    "        predicted_rating = 0.0\n",
    "        \n",
    "        for factor_index in range(num_factors):\n",
    "            predicted_rating += user_factors[user_id, factor_index] * item_factors[item_id, factor_index]\n",
    " \n",
    "        # Compute prediction error, or gradient\n",
    "        prediction_error = rating - predicted_rating\n",
    "        loss += prediction_error**2\n",
    "\n",
    "        # Copy original value to avoid messing up the updates\n",
    "        for factor_index in range(num_factors):\n",
    "            \n",
    "            H_i = item_factors[item_id,factor_index]\n",
    "            W_u = user_factors[user_id,factor_index]  \n",
    "\n",
    "            user_factors_update = prediction_error * H_i - regularization * W_u\n",
    "            item_factors_update = prediction_error * W_u - regularization * H_i\n",
    "\n",
    "            user_factors[user_id,factor_index] += learning_rate * user_factors_update \n",
    "            item_factors[item_id,factor_index] += learning_rate * item_factors_update    \n",
    "\n",
    "        # Print some stats\n",
    "        if (sample_num +1)% 500000 == 0:\n",
    "            elapsed_time = time.time() - start_time\n",
    "            samples_per_second = sample_num/elapsed_time\n",
    "            print(\"Iteration {} in {:.2f} seconds, loss is {:.2f}. Samples per second {:.2f}\".format(sample_num+1, elapsed_time, loss/sample_num, samples_per_second))\n",
    "\n",
    "    return loss, samples_per_second"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 500000 in 0.26 seconds, loss is 1.95. Samples per second 1920257.52\n",
      "Iteration 1000000 in 0.36 seconds, loss is 1.66. Samples per second 2806649.96\n",
      "Iteration 1500000 in 0.46 seconds, loss is 1.51. Samples per second 3252549.38\n",
      "Iteration 2000000 in 0.56 seconds, loss is 1.40. Samples per second 3583733.98\n",
      "Iteration 2500000 in 0.66 seconds, loss is 1.33. Samples per second 3782231.25\n",
      "Iteration 3000000 in 0.76 seconds, loss is 1.27. Samples per second 3963662.37\n",
      "Iteration 3500000 in 0.86 seconds, loss is 1.22. Samples per second 4085038.53\n",
      "Iteration 4000000 in 0.96 seconds, loss is 1.18. Samples per second 4185541.57\n",
      "Iteration 4500000 in 1.06 seconds, loss is 1.15. Samples per second 4263104.84\n",
      "Iteration 5000000 in 1.15 seconds, loss is 1.12. Samples per second 4330999.02\n",
      "Iteration 5500000 in 1.26 seconds, loss is 1.09. Samples per second 4358844.88\n",
      "Iteration 6000000 in 1.35 seconds, loss is 1.07. Samples per second 4435524.89\n",
      "Iteration 6500000 in 1.46 seconds, loss is 1.05. Samples per second 4462447.72\n",
      "Iteration 7000000 in 1.55 seconds, loss is 1.03. Samples per second 4505940.70\n",
      "Iteration 7500000 in 1.65 seconds, loss is 1.01. Samples per second 4538828.35\n",
      "Iteration 8000000 in 1.75 seconds, loss is 1.00. Samples per second 4573244.83\n"
     ]
    }
   ],
   "source": [
    "loss, samples_per_second = do_some_training(URM_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated time with the previous training speed is 17.49 seconds, or 0.29 minutes\n"
     ]
    }
   ],
   "source": [
    "estimated_seconds = 8e6 * 10 / samples_per_second\n",
    "print(\"Estimated time with the previous training speed is {:.2f} seconds, or {:.2f} minutes\".format(estimated_seconds, estimated_seconds/60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Works nicely, let's put an additional for loop to do multiple epochs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "from libc.stdlib cimport rand, srand, RAND_MAX\n",
    "\n",
    "def train_multiple_epochs(URM_train, learning_rate_input, regularization_input, n_epochs):\n",
    "    \n",
    "    URM_train_coo = URM_train.tocoo()\n",
    "    n_users, n_items = URM_train_coo.shape\n",
    "    cdef int n_interactions = URM_train.nnz\n",
    "    \n",
    "    cdef int sample_num, sample_index, user_id, item_id, factor_index\n",
    "    cdef double rating, predicted_rating, prediction_error\n",
    "\n",
    "    cdef int num_factors = 10\n",
    "    cdef double learning_rate = learning_rate_input\n",
    "    cdef double regularization = regularization_input\n",
    "    \n",
    "    cdef int[:] URM_train_coo_row = URM_train_coo.row\n",
    "    cdef int[:] URM_train_coo_col = URM_train_coo.col\n",
    "    cdef double[:] URM_train_coo_data = URM_train_coo.data\n",
    "\n",
    "    cdef double[:,:] user_factors = np.random.random((n_users, num_factors))\n",
    "    cdef double[:,:] item_factors = np.random.random((n_items, num_factors))\n",
    "    cdef double H_i, W_u\n",
    "    cdef double item_factors_update, user_factors_update\n",
    "                \n",
    "    cdef double loss = 0.0\n",
    "    cdef long start_time = time.time()\n",
    "    \n",
    "    for n_epoch in range(n_epochs):\n",
    "\n",
    "        loss = 0.0\n",
    "        start_time = time.time()\n",
    "\n",
    "        for sample_num in range(URM_train.nnz):\n",
    "\n",
    "            # Randomly pick sample\n",
    "            sample_index = rand() % n_interactions\n",
    "\n",
    "            user_id = URM_train_coo_row[sample_index]\n",
    "            item_id = URM_train_coo_col[sample_index]\n",
    "            rating = URM_train_coo_data[sample_index]\n",
    "\n",
    "            # Compute prediction\n",
    "            predicted_rating = 0.0\n",
    "\n",
    "            for factor_index in range(num_factors):\n",
    "                predicted_rating += user_factors[user_id, factor_index] * item_factors[item_id, factor_index]\n",
    "\n",
    "            # Compute prediction error, or gradient\n",
    "            prediction_error = rating - predicted_rating\n",
    "            loss += prediction_error**2\n",
    "\n",
    "            # Copy original value to avoid messing up the updates\n",
    "            for factor_index in range(num_factors):\n",
    "\n",
    "                H_i = item_factors[item_id,factor_index]\n",
    "                W_u = user_factors[user_id,factor_index]  \n",
    "\n",
    "                user_factors_update = prediction_error * H_i - regularization * W_u\n",
    "                item_factors_update = prediction_error * W_u - regularization * H_i\n",
    "\n",
    "                user_factors[user_id,factor_index] += learning_rate * user_factors_update \n",
    "                item_factors[item_id,factor_index] += learning_rate * item_factors_update    \n",
    "            \n",
    "        elapsed_time = time.time() - start_time\n",
    "        samples_per_second = sample_num/elapsed_time\n",
    "     \n",
    "        print(\"Epoch {} complete in in {:.2f} seconds, loss is {:.3E}. Samples per second {:.2f}\".format(n_epoch+1, time.time() - start_time, loss/sample_num, samples_per_second))\n",
    "\n",
    "    return np.array(user_factors), np.array(item_factors), loss, samples_per_second    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1 complete in in 1.78 seconds, loss is 5.799E-01. Samples per second 4493457.75\n",
      "Epoch 2 complete in in 2.50 seconds, loss is 3.020E-01. Samples per second 3200270.12\n",
      "Epoch 3 complete in in 2.09 seconds, loss is 2.202E-01. Samples per second 3829945.92\n",
      "Epoch 4 complete in in 1.74 seconds, loss is 1.884E-01. Samples per second 4590142.63\n",
      "Epoch 5 complete in in 2.48 seconds, loss is 1.716E-01. Samples per second 3220907.72\n",
      "Epoch 6 complete in in 2.20 seconds, loss is 1.611E-01. Samples per second 3636394.87\n",
      "Epoch 7 complete in in 1.84 seconds, loss is 1.538E-01. Samples per second 4337674.22\n",
      "Epoch 8 complete in in 2.40 seconds, loss is 1.483E-01. Samples per second 3329157.48\n",
      "Epoch 9 complete in in 2.04 seconds, loss is 1.442E-01. Samples per second 3919484.42\n",
      "Epoch 10 complete in in 1.78 seconds, loss is 1.407E-01. Samples per second 4484293.32\n"
     ]
    }
   ],
   "source": [
    "n_items = URM_train.shape[1]\n",
    "learning_rate = 1e-3\n",
    "regularization = 1e-5\n",
    "    \n",
    "user_factors, item_factors, loss, samples_per_second =  train_multiple_epochs(URM_train, learning_rate, regularization, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### From 20 minutes of training time to a few seconds..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}