task_02_rating_prediction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project task 02: Restaurant recommendation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy.sparse as sp\n",
    "import numpy as np\n",
    "import scipy\n",
    "from scipy.sparse.linalg import svds\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from sklearn.linear_model import LinearRegression, Ridge\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this task is to recommend restaurants to users based on the rating data in the Yelp dataset. For this, we try to predict the rating a user will give to a restaurant they have not been to yet based on a latent factor model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Load data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First download `ratings.npy` from Piazza ([download link](https://syncandshare.lrz.de/dl/fiKMoxRNusLoFpFHkXXEgvdZ/ratings.npy))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ratings = np.load(\"ratings.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[101968,   1880,      1],\n",
       "       [101968,    284,      5],\n",
       "       [101968,   1378,      2],\n",
       "       ..., \n",
       "       [ 72452,   2100,      4],\n",
       "       [ 72452,   2050,      5],\n",
       "       [ 74861,   3979,      5]], dtype=uint32)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We have triplets of (user, restaurant, rating).\n",
    "ratings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we transform the data into a matrix of dimension [N, D], where N is the number of users and D is the number of restaurants in the dataset.  \n",
    "We **strongly recommend** to load the data as a sparse matrix to avoid out-of-memory issues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<337867x5899 sparse matrix of type '<class 'numpy.uint32'>'\n",
       "\twith 929606 stored elements in Compressed Sparse Row format>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Store the matrix into the variable M\n",
    "\n",
    "### YOUR CODE HERE ###\n",
    "ratings_matrix = sp.csr_matrix((ratings[:,2], (ratings[:,0], ratings[:,1])))\n",
    "ratings_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Preprocess the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the preprocessing step, we recursively remove all users and restaurants with 10 or less ratings. \n",
    "\n",
    "Then, we randomly select 200 data points for the validation and test sets, respectively.\n",
    "\n",
    "After this, we subtract the mean rating for each users to account for this global effect.   \n",
    "**Hint**: Some entries might become zero in this process -- but these entries are different than the 'unknown' zeros in the matrix. Store the indices of which we have data available in a separate variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#function to remove rows from a sparse matrix given indices\n",
    "def delete_rows_csr1(mat, indices):\n",
    "    \"\"\"\n",
    "    Remove the rows denoted by ``indices`` form the CSR sparse matrix ``mat``.\n",
    "    \"\"\"\n",
    "    if not isinstance(mat, sp.csr_matrix):\n",
    "        raise ValueError(\"works only for CSR format -- use .tocsr() first\")\n",
    "    indices = list(indices)\n",
    "    mask = np.ones(mat.shape[0], dtype=bool)\n",
    "    mask[indices] = False\n",
    "    return mat[mask]\n",
    "\n",
    "def delete_cols_csr1(mat, indices):\n",
    "    \"\"\"\n",
    "    Remove the cols denoted by ``indices`` form the CSR sparse matrix ``mat``.\n",
    "    \"\"\"\n",
    "    if not isinstance(mat, sp.csr_matrix):\n",
    "        raise ValueError(\"works only for CSR format -- use .tocsr() first\")\n",
    "    indices = list(indices)\n",
    "    mask = np.ones(mat.shape[1], dtype=bool)\n",
    "    mask[indices] = False\n",
    "    return mat[:,mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cold_start_preprocessing(matrix, min_entries):\n",
    "    \"\"\"\n",
    "    Recursively removes rows and columns from the input matrix which have less than min_entries nonzero entries.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    matrix      : sp.spmatrix, shape [N, D]\n",
    "                  The input matrix to be preprocessed.\n",
    "    min_entries : int\n",
    "                  Minimum number of nonzero elements per row and column.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    matrix      : sp.spmatrix, shape [N', D']\n",
    "                  The pre-processed matrix, where N' <= N and D' <= D\n",
    "        \n",
    "    \"\"\"\n",
    "    print(\"Shape before: {}\".format(matrix.shape))\n",
    "    \n",
    "    #matrix = np.asmatrix(matrix.toarray())\n",
    "#     while(~((((matrix>0).sum(0).A1 > min_entries).all())&(((matrix>0).sum(1).A1 > min_entries).all()))):\n",
    "#     rowise = np.delete(matrix, np.where(((matrix>0).sum(1))<=min_entries), 0)\n",
    "#     colwise = np.delete(rowise, np.where(((rowise>0).sum(0))<=min_entries), 1)\n",
    "#     matrix = colwise\n",
    "    #If above while loop doesnt work please run the below commented while loop\n",
    "    while(~((((matrix>0).sum(0) > min_entries).all())&(((matrix>0).sum(1) > min_entries).all()))):\n",
    "        rr =  delete_rows_csr1(matrix,np.where((matrix>0).sum(1) <=min_entries)[0] )\n",
    "        cc1 = delete_cols_csr1(rr,np.where((rr>0).sum(0) <=min_entries)[1])\n",
    "        matrix = cc1\n",
    "        \n",
    "    #matrix = sp.csr_matrix(matrix)\n",
    "    new_shape = matrix.shape\n",
    "    \n",
    "    print(\"Shape after: {}\".format(matrix.shape))\n",
    "    \n",
    "    nnz = matrix>0\n",
    "    \n",
    "    assert (nnz.sum(0).A1 > min_entries).all()\n",
    "    assert (nnz.sum(1).A1 > min_entries).all()\n",
    "    return matrix, new_shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def shift_user_mean(matrix):\n",
    "    \"\"\"\n",
    "    Subtract the mean rating per user from the non-zero elements in the input matrix.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    matrix : sp.spmatrix, shape [N, D]\n",
    "             Input sparse matrix.\n",
    "    Returns\n",
    "    -------\n",
    "    matrix : sp.spmatrix, shape [N, D]\n",
    "             The modified input matrix.\n",
    "    \n",
    "    user_means : np.array, shape [N, 1]\n",
    "                 The mean rating per user that can be used to recover the absolute ratings from the mean-shifted ones.\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    ### YOUR CODE HERE ###\n",
    "    reduced_matrix = matrix\n",
    "    #if np.diff(reduced_matrix.indptr)>0:\n",
    "    user_means = np.array(reduced_matrix.sum(axis=1).squeeze())[0]/np.diff(reduced_matrix.indptr)\n",
    "    user_means[ ~ np.isfinite( user_means )] = 0 \n",
    "    #print(type(user_means))\n",
    "    #print(user_means)\n",
    "    #else:\n",
    "    #user_means = 0\n",
    "    #matrix_mean_diag = sp.diags(user_means,0)\n",
    "    matrix_mean_diag = np.diag(user_means)\n",
    "    #print(matrix_mean_diag.shape)\n",
    "    reduced_matrix_copy = reduced_matrix.copy()\n",
    "    reduced_matrix_copy.data = np.ones_like(reduced_matrix_copy.data)\n",
    "    matrix_mean_diag = sp.csr_matrix(matrix_mean_diag)\n",
    "    M_demean = (reduced_matrix-(matrix_mean_diag*reduced_matrix_copy).todense())\n",
    "    assert np.all(np.isclose(M_demean.mean(1), 0))\n",
    "    return M_demean, user_means\n",
    "    #assert np.all(np.isclose(matrix.mean(1), 0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def split_data(matrix, n_validation, n_test):\n",
    "    \"\"\"\n",
    "    Extract validation and test entries from the input matrix. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    matrix          : sp.spmatrix, shape [N, D]\n",
    "                      The input data matrix.\n",
    "    n_validation    : int\n",
    "                      The number of validation entries to extract.\n",
    "    n_test          : int\n",
    "                      The number of test entries to extract.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    matrix_split    : sp.spmatrix, shape [N, D]\n",
    "                      A copy of the input matrix in which the validation and test entries have been set to zero.\n",
    "    \n",
    "    val_idx         : tuple, shape [2, n_validation]\n",
    "                      The indices of the validation entries.\n",
    "    \n",
    "    test_idx        : tuple, shape [2, n_test]\n",
    "                      The indices of the test entries.\n",
    "    \n",
    "    val_values      : np.array, shape [n_validation, ]\n",
    "                      The values of the input matrix at the validation indices.\n",
    "                      \n",
    "    test_values     : np.array, shape [n_test, ]\n",
    "                      The values of the input matrix at the test indices.\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    ### YOUR CODE HERE ###\n",
    "    #matrix = matrix[0]\n",
    "    nontrain_ind = np.random.choice(matrix.shape[0], size=(n_validation+n_test), replace=False)\n",
    "\n",
    "    nontrain_mat = matrix[nontrain_ind, :]\n",
    "    val_rows = np.random.choice(nontrain_ind.shape[0], size=n_validation, replace=False)\n",
    "    val_mat = nontrain_mat[val_rows, :]\n",
    "    val_ind = nontrain_ind[val_rows]\n",
    "    \n",
    "    test_mat = delete_rows_csr1(nontrain_mat.tocsr(), val_rows)\n",
    "    test_ind = np.delete(nontrain_ind,val_rows,0)\n",
    "    test_idx = []\n",
    "#     t = test_mat.nonzero()\n",
    "#     test_idx = list(zip(t[0],t[1]))\n",
    "    test_values = []\n",
    "    \n",
    "    for i in range(len(test_ind)):\n",
    "        for j in range(matrix.shape[1]):\n",
    "            #print(test_ind[i])\n",
    "            test_idx.append((test_ind[i], j))\n",
    "\n",
    "    test_idx = tuple(test_idx)\n",
    "    \n",
    "    for i in range(len(test_idx)):\n",
    "        test_values.append(matrix[test_idx[i]])\n",
    "    test_values = np.array(test_values)\n",
    "\n",
    "    val_idx = []\n",
    "    for i in range(len(val_ind)):\n",
    "        for j in range(matrix.shape[1]):\n",
    "            val_idx.append((val_ind[i], j))\n",
    "\n",
    "    val_idx = tuple(val_idx)\n",
    "#     v = val_mat.nonzero()\n",
    "#     val_idx = list(zip(v[0],v[1]))\n",
    "    val_values = []\n",
    "    for i in range(len(val_idx)):\n",
    "        val_values.append(matrix[val_idx[i]])\n",
    "    val_values = np.array(val_values)\n",
    "\n",
    "    matrix[nontrain_ind, :] = 0\n",
    "    \n",
    "    matrix_split = sp.csr_matrix(matrix)\n",
    "\n",
    "    matrix_split.eliminate_zeros()\n",
    "    return matrix_split, val_idx, test_idx, val_values, test_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape before: (337867, 5899)\n",
      "Shape after: (11275, 3531)\n"
     ]
    }
   ],
   "source": [
    "M = ratings_matrix\n",
    "M = cold_start_preprocessing(M, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\User\\Anaconda3\\lib\\site-packages\\scipy\\sparse\\compressed.py:774: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.\n",
      "  SparseEfficiencyWarning)\n"
     ]
    }
   ],
   "source": [
    "from copy import deepcopy\n",
    "n_validation = 200\n",
    "n_test = 200\n",
    "# Split data\n",
    "M1 = deepcopy(M[0])\n",
    "M_train, val_idx, test_idx, val_values, test_values = split_data(M[0], n_validation, n_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "validation_matrix = np.array([M1[lst] for lst in val_idx])\n",
    "validation_matrix = np.reshape(validation_matrix,(n_validation,M1.shape[1]))\n",
    "test_matrix = np.array([M1[lst] for lst in test_idx])\n",
    "test_matrix = np.reshape(test_matrix,(n_test,M1.shape[1]))\n",
    "#print(validation_matrix.shape,test_matrix.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\User\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:22: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    }
   ],
   "source": [
    "# Store away the nonzero indices of M before subtracting the row means.\n",
    "v = M[0].nonzero()\n",
    "nonzero_indices =  list(zip(v[0],v[1]))\n",
    "# Remove user means.\n",
    "M_train = sp.csr_matrix(M_train)\n",
    "M_shifted, train_user_means = shift_user_mean(M_train)\n",
    "# Apply the same shift to the validation and test data.\n",
    "validation_matrix = sp.csr_matrix(validation_matrix)\n",
    "test_matrix = sp.csr_matrix(test_matrix)\n",
    "val_values_shifted,val_user_means = shift_user_mean(validation_matrix)\n",
    "test_values_shifted,test_user_means = shift_user_mean(test_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Alternating optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the first step, we will approach the problem via alternating optimization, as learned in the lecture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def initialize_Q_P(matrix, k, init='random'):\n",
    "    \"\"\"\n",
    "    Initialize the matrices Q and P for a latent factor model.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    matrix : sp.spmatrix, shape [N, D]\n",
    "             The matrix to be factorized.\n",
    "    k      : int\n",
    "             The number of latent dimensions.\n",
    "    init   : str in ['svd', 'random'], default: 'random'\n",
    "             The initialization strategy. 'svd' means that we use SVD to initialize P and Q, 'random' means we initialize\n",
    "             the entries in P and Q randomly in the interval [0, 1).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Q : np.array, shape [N, k]\n",
    "        The initialized matrix Q of a latent factor model.\n",
    "\n",
    "    P : np.array, shape [k, D]\n",
    "        The initialized matrix P of a latent factor model.\n",
    "    \"\"\"\n",
    "\n",
    "    \n",
    "    if init == 'svd':\n",
    "        Q, s, Vh = np.linalg.svd(matrix, full_matrices=False)\n",
    "        P = np.dot(Q,np.dot(np.diag(s), Vh))\n",
    "    elif init == 'random':\n",
    "    ### YOUR CODE HERE ###\n",
    "        Q = np.random.normal(scale = 1/k,size=(matrix.shape[0], k))\n",
    "        P = np.random.normal(scale = 1/k,size=(k, matrix.shape[1]))\n",
    "    else:\n",
    "        raise ValueError\n",
    "        \n",
    "    assert Q.shape == (matrix.shape[0], k)\n",
    "    assert P.shape == (k, matrix.shape[1])\n",
    "    return Q, P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def latent_factor_alternating_optimization(M, non_zero_idx, k, val_idx, val_values,\n",
    "                                           reg_lambda, max_steps=100, init='random',\n",
    "                                           log_every=1, patience=10, eval_every=1):\n",
    "    \"\"\"\n",
    "    Perform matrix factorization using alternating optimization. Training is done via patience,\n",
    "    i.e. we stop training after we observe no improvement on the validation loss for a certain\n",
    "    amount of training steps. We then return the best values for Q and P oberved during training.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    M                 : sp.spmatrix, shape [N, D]\n",
    "                        The input matrix to be factorized.\n",
    "                      \n",
    "    non_zero_idx      : np.array, shape [nnz, 2]\n",
    "                        The indices of the non-zero entries of the un-shifted matrix to be factorized. \n",
    "                        nnz refers to the number of non-zero entries. Note that this may be different\n",
    "                        from the number of non-zero entries in the input matrix M, e.g. in the case\n",
    "                        that all ratings by a user have the same value.\n",
    "    \n",
    "    k                 : int\n",
    "                        The latent factor dimension.\n",
    "    \n",
    "    val_idx           : tuple, shape [2, n_validation]\n",
    "                        Tuple of the validation set indices.\n",
    "                        n_validation refers to the size of the validation set.\n",
    "                      \n",
    "    val_values        : np.array, shape [n_validation, ]\n",
    "                        The values in the validation set.\n",
    "                      \n",
    "    reg_lambda        : float\n",
    "                        The regularization strength.\n",
    "                      \n",
    "    max_steps         : int, optional, default: 100\n",
    "                        Maximum number of training steps. Note that we will stop early if we observe\n",
    "                        no improvement on the validation error for a specified number of steps\n",
    "                        (see \"patience\" for details).\n",
    "                      \n",
    "    init              : str in ['random', 'svd'], default 'random'\n",
    "                        The initialization strategy for P and Q. See function initialize_Q_P for details.\n",
    "    \n",
    "    log_every         : int, optional, default: 1\n",
    "                        Log the training status every X iterations.\n",
    "                    \n",
    "    patience          : int, optional, default: 10\n",
    "                        Stop training after we observe no improvement of the validation loss for X evaluation\n",
    "                        iterations (see eval_every for details). After we stop training, we restore the best \n",
    "                        observed values for Q and P (based on the validation loss) and return them.\n",
    "                      \n",
    "    eval_every        : int, optional, default: 1\n",
    "                        Evaluate the training and validation loss every X steps. If we observe no improvement\n",
    "                        of the validation error, we decrease our patience by 1, else we reset it to *patience*.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    best_Q            : np.array, shape [N, k]\n",
    "                        Best value for Q (based on validation loss) observed during training\n",
    "                      \n",
    "    best_P            : np.array, shape [k, D]\n",
    "                        Best value for P (based on validation loss) observed during training\n",
    "                      \n",
    "    validation_losses : list of floats\n",
    "                        Validation loss for every evaluation iteration, can be used for plotting the validation\n",
    "                        loss over time.\n",
    "                        \n",
    "    train_losses      : list of floats\n",
    "                        Training loss for every evaluation iteration, can be used for plotting the training\n",
    "                        loss over time.                     \n",
    "    \n",
    "    converged_after   : int\n",
    "                        it - patience*eval_every, where it is the iteration in which patience hits 0,\n",
    "                        or -1 if we hit max_steps before converging. \n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    ### YOUR CODE HERE ###\n",
    "    \n",
    "    Q,P = initialize_Q_P(M, k, init='random')\n",
    "    print(Q.shape,P.shape)\n",
    "    v = M.nonzero()\n",
    "    nonzero_train_indices =  list(zip(v[0],v[1]))\n",
    "    M_train_binary = np.zeros((M.shape[0],M.shape[1]))\n",
    "    for lst in nonzero_train_indices:\n",
    "        M_train_binary[lst] = 1\n",
    "        \n",
    "    def get_error(M, Q, P, M_train_binary):\n",
    "        # This calculates the MSE of nonzero elements\n",
    "        return np.sum((M_train_binary * np.asarray(M - np.dot(Q, P))) ** 2) / np.sum(M_train_binary)\n",
    "    \n",
    "    \n",
    "    print (\"Starting Iterations\")\n",
    "    \n",
    "    MSE_list = []\n",
    "    \n",
    "    for iter in range(max_steps):\n",
    "        for i, Mi in enumerate(M_train_binary):\n",
    "            print (i)\n",
    "            \n",
    "            temp = np.dot(P, np.dot(np.diag(Mi), P.T)) + reg_lambda * np.eye(k)\n",
    "            #print(temp.shape)\n",
    "            TEMP2 = (np.diag(Mi)).shape\n",
    "            TEMP3 = M[i].T.shape\n",
    "            #print(TEMP2,TEMP3)\n",
    "            temp1 = np.dot(P, np.dot((np.diag(Mi)), M[i].T))\n",
    "            #print(temp1.shape)\n",
    "            Q[i] = np.linalg.solve(np.dot(P, np.dot(np.diag(Mi), P.T)) + reg_lambda * np.eye(k),\n",
    "                                       np.dot(P,np.dot(np.asmatrix(np.diag(Mi)), M[i].T))).T\n",
    "            print (\"Error after solving for Q:\", get_error(M, Q, P, M_train_binary))\n",
    "\n",
    "        for j, Mj in enumerate(M_train_binary.T):\n",
    "            print(j)\n",
    "            print((np.diag(Mj)).shape)\n",
    "            P[:,j] = list(np.linalg.solve(np.dot(Q.T, np.dot(np.diag(Mj), Q)) + reg_lambda * np.eye(k),\n",
    "                                     np.dot(Q.T, np.dot(np.diag(Mj), M[:, j]))))\n",
    "            print (\"Error after solving for P:\", get_error(M, Q, P, M_train_binary))\n",
    "\n",
    "        MSE_List.append(get_error(M, Q, P, M_train_binary))\n",
    "        print ('%sth iteration is complete...' % iter)\n",
    "\n",
    "    print (MSE_List)\n",
    "    \n",
    "    \n",
    "    #return best_Q, best_P, validation_losses, train_losses, converged_after\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Train the latent factor model with alternating optimization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a) Learn the optimal $P$ and $Q$ using alternating optimization. That is, during each iteration you first update $Q$ while having $P$ fixed and then vice versa. Run the alternating optimization algorithm with $k=100$ and $\\lambda=1$. Plot the training and validation losses over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(11275, 100) (100, 3531)\n",
      "Starting Iterations\n",
      "0\n",
      "Error after solving for User Matrix: 1.20150091507\n",
      "1\n",
      "Error after solving for User Matrix: 1.20149978438\n",
      "2\n",
      "Error after solving for User Matrix: 1.201498618\n",
      "3\n",
      "Error after solving for User Matrix: 1.201498618\n",
      "4\n",
      "Error after solving for User Matrix: 1.20149811623\n",
      "5\n",
      "Error after solving for User Matrix: 1.20149771081\n",
      "6\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-54-32ae7ebfccdb>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      3\u001b[0m                                                                            \u001b[0mval_values\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mval_values_shifted\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m                                                                            \u001b[0mreg_lambda\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0minit\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'random'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 5\u001b[1;33m                                                                            max_steps=100, patience=10)\n\u001b[0m",
      "\u001b[1;32m<ipython-input-49-e9001a7fd7a8>\u001b[0m in \u001b[0;36mlatent_factor_alternating_optimization\u001b[1;34m(M, non_zero_idx, k, val_idx, val_values, reg_lambda, max_steps, init, log_every, patience, eval_every)\u001b[0m\n\u001b[0;32m    105\u001b[0m             Q[i] = np.linalg.solve(np.dot(P, np.dot(np.diag(Mi), P.T)) + reg_lambda * np.eye(k),\n\u001b[0;32m    106\u001b[0m                                        np.dot(P,np.dot(np.asmatrix(np.diag(Mi)), M[i].T))).T\n\u001b[1;32m--> 107\u001b[1;33m             \u001b[0mprint\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;34m\"Error after solving for User Matrix:\"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mget_error\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mM\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mQ\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mP\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mM_train_binary\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    108\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    109\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mMj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mM_train_binary\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-49-e9001a7fd7a8>\u001b[0m in \u001b[0;36mget_error\u001b[1;34m(M, Q, P, M_train_binary)\u001b[0m\n\u001b[0;32m     85\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mget_error\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mM\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mQ\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mP\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mM_train_binary\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     86\u001b[0m         \u001b[1;31m# This calculates the MSE of nonzero elements\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 87\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mM_train_binary\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mM\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mQ\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mP\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m**\u001b[0m \u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m/\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mM_train_binary\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     88\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     89\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py\u001b[0m in \u001b[0;36msum\u001b[1;34m(a, axis, dtype, out, keepdims)\u001b[0m\n\u001b[0;32m   1832\u001b[0m             \u001b[1;32mreturn\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0maxis\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mout\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1833\u001b[0m     return _methods._sum(a, axis=axis, dtype=dtype,\n\u001b[1;32m-> 1834\u001b[1;33m                          out=out, **kwargs)\n\u001b[0m\u001b[0;32m   1835\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1836\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\numpy\\core\\_methods.py\u001b[0m in \u001b[0;36m_sum\u001b[1;34m(a, axis, dtype, out, keepdims)\u001b[0m\n\u001b[0;32m     30\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     31\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_sum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 32\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0mumr_sum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     33\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     34\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_prod\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "latent_factor_alternating_optimization(M_shifted, nonzero_indices, \n",
    "                                                                           k=100, val_idx=val_idx,\n",
    "                                                                           val_values=val_values_shifted, \n",
    "                                                                           reg_lambda=1, init='random',\n",
    "                                                                           max_steps=100, patience=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the validation and training losses over (training) time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "### YOUR PLOTTING CODE HERE ###\n",
    "np.linalg.solve(np.dot(P, np.dot(np.diag(Mi), P.T)) + reg_lambda * np.eye(k),\n",
    "                                       np.dot(P, np.dot(np.asmatrix(np.diag(M_train_binary)).T, M[i]))).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b) (**Optional**): Try some different latent dimensions $k$ in the range [5, 100]. What do you observe (convergence time, final training/validation losses)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Q_a_2, P_a_2, val_l_a_2, tr_l_a_2, conv_a_2 = latent_factor_alternating_optimization(M_shifted, nonzero_indices, \n",
    "                                                                           k=20, val_idx=val_idx,\n",
    "                                                                           val_values=val_values_shifted, \n",
    "                                                                           reg_lambda=0.1, init='random',\n",
    "                                                                           max_steps=100, patience=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the validation and training losses over (training) time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "### YOUR PLOTTING CODE HERE ###\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Latent factorization using gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now use gradient descent to factorize our ratings matrix. We will try both (mini-) batch and stochastic gradient descent. You can use the following equations for your implementation.\n",
    "\n",
    "Recall that the objective function (loss) we wanted to optimize was:\n",
    "$$\n",
    "\\mathcal{L} = \\min_{P, Q} \\sum_{(x, i) \\in W} (r_{xi} - \\mathbf{q}_i^T\\mathbf{p}_x)^2 + \\lambda_1\\sum_x{\\left\\lVert \\mathbf{p}_x  \\right\\rVert}^2 + \\lambda_2\\sum_i {\\left\\lVert\\mathbf{q}_i  \\right\\rVert}^2\n",
    "$$\n",
    "\n",
    "where $W$ is the set of $(x, i)$ pairs for which $r_{xi}$ is known (in this case our known play counts). Here we have also introduced two regularization terms to help us with overfitting where $\\lambda_1$ and $\\lambda_2$ are hyper-parameters that control the strength of the regularization.\n",
    "\n",
    "Naturally optimizing with gradient descent involves computing the gradient of the loss function $\\mathcal{L}$ w.r.t. to the parameters. To help you solve the task we provide the following:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial ((r_{xi} - \\mathbf{q}_i^T\\mathbf{p}_x)^2)}{\\partial \\mathbf{p}_x} = -2(r_{xi} - \\mathbf{q}_i^T\\mathbf{p}_x)\\mathbf{q}_i\\;, ~~~\n",
    "\\frac{\\partial ((r_{xi} - \\mathbf{q}_i^T\\mathbf{p}_x)^2)}{\\partial \\mathbf{q}_i} = -2(r_{xi} - \\mathbf{q}_i^T\\mathbf{p}_x)\\mathbf{p}_x \n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial(\\lambda_1{\\left\\lVert \\mathbf{p}_x \\right\\rVert}^2)}{\\partial \\mathbf{p}_x} = 2 \\lambda_1 \\mathbf{p_x} \\;, ~~~\n",
    "\\frac{\\partial(\\lambda_2{\\left\\lVert \\mathbf{q}_i \\right\\rVert}^2)}{\\partial \\mathbf{q}_i} = 2 \\lambda_2 \\mathbf{q_i}\n",
    "$$\n",
    "\n",
    "**Hint**: You have to carefully consider how to combine the given partial gradients depending\n",
    "on which variants of gradient descent you are using.  \n",
    "**Hint 2**: It may be useful to scale the updates to $P$ and $Q$ by $\\frac{1}{batch\\_size}$ (in the case of full-sweep updates, this would be $\\frac{1}{n\\_users}$ for $Q$ and $\\frac{1}{n\\_restaurants}$ for $P$).\n",
    "\n",
    "\n",
    "For each of the gradients descent variants you try report and compare the following:\n",
    "* How many iterations do you need for convergence.\n",
    "* Plot the loss (y axis) for each iteration (x axis).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def latent_factor_gradient_descent(M, non_zero_idx, k, val_idx, val_values, \n",
    "                                   reg_lambda, learning_rate, batch_size=-1,\n",
    "                                   max_steps=50000, init='random',\n",
    "                                   log_every=1000, patience=20,\n",
    "                                   eval_every=50):\n",
    "    \"\"\"\n",
    "    Perform matrix factorization using gradient descent. Training is done via patience,\n",
    "    i.e. we stop training after we observe no improvement on the validation loss for a certain\n",
    "    amount of training steps. We then return the best values for Q and P oberved during training.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    M                 : sp.spmatrix, shape [N, D]\n",
    "                        The input matrix to be factorized.\n",
    "                      \n",
    "    non_zero_idx      : np.array, shape [nnz, 2]\n",
    "                        The indices of the non-zero entries of the un-shifted matrix to be factorized. \n",
    "                        nnz refers to the number of non-zero entries. Note that this may be different\n",
    "                        from the number of non-zero entries in the input matrix M, e.g. in the case\n",
    "                        that all ratings by a user have the same value.\n",
    "    \n",
    "    k                 : int\n",
    "                        The latent factor dimension.\n",
    "    \n",
    "    val_idx           : tuple, shape [2, n_validation]\n",
    "                        Tuple of the validation set indices.\n",
    "                        n_validation refers to the size of the validation set.\n",
    "                      \n",
    "    val_values        : np.array, shape [n_validation, ]\n",
    "                        The values in the validation set.\n",
    "                      \n",
    "    reg_lambda        : float\n",
    "                        The regularization strength.\n",
    "\n",
    "    learning_rate     : float\n",
    "                        Step size of the gradient descent updates.\n",
    "                        \n",
    "    batch_size        : int, optional, default: -1\n",
    "                        (Mini-) batch size. -1 means we perform standard full-sweep gradient descent.\n",
    "                        If the batch size is >0, use mini batches of this given size.\n",
    "                        \n",
    "    max_steps         : int, optional, default: 100\n",
    "                        Maximum number of training steps. Note that we will stop early if we observe\n",
    "                        no improvement on the validation error for a specified number of steps\n",
    "                        (see \"patience\" for details).\n",
    "                      \n",
    "    init              : str in ['random', 'svd'], default 'random'\n",
    "                        The initialization strategy for P and Q. See function initialize_Q_P for details.\n",
    "    \n",
    "    log_every         : int, optional, default: 1\n",
    "                        Log the training status every X iterations.\n",
    "                    \n",
    "    patience          : int, optional, default: 10\n",
    "                        Stop training after we observe no improvement of the validation loss for X evaluation\n",
    "                        iterations (see eval_every for details). After we stop training, we restore the best \n",
    "                        observed values for Q and P (based on the validation loss) and return them.\n",
    "                      \n",
    "    eval_every        : int, optional, default: 1\n",
    "                        Evaluate the training and validation loss every X steps. If we observe no improvement\n",
    "                        of the validation error, we decrease our patience by 1, else we reset it to *patience*.\n",
    "                        \n",
    "    Returns\n",
    "    -------\n",
    "    best_Q            : np.array, shape [N, k]\n",
    "                        Best value for Q (based on validation loss) observed during training\n",
    "                      \n",
    "    best_P            : np.array, shape [k, D]\n",
    "                        Best value for P (based on validation loss) observed during training\n",
    "                      \n",
    "    validation_losses : list of floats\n",
    "                        Validation loss for every evaluation iteration, can be used for plotting the validation\n",
    "                        loss over time.\n",
    "                        \n",
    "    train_losses      : list of floats\n",
    "                        Training loss for every evaluation iteration, can be used for plotting the training\n",
    "                        loss over time.                     \n",
    "    \n",
    "    converged_after   : int\n",
    "                        it - patience*eval_every, where it is the iteration in which patience hits 0,\n",
    "                        or -1 if we hit max_steps before converging. \n",
    "\n",
    "    \"\"\"\n",
    "       \n",
    "    ### YOUR CODE HERE ###\n",
    "    num_users, num_restaurants = M.shape \n",
    "    Q, s,P = svds(M, k)\n",
    "    print(P.shape)\n",
    "    print(Q.shape)\n",
    "    #Q = Q/num_users\n",
    "    #P= P/num_restaurants\n",
    "    b_q = np.zeros(num_users)\n",
    "    b_p = np.zeros(num_restaurants)\n",
    "\n",
    "    train_losses = []\n",
    "    for steps in range(max_steps):\n",
    "        #np.random.shuffle(samples)\n",
    "        prediction1 = np.dot(Q,P)\n",
    "        e = (M - prediction1)\n",
    "        print((reg_lambda*P).shape)\n",
    "        Q +=learning_rate*2*(-np.dot(e,P.T) + (reg_lambda*Q))   \n",
    "        P += learning_rate*2*(-np.dot(Q.T,e)+(reg_lambda*P))\n",
    "        predicted = np.dot(Q,P)\n",
    "        error = np.sum(np.square(np.subtract(M, predicted))) + reg_lambda * np.linalg.norm(P) + reg_lambda * np.linalg.norm(Q)\n",
    "        mse = np.sqrt(error)\n",
    "        train_losses.append((steps, mse))\n",
    "        print(\"Iteration: %d ; error = %.4f\" % (steps, mse))\n",
    "        converged_after = steps                \n",
    "    best_Q = Q\n",
    "    best_P = P\n",
    "    validation_losses = []\n",
    "    return best_Q, best_P, validation_losses, train_losses, converged_after\n",
    "            \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Train the latent factor model with alternating optimization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a) Learn the optimal $P$ and $Q$ using standard gradient descent. That is, during each iteration you have to use all of the training examples and update $Q$ and $P$ for all users and songs at once. Try the algorithm with $k=30$, $\\lambda=1$, and learning rate of 0.1. Initialize $Q$ and $P$ with SVD.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(30, 3531)\n",
      "(11275, 30)\n",
      "(30, 3531)\n",
      "Iteration: 0 ; error = 1142.7737\n",
      "(30, 3531)\n",
      "Iteration: 1 ; error = 6294126537.8894\n",
      "(30, 3531)\n",
      "Iteration: 2 ; error = 1386197863961803928718641535171203956736.0000\n",
      "(30, 3531)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\User\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:103: RuntimeWarning: overflow encountered in square\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration: 3 ; error = inf\n",
      "(30, 3531)\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-53-bfb997e21f77>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      5\u001b[0m                                                                                                    \u001b[0minit\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'svd'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbatch_size\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      6\u001b[0m                                                                                                    \u001b[0mmax_steps\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m10000\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlog_every\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m20\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 7\u001b[1;33m                                                                                                    eval_every=20)\n\u001b[0m",
      "\u001b[1;32m<ipython-input-51-fa44657aa80c>\u001b[0m in \u001b[0;36mlatent_factor_gradient_descent\u001b[1;34m(M, non_zero_idx, k, val_idx, val_values, reg_lambda, learning_rate, batch_size, max_steps, init, log_every, patience, eval_every)\u001b[0m\n\u001b[0;32m     98\u001b[0m         \u001b[0me\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mM\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mprediction1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     99\u001b[0m         \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mreg_lambda\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mP\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 100\u001b[1;33m         \u001b[0mQ\u001b[0m \u001b[1;33m+=\u001b[0m\u001b[0mlearning_rate\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0me\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mP\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mreg_lambda\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mQ\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    101\u001b[0m         \u001b[0mP\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0mlearning_rate\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mQ\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0me\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mreg_lambda\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mP\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    102\u001b[0m         \u001b[0mpredicted\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mQ\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mP\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "Q_g_sweep, P_g_sweep, val_l_g_sweep, tr_l_g_sweep, conv_g_sweep =  latent_factor_gradient_descent(M_shifted, nonzero_indices, \n",
    "                                                                                                   k=30, val_idx=val_idx,\n",
    "                                                                                                   val_values=val_values_shifted, \n",
    "                                                                                                   reg_lambda=1, learning_rate=1e-1,\n",
    "                                                                                                   init='svd', batch_size=-1,\n",
    "                                                                                                   max_steps=10000, log_every=20, \n",
    "                                                                                                   eval_every=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the validation and training losses over (training) time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "### YOUR PLOTTING CODE HERE ###\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b) Learn the optimal $P$ and $Q$ using the original stochastic gradient descent (mini-batches of size 1). That is, during each iteration you sample a single random training example $r_{xi}$ and update only the respective affected parameters $\\mathbf{p_x}$ and $\\mathbf{q}_i$. Set the learning rate to 0.01 and keep the other parameters as in a)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Q_g_st, P_g_st, val_l_g_st, tr_l_g_st, conv_g_st = latent_factor_gradient_descent(M_shifted, nonzero_indices, \n",
    "                                                                                   k=30, val_idx=val_idx,\n",
    "                                                                                   val_values=val_values_shifted, \n",
    "                                                                                   reg_lambda=1, learning_rate=1e-2,\n",
    "                                                                                   init='svd', batch_size=1,\n",
    "                                                                                   max_steps=20000, log_every=500, \n",
    "                                                                                   eval_every=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the validation and training losses over (training) time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "### YOUR PLOTTING CODE HERE ###\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "c) (**Optional**) Learn the optimal $P$ and $Q$ similarly to b) this time using larger mini-batches of size $S$, e.g. 32."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Q_g_mb, P_g_mb, val_l_g_mb, tr_l_g_mb, conv_g_mb = latent_factor_gradient_descent(M_shifted, nonzero_indices, \n",
    "                                                                                   k=30, val_idx=val_idx,\n",
    "                                                                                   val_values=val_values_shifted, \n",
    "                                                                                   reg_lambda=1, learning_rate=1e-1,\n",
    "                                                                                   init='svd', batch_size=32,\n",
    "                                                                                   max_steps=10000, log_every=100, \n",
    "                                                                                   eval_every=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the validation and training losses over (training) time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "### YOUR PLOTTING CODE HERE ###\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.5 Hyperparameter search"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Machine learning models are often heavily dependent on the hyperparameter settings, e.g. the learning rate. Here, we will try a simple random search to find good values of the latent factor dimension $k$, the batch size, learning rate, and regularization.  \n",
    "\n",
    "### Tasks:\n",
    "\n",
    "Perform a hyperparameter search to find good values for the batch size, lambda, learning rate, and latent dimension. \n",
    "\n",
    "* For the batch size, evaluate all values in [1, 32, 512, -1] (-1 corresponds to full-sweep gradient descent).\n",
    "* For $\\lambda$, randomly sample three values in the interval [0, 1).\n",
    "* For the learning rate, evaluate all values in [1, 0.1, 0.01].\n",
    "* For the latent dimension, uniformly sample three values in the interval [5,30].\n",
    "\n",
    "Perform an exhaustive search among all combinations of these values;\n",
    "\n",
    "**Hint**: This may take a while to compute. **You don't have to wait for all the models to train** -- simply use \"dummy\" code instead of actual model training (or let it train, e.g., for only one iteration) if you don't want to wait. Note that the signature of this dummy code has to match the function 'latent_factor_gradient_descent' so that we could simply plug in the actual function.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def parameter_search(M_train, val_idx, val_values):\n",
    "    \"\"\"\n",
    "    Hyperparameter search using random search.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    M_train     : sp.spmatrix, shape [N, D]\n",
    "                  Input sparse matrix where the user means have not\n",
    "                  been subtracted yet. \n",
    "                  \n",
    "    val_idx     : tuple, shape [2, n_validation]\n",
    "                  The indices used for validation, where n_validation\n",
    "                  is the size of the validation set.\n",
    "                  \n",
    "    val_values  : np.array, shape [n_validation, ]\n",
    "                  Validation set values, where n_validation is the\n",
    "                  size of the validation set.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    best_conf   : tuple, (batch_size, lambda, learning_rate, latent_dimension)\n",
    "                  The best-performing hyperparameters.\n",
    "                  \n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    ### YOUR CODE HERE ###\n",
    "        \n",
    "    print(\"Best configuration is {}\").format(best_conf)\n",
    "    return best_conf\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "best_configuration = parameter_search(M_train, val_idx, val_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Output the best hyperparameter optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "### YOUR CODE HERE ###\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Comparison of gradient descent and alternating optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After training the latent factor model with both alternating optimization and gradient descent, we now compare their results on the training, validation, and test set.\n",
    "\n",
    "### Tasks\n",
    "\n",
    "* Compare the root mean square errors (RMSE) for the training, validation, and test sets different settings of $k$ for both alternating optimization and gradient descent. What do you observe?\n",
    "* Compare the test RMSE for the alternating optimization model and the gradient descent model. Which performs better?\n",
    "* Plot the predicted ratings\n",
    "\n",
    "**Hint**: The output values and plots below are the ones we got when testing this sheet. Yours may be different, but if your validation or test RMSE values are larger than 1.5 or 2, it is likely that you have a bug in your implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "### YOUR CODE HERE ###\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "### YOUR CODE HERE ###\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plots: Prediction vs. ground truth ratings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "### YOUR PLOTTING CODE HERE ###\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "### YOUR PLOTTING CODE HERE ###\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "### YOUR PLOTTING CODE HERE ###\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "### YOUR PLOTTING CODE HERE ###\n"
   ]
  }
 ],
 "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
}