Practice 07b - MF AsySVD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recommender Systems 2022/23\n",
    "\n",
    "### Practice - AsySVD implemented with Python\n",
    "\n",
    "AsymmetricSVD is a model-based matrix factorization algorithm in which the user latent factors are represented as a function of their user profile and of a second item factor matrix."
   ]
  },
  {
   "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: 1.92E+03 / 2.01E+03, Num features: 1, feature occurrences: 10681, density 1.00E+00\n",
      "\n",
      "\n",
      "Warning: 69 (0.10 %) 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 AsySVD?\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, V \\in R^{I \\times E}$ with E the embedding size, and biases\n",
    "\n",
    "#### How to compute the predictions\n",
    "$ \\hat{r}_{ui} = \\sum_{k=0}^{E}\\sum_{j=0}^{I} r_{uj}W_{jk}H_{ki}$\n",
    "\n",
    "\n",
    "#### The loss function we are interested in minimizing is\n",
    "$L = ||R - RWH||_2 + \\alpha||W||_2 + \\beta||H||_2$\n",
    "\n",
    "#### Gradients\n",
    "\n",
    "$\\frac{\\partial}{\\partial W} L = -2(R - RWH)RH + 2\\alpha W $\n",
    "\n",
    "$\\frac{\\partial}{\\partial H} L = -2(R - RWH)RW + 2\\alpha H $\n",
    "\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 - RWH)RH - \\alpha W)$, with $l$ the learning rate\n"
   ]
  },
  {
   "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_profile_factors = np.random.random((n_items, 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.5157391 , 0.53619616, 0.40301611, ..., 0.12244558, 0.92860591,\n",
       "        0.79513609],\n",
       "       [0.07758103, 0.55094332, 0.10471291, ..., 0.61468311, 0.88315831,\n",
       "        0.890861  ],\n",
       "       [0.3762739 , 0.36680974, 0.9742875 , ..., 0.50090516, 0.37604545,\n",
       "        0.05379376],\n",
       "       ...,\n",
       "       [0.74352198, 0.71944669, 0.76706312, ..., 0.47981498, 0.54304556,\n",
       "        0.25849493],\n",
       "       [0.05371675, 0.67487388, 0.51322897, ..., 0.93147721, 0.50062325,\n",
       "        0.13569913],\n",
       "       [0.6590606 , 0.21083673, 0.20440548, ..., 0.1573705 , 0.20119   ,\n",
       "        0.06995089]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "USER_profile_factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.8499208 , 0.86941327, 0.23763236, ..., 0.71930711, 0.47732552,\n",
       "        0.01173423],\n",
       "       [0.85122788, 0.05495773, 0.51892409, ..., 0.89710259, 0.89205743,\n",
       "        0.02220628],\n",
       "       [0.64667339, 0.05925446, 0.07955873, ..., 0.42228871, 0.42887448,\n",
       "        0.06674025],\n",
       "       ...,\n",
       "       [0.41285108, 0.93938197, 0.84408214, ..., 0.9657658 , 0.10726293,\n",
       "        0.88987672],\n",
       "       [0.9882951 , 0.98591121, 0.21787309, ..., 0.16698451, 0.8014681 ,\n",
       "        0.21040748],\n",
       "       [0.17710744, 0.96545593, 0.93343823, ..., 0.63760433, 0.2530641 ,\n",
       "        0.52026922]])"
      ]
     },
     "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 model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7510999"
      ]
     },
     "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": [
       "(65811, 6959, 2.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_profile = URM_train[user_id]\n",
    "\n",
    "(user_id, item_id, rating)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.82478744, 1.75239499, 1.80015999, 1.79117253, 1.78692517,\n",
       "       1.73312823, 1.77264159, 1.71399974, 1.820476  , 1.82034022])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The estimated user factors may be divided by the square root of the profile length or the length itself\n",
    "# to improve learning stability (otherwise the dot product produces an embedding vector with very large numbers)\n",
    "USER_estimated_factors = user_profile.dot(USER_profile_factors).ravel()/user_profile.nnz\n",
    "USER_estimated_factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11.47797561031841"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predicted_rating = np.dot(USER_estimated_factors, 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": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-9.47797561031841"
      ]
     },
     "execution_count": 12,
     "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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copy original value to avoid messing up the updates\n",
    "H_all = ITEM_factors.copy()\n",
    "W_all = USER_profile_factors.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.8499208 , 0.86941327, 0.23763236, ..., 0.71930711, 0.47732552,\n",
       "        0.01173423],\n",
       "       [0.85122788, 0.05495773, 0.51892409, ..., 0.89710259, 0.89205743,\n",
       "        0.02220628],\n",
       "       [0.64667339, 0.05925446, 0.07955873, ..., 0.42228871, 0.42887448,\n",
       "        0.06674025],\n",
       "       ...,\n",
       "       [0.41285108, 0.93938197, 0.84408214, ..., 0.9657658 , 0.10726293,\n",
       "        0.88987672],\n",
       "       [0.9882951 , 0.98591121, 0.21787309, ..., 0.16698451, 0.8014681 ,\n",
       "        0.21040748],\n",
       "       [0.17710744, 0.96545593, 0.93343823, ..., 0.63760433, 0.2530641 ,\n",
       "        0.52026922]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.5157391 , 0.53619616, 0.40301611, ..., 0.12244558, 0.92860591,\n",
       "        0.79513609],\n",
       "       [0.07758103, 0.55094332, 0.10471291, ..., 0.61468311, 0.88315831,\n",
       "        0.890861  ],\n",
       "       [0.3762739 , 0.36680974, 0.9742875 , ..., 0.50090516, 0.37604545,\n",
       "        0.05379376],\n",
       "       ...,\n",
       "       [0.74352198, 0.71944669, 0.76706312, ..., 0.47981498, 0.54304556,\n",
       "        0.25849493],\n",
       "       [0.05371675, 0.67487388, 0.51322897, ..., 0.93147721, 0.50062325,\n",
       "        0.13569913],\n",
       "       [0.6590606 , 0.21083673, 0.20440548, ..., 0.1573705 , 0.20119   ,\n",
       "        0.06995089]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W_all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Apply the update rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "learning_rate = 1e-9    # Notice the low learning rate\n",
    "regularization = 1e-1   # Notice the high regularization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-8596.32279419, -8357.79801709, -8855.39794167, ...,\n",
       "        -8631.05266648, -8173.88499546, -8131.59700695],\n",
       "       [-8596.27897838, -8357.7994918 , -8855.36811135, ...,\n",
       "        -8631.10189024, -8173.8804507 , -8131.60657944],\n",
       "       [-8596.30884767, -8357.78107844, -8855.45506881, ...,\n",
       "        -8631.09051244, -8173.82973941, -8131.52287272],\n",
       "       ...,\n",
       "       [-8596.34557248, -8357.81634214, -8855.43434637, ...,\n",
       "        -8631.08840342, -8173.84643942, -8131.54334284],\n",
       "       [-8596.27659196, -8357.81188486, -8855.40896295, ...,\n",
       "        -8631.13356965, -8173.84219719, -8131.53106326],\n",
       "       [-8596.33712634, -8357.76548114, -8855.3780806 , ...,\n",
       "        -8631.05615897, -8173.81225387, -8131.52448843]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "user_factors_update = prediction_error * user_profile.dot(H_all) - regularization * W_all\n",
    "user_factors_update"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-8595.84455157, -8254.83796846, -8479.77437666, ...,\n",
       "        -8073.96004377, -8575.49801712, -8574.81185495],\n",
       "       [-8595.84468228, -8254.75652291, -8479.80250583, ...,\n",
       "        -8073.97782332, -8575.53949031, -8574.81290216],\n",
       "       [-8595.82422683, -8254.75695258, -8479.7585693 , ...,\n",
       "        -8073.93034193, -8575.49317201, -8574.81735556],\n",
       "       ...,\n",
       "       [-8595.8008446 , -8254.84496533, -8479.83502164, ...,\n",
       "        -8073.98468964, -8575.46101086, -8574.8996692 ],\n",
       "       [-8595.858389  , -8254.84961826, -8479.77240073, ...,\n",
       "        -8073.90481151, -8575.53043138, -8574.83172228],\n",
       "       [-8595.77727024, -8254.84757273, -8479.84395725, ...,\n",
       "        -8073.95187349, -8575.47559098, -8574.86270845]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_factors_update = prediction_error * user_profile.dot(W_all) - regularization * ITEM_factors\n",
    "item_factors_update"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "USER_profile_factors += learning_rate * user_factors_update \n",
    "ITEM_factors += 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": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "255.87654424685758"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "USER_estimated_factors = user_profile.dot(USER_profile_factors).ravel()/np.sqrt(user_profile.nnz)\n",
    "\n",
    "predicted_rating = np.dot(USER_estimated_factors, ITEM_factors[item_id,:])\n",
    "predicted_rating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 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": [
    "### 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 111.45 seconds, loss is 499.39. Samples per second 897.23\n",
      "Iteration 200000 in 226.92 seconds, loss is 262.08. Samples per second 881.37\n",
      "Iteration 300000 in 332.05 seconds, loss is 182.15. Samples per second 903.47\n",
      "Iteration 400000 in 437.68 seconds, loss is 142.12. Samples per second 913.92\n",
      "Iteration 500000 in 545.51 seconds, loss is 117.99. Samples per second 916.57\n",
      "Iteration 600000 in 651.23 seconds, loss is 101.91. Samples per second 921.33\n",
      "Iteration 700000 in 757.03 seconds, loss is 90.42. Samples per second 924.66\n",
      "Iteration 800000 in 863.30 seconds, loss is 81.79. Samples per second 926.67\n",
      "Iteration 900000 in 969.82 seconds, loss is 75.08. Samples per second 928.01\n",
      "Iteration 1000000 in 1075.50 seconds, loss is 69.67. Samples per second 929.79\n"
     ]
    }
   ],
   "source": [
    "URM_train_coo = URM_train.tocoo()\n",
    "\n",
    "num_factors = 10\n",
    "learning_rate = 1e-9    # Notice the low learning rate\n",
    "regularization = 1e-1   # Notice the high regularization\n",
    "\n",
    "USER_profile_factors = np.random.random((n_items, 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",
    "    user_profile = URM_train[user_id]\n",
    "    \n",
    "    if user_profile.nnz == 0:\n",
    "        continue \n",
    "        \n",
    "    USER_estimated_factors = user_profile.dot(USER_profile_factors).ravel()/np.sqrt(user_profile.nnz)\n",
    "    predicted_rating = np.dot(USER_estimated_factors, ITEM_factors[item_id,:])\n",
    "        \n",
    "    # Compute prediction error, or gradient\n",
    "    prediction_error = rating - predicted_rating\n",
    "    loss += prediction_error**2\n",
    "\n",
    "    if np.isnan(loss):\n",
    "        break \n",
    "        \n",
    "    # Copy original value to avoid messing up the updates\n",
    "    H_all = ITEM_factors.copy()\n",
    "    W_all = USER_profile_factors.copy()\n",
    "    \n",
    "    # Apply the updates\n",
    "    user_factors_update = prediction_error * user_profile.dot(H_all) - regularization * W_all\n",
    "    item_factors_update = prediction_error * user_profile.dot(W_all) - regularization * H_all\n",
    "    \n",
    "    USER_profile_factors += learning_rate * user_factors_update \n",
    "    ITEM_factors += 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",
    "### With higher learning rates or lower regularization you may see numerical instability (i.e., the loss suddendly explodes and then becomes nan, at which point some model parameters will also become none and the model is ruined)\n",
    "\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 86040.48 seconds (1434.01 minutes, 23.90 hours)\n"
     ]
    }
   ],
   "source": [
    "estimated_seconds = 8e6 * 10 / samples_per_second\n",
    "print(\"Estimated time with the previous training speed is {:.2f} seconds ({:.2f} minutes, {:.2f} hours)\".format(estimated_seconds, estimated_seconds/60, estimated_seconds/3600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### AsySVD can be very slow. Each sample requires to compute dozens (or hundreds) of dot products. Cython does not help in this case because most of the computational cost is already vectorized by numpy. Tools such as PyTorch may become useful in this case because they allow to better parallelize these operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}