recommender system.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recommender System\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Recommender system is a system widely used in Amazon, Netflix to predict user rating for a given item. Usually the system involves collaborative filtering, content-based filtering, session-based filtering or even a mixture.\n",
    "\n",
    "* Collaborative filtering seeks connections across different users and items to predict the rating. It is more related to unsupervised learning. Even inside collaborative filtering, there are a few different techniques such as model-based (matrix factorization, latent variables), memory-based (KNN, KNN with z score) and even a hybrid of both. Collaborative filtering will be the main focus of this script. \n",
    "* Content-based filtering collects user and item profile. Based upon the features in the profile, it forecasts the user preference towards different items. It is more related to supervised learning.\n",
    "* Session-based filtering monitors the user interaction within a session to create recommendations. It is quite helpful to navigate through the cold start problem where a new user has not much available information for modelling. It is more related to deep learning, in particular, Recurrent Neural Network.\n",
    "\n",
    "In a recommender system, not every customer rates every item so it effectively forms a partially filled customer vs item matrix. To recommend anything to the existing customer, a fully filled matrix must be inferred from the dataset. This falls into the spectrum of matrix completion. The most popular sub-problem of matrix completion is to find a low rank matrix via convex optimization (Candès and Recht, 2008). The sub-problem assumes there must be latent variables influencing the users and the items in the matrix. Thus, the matrix must be low rank.\n",
    "\n",
    "Reference to image recovery style matrix completion\n",
    "\n",
    "https://github.com/je-suis-tm/machine-learning/blob/master/matrix%20completion.ipynb\n",
    "\n",
    "Reference to surprise library for recommender system (this script draws a lot of inspirations from it though much faster convergence due to the intensive usage of linear algebra rather than iterations)\n",
    "\n",
    "https://surprise.readthedocs.io/en/stable/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.chdir('K:/ecole/github/televerser/matrix completion')\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import numba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#data comes from the link below\n",
    "# http://files.grouplens.org/datasets/movielens/ml-100k/u.data\n",
    "# https://github.com/je-suis-tm/machine-learning/blob/master/data/movielens.csv\n",
    "df=pd.read_csv('movielens.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:2: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "#shrink data size for better performance\n",
    "df=df[df['user']<=100][df['item']<=500].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert array into matrix form\n",
    "matrix=df.pivot(index='item',columns='user',values='rating')\n",
    "arr=np.array(matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#filled elements\n",
    "known_values=list(zip(np.where(~np.isnan(arr))[0],\n",
    "        np.where(~np.isnan(arr))[1]))\n",
    "\n",
    "#randomly select 30% as testing set\n",
    "lucky_draw=set(np.random.choice(range(len(known_values)),\n",
    "                 size=int(len(known_values)*0.3),\n",
    "                 p=[1/len(known_values)]*len(known_values)))\n",
    "\n",
    "unlucky_draw=[i for i in range(len(known_values)) if i not in lucky_draw]\n",
    "testing_idx=(np.where(~np.isnan(arr))[0][list(lucky_draw)],\n",
    "np.where(~np.isnan(arr))[1][list(lucky_draw)])\n",
    "training_idx=(np.where(~np.isnan(arr))[0][list(unlucky_draw)],\n",
    "np.where(~np.isnan(arr))[1][list(unlucky_draw)])\n",
    "\n",
    "#train test split\n",
    "mask=np.ones(arr.shape)\n",
    "mask[testing_idx]=np.nan\n",
    "arr_train=np.multiply(arr,mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model-based\n",
    "\n",
    "#### Funk SVD\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Simon Funk developed a SVD-like latent factor model to win 3rd prize in 2006 Netflix problem. Conventionally people call it Funk SVD, although it is merely inspired by Singular Value Decomposition without explicitly using SVD. Funk SVD is really easy to implement and quick to converge with high accuracy. It is a type of collaborative filtering focusing on matrix factorization. The official optimization problem is formulated as below.\n",
    "\n",
    "$$ \\min_{p_*,q_*,b_*}\\,\\sum_{r_{ui}\\,\\in\\,\\mathcal{K}} \\left(r_{ui} - \\hat{r}_{ui} \\right)^2 + \\lambda( ||p_u||^2 + ||q_i||^2 + b_u^2 + b_i^2 )$$\n",
    "\n",
    "where \n",
    "\n",
    "${r}_{ui}$ denotes the rating of item $i$ by user $u$\n",
    "\n",
    "$\\hat{r}_{ui}$ denotes the estimated rating of item $i$ by user $u$, it can be decomposed into the form of $\\mu + b_u + b_i + q_i^Tp_u$\n",
    "\n",
    "$\\mathcal{K}$ denotes the observed user vs item matrix (partially filled)\n",
    "\n",
    "$\\mu$ denotes overall average rating\n",
    "\n",
    "$b_u$ denotes the deviation from overall average rating caused by user $u$, can be seen as user baseline where some users are more critical so they give lower ratings in general\n",
    "\n",
    "$b_i$ denotes the deviation from overall average rating caused by item $i$, can be seen as item baseline where some items are more appealing so they attract higher ratings in general\n",
    "\n",
    "$p_u$ denotes the latent factors of user $u$, in translation, user preference\n",
    "\n",
    "$q_i$ denotes the latent factors of item $i$, in translation, item attributes\n",
    "\n",
    "$\\lambda$ denotes LaGrange multiplier which is the coefficient of L2 penalty function\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "To explain a bit, Funk SVD assumes that both user and item are affected by $m$ number of latent factors ($p_u \\in \\mathbb{R}^{m}\\,\\&\\,q_i \\in \\mathbb{R}^{m}$). These latent factors could be percentage of action in the movie, number of tier 1 hollywood stars, etc. The rating $r_{ui}$ is merely the linear combination $q_i^Tp_u$ of user preference and item attributes. This decomposition is similar to SVD in the form of $U\\Sigma V^T$ without the eigenvalue diagonal $\\Sigma$ as scaling factors. Hence, the first bit of the objective function is an ordinary least square to minimize the sum of squared error between existing rating and estimated rating. Since the matrix is partially filled, the second bit of the objective function is L2 norm regularization on user preference $p_u$, item attributes $q_i$, user deviation $b_u$ and item deviation $b_i$ to avoid overfit problem.\n",
    "\n",
    "The actual optimization solver is inspired by gradient descent $\\theta:=\\theta-\\alpha\\frac{\\partial J(\\theta)}{\\partial \\theta}$. Applying the same logic to every unknown parameters $p_u$, $q_i$, $b_u$ and $b_i$, we end up with the following updates.\n",
    "\n",
    "$$b_u := b_u + \\alpha (\\epsilon_{ui} - \\lambda b_u)$$\n",
    "$$b_i := b_i + \\alpha (\\epsilon_{ui} - \\lambda b_i)$$\n",
    "$$p_u := p_u + \\alpha (\\epsilon_{ui} \\cdot q_i - \\lambda p_u)$$\n",
    "$$q_i := q_i + \\alpha (\\epsilon_{ui} \\cdot p_u - \\lambda q_i)$$\n",
    "\n",
    "where\n",
    "\n",
    "$\\epsilon_{ui}$ denotes the error between existing rating and estimated rating in the form of $r_{ui} - \\hat{r}_{ui}$\n",
    "\n",
    "$\\alpha$ denotes the learning rate of gradient descent which dictates how soon the solver reaches the local optima\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Reference to Singular Value Decomposition\n",
    "\n",
    "https://github.com/je-suis-tm/machine-learning/blob/master/principal%20component%20analysis.ipynb\n",
    "\n",
    "Reference to Gradient Descent\n",
    "\n",
    "https://github.com/je-suis-tm/machine-learning/blob/master/gradient%20descent.ipynb\n",
    "\n",
    "Reference to Simon Funk's blog\n",
    "\n",
    "https://sifter.org/~simon/journal/20061211.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use numba to dramatically boost the speed of linear algebra\n",
    "#this will be a lot faster than surprise library\n",
    "@numba.njit\n",
    "def funk_svd_epoch(arr,miu,b_u,b_i,p_u,q_i,alpha,lambda_):\n",
    "    \n",
    "    #initialize\n",
    "    error=0\n",
    "    \n",
    "    #only iterate known ratings\n",
    "    for i in range(arr.shape[0]):\n",
    "        for u in range(arr.shape[1]):\n",
    "            r_ui=arr[i,u]\n",
    "            \n",
    "            #another way to identify nan\n",
    "            #r_ui!=r_ui\n",
    "            if np.isnan(r_ui):\n",
    "                continue\n",
    "\n",
    "            #compute error\n",
    "            epsilon_ui=r_ui-miu-b_u[u]-b_i[i]-q_i[i].T@p_u[u]\n",
    "            error+=epsilon_ui**2\n",
    "\n",
    "            #update\n",
    "            b_u[u]+=alpha*(epsilon_ui-lambda_*b_u[u])\n",
    "            b_i[i]+=alpha*(epsilon_ui-lambda_*b_i[i])\n",
    "            p_u[u]+=alpha*(epsilon_ui*q_i[i]-lambda_*p_u[u])\n",
    "            q_i[i]+=alpha*(epsilon_ui*p_u[u]-lambda_*q_i[i])\n",
    "    \n",
    "    return error,b_u,b_i,p_u,q_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#svd inspired latent factor model by simon funk\n",
    "def funk_svd(arr,miu_init=None,b_u_init=[],b_i_init=[],\n",
    "             p_u_init=[],q_i_init=[],num_of_rank=40,\n",
    "             alpha=0.005,lambda_=0.02,tau=0.0001,\n",
    "             max_iter=20,diagnosis=True\n",
    "             ):\n",
    "\n",
    "    #initialize\n",
    "    stop=False\n",
    "    counter=0\n",
    "    sse=None\n",
    "    \n",
    "    #global mean\n",
    "    if not miu_init:       \n",
    "        miu=arr[~np.isnan(arr)].mean()\n",
    "    else:\n",
    "        miu=miu_init\n",
    "        \n",
    "    #user baseline\n",
    "    if len(b_u_init)==0:\n",
    "        b_u=np.zeros(arr.shape[1])\n",
    "    else:\n",
    "        b_u=b_u_init\n",
    "    \n",
    "    #item baseline\n",
    "    if len(b_i_init)==0:\n",
    "        b_i=np.zeros(arr.shape[0])\n",
    "    else:\n",
    "        b_i=b_i_init\n",
    "        \n",
    "    #user latent factors\n",
    "    if len(p_u_init)==0:\n",
    "        p_u=np.zeros((arr.shape[1],num_of_rank))\n",
    "        p_u.fill(0.1)\n",
    "    else:\n",
    "        p_u=p_u_init\n",
    "    \n",
    "    #item latent factors\n",
    "    if len(q_i_init)==0:\n",
    "        q_i=np.zeros((arr.shape[0],num_of_rank))\n",
    "        q_i.fill(0.1)\n",
    "    else:\n",
    "        q_i=q_i_init\n",
    "    \n",
    "    #gradient descent\n",
    "    while not stop:\n",
    "        \n",
    "        error,b_u,b_i,p_u,q_i=funk_svd_epoch(arr,miu,\n",
    "                                             b_u,b_i,\n",
    "                                             p_u,q_i,\n",
    "                                             alpha,lambda_)\n",
    "\n",
    "        counter+=1\n",
    "\n",
    "        #maximum number of epoch\n",
    "        if counter>=max_iter:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print('Not converged.',\n",
    "                      'Consider increase number of iterations or tolerance')\n",
    "                \n",
    "        #use sum of squared error to determine if converged\n",
    "        sse_prev=sse\n",
    "        sse=error\n",
    "        if sse_prev and abs(sse/sse_prev-1)<=tau:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print(f'{counter} iterations to reach convergence\\n')\n",
    "\n",
    "    return b_u,b_i,p_u,q_i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize\n",
    "num_of_latent_factors=40\n",
    "max_num_of_epoch=500\n",
    "learning_rate=0.01\n",
    "lagrange_multiplier=0.02\n",
    "tolerance=0.005"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30 iterations to reach convergence\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#funk svd\n",
    "b_u,b_i,p_u,q_i=funk_svd(arr_train,num_of_rank=num_of_latent_factors,\n",
    "         alpha=learning_rate,\n",
    "         lambda_=lagrange_multiplier,tau=tolerance,\n",
    "         max_iter=max_num_of_epoch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute global mean\n",
    "miu=arr_train[~np.isnan(arr_train)].mean()\n",
    "\n",
    "#matrix completion\n",
    "output=miu+np.repeat(\n",
    "            b_u.reshape(1,-1),\n",
    "            arr_train.shape[0],axis=0)+np.repeat(\n",
    "            b_i.reshape(-1,1),arr_train.shape[1],axis=1)+q_i@p_u.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Funk SVD Mean Squared Error: 0.927\n"
     ]
    }
   ],
   "source": [
    "#use mse as benchmark for comparison\n",
    "mse_funk_svd=np.square((\n",
    "    output-arr)[testing_idx]).sum()/len(arr[testing_idx])\n",
    "print('Funk SVD Mean Squared Error:',round(mse_funk_svd,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### SVD++\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Inspired by Funk SVD, more and more variations have been invented to enhance the prediction accuracy and the convergence speed. One of the most famous variation is SVD++ (Koren, 2008). It can be seen as an augmented Funk SVD with user implicit feedback. The accuracy is further improved at the cost of implicit feedback integration. The official optimization problem is formulated as below. The notations are almost the same with only a few extra from implicit feedback.\n",
    "\n",
    "$$ \\min_{p_*,q_*,b_*,y_*}\\,\\sum_{r_{ui}\\,\\in\\,\\mathcal{K}} \\left(r_{ui} - \\mu - b_u - b_i - q_i^T\\left(p_u+|N(u)|^{-\\frac{1}{2}} \\sum_{j \\in N(u)}y_j\\right) \\right)^2 + \\lambda\\left( ||p_u||^2 + ||q_i||^2 + \\sum_{j \\in N(u)}||y_j||^2 + b_u^2 + b_i^2 \\right)$$\n",
    "\n",
    "where \n",
    "\n",
    "$N(u)$ denotes the implicit preference list of user $u$\n",
    "\n",
    "$y_j$ denotes the item $j$ that user $u$ implicitly prefers\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "To explain a bit, SVD++ incorporates user implicit feedback into the model. Users may have watched a lot of movies on Netflix yet they choose to rate items in $N(u)$ for a particular reason. For instance, user $u$ may have watched six seasons of House of Cards but he strongly resents the last season. Hence, he gives one star rating to the last season. This kind of implicit feedback is not captured in any latent factors. Nonetheless, the implicit feedback $y_j$ is not independent of any latent factors. Since certain combination $q_i^Tp_u$ triggers the user response and amplifies the extreme rating, $y_j$ is included in both the objective function and the constraint to avoid overfit problem. FYI, $y_j$ can be seen as another latent factor vector with the exact dimension as $q_i$.\n",
    "\n",
    "The actual optimization solver is similar to Funk SVD with one more parameter $y_j$. \n",
    "\n",
    "$$b_u := b_u + \\alpha (\\epsilon_{ui} - \\lambda b_u)$$\n",
    "$$b_i := b_i + \\alpha (\\epsilon_{ui} - \\lambda b_i)$$\n",
    "$$p_u := p_u + \\alpha (\\epsilon_{ui} \\cdot q_i - \\lambda p_u)$$\n",
    "$$q_i := q_i + \\alpha (\\epsilon_{ui} \\cdot (p_u+|N(u)|^{-\\frac{1}{2}} \\sum_{j \\in N(u)}y_j) - \\lambda q_i)$$\n",
    "$$y_j := y_j + \\alpha (|N(u)|^{-\\frac{1}{2}} \\cdot \\epsilon_{ui} \\cdot q_i - \\lambda y_j)$$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Reference to the original paper of SVD++ (equation 15)\n",
    "\n",
    "https://people.engr.tamu.edu/huangrh/Spring16/papers_course/matrix_factorization.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use numba to dramatically boost the speed of linear algebra\n",
    "#this will be a lot faster than surprise library\n",
    "@numba.njit\n",
    "def svd_plus_plus_epoch(arr,miu,b_u,b_i,p_u,q_i,y_j,alpha,lambda_):\n",
    "    \n",
    "    #initialize\n",
    "    error=0\n",
    "    \n",
    "    #only iterate known ratings\n",
    "    for i in range(arr.shape[0]):\n",
    "        for u in range(arr.shape[1]):\n",
    "            r_ui=arr[i,u]\n",
    "            \n",
    "            #another way to identify nan\n",
    "            #r_ui!=r_ui\n",
    "            if np.isnan(r_ui):\n",
    "                continue\n",
    "                \n",
    "            #compute implicit feedback\n",
    "            N_u=np.where(~np.isnan(arr[:,u]))[0]\n",
    "            N_u_norm_sqrt=arr[:,u][~np.isnan(arr[:,u])].shape[0]**0.5\n",
    "            feedback=(y_j[N_u,:]/N_u_norm_sqrt).sum(axis=0)\n",
    "\n",
    "            #compute error\n",
    "            epsilon_ui=(r_ui-miu-b_u[u]-b_i[i]-q_i[i].T@(\n",
    "                feedback+p_u[u]).reshape(-1,1)).item()\n",
    "            error+=epsilon_ui**2\n",
    "\n",
    "            #update\n",
    "            b_u[u]+=alpha*(epsilon_ui-lambda_*b_u[u])\n",
    "            b_i[i]+=alpha*(epsilon_ui-lambda_*b_i[i])\n",
    "            p_u[u]+=alpha*(epsilon_ui*q_i[i]-lambda_*p_u[u])\n",
    "            q_i[i]+=alpha*(epsilon_ui*(p_u[u]+feedback)-lambda_*q_i[i])\n",
    "            y_j[N_u]+=alpha*(epsilon_ui*q_i[N_u]/N_u_norm_sqrt-lambda_*y_j[N_u])\n",
    "               \n",
    "    return error,b_u,b_i,p_u,q_i,y_j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#svdpp latent factor model\n",
    "def svd_plus_plus(arr,miu_init=None,b_u_init=[],b_i_init=[],\n",
    "             p_u_init=[],q_i_init=[],y_j_init=[],num_of_rank=40,\n",
    "             alpha=0.005,lambda_=0.02,tau=0.0001,\n",
    "             max_iter=20,diagnosis=True\n",
    "             ):\n",
    "\n",
    "    #initialize\n",
    "    stop=False\n",
    "    counter=0\n",
    "    sse=None\n",
    "    \n",
    "    #global mean\n",
    "    if not miu_init:       \n",
    "        miu=arr[~np.isnan(arr)].mean()\n",
    "    else:\n",
    "        miu=miu_init\n",
    "        \n",
    "    #user baseline\n",
    "    if len(b_u_init)==0:\n",
    "        b_u=np.zeros(arr.shape[1])\n",
    "    else:\n",
    "        b_u=b_u_init\n",
    "    \n",
    "    #item baseline\n",
    "    if len(b_i_init)==0:\n",
    "        b_i=np.zeros(arr.shape[0])\n",
    "    else:\n",
    "        b_i=b_i_init\n",
    "        \n",
    "    #user latent factors\n",
    "    if len(p_u_init)==0:\n",
    "        p_u=np.zeros((arr.shape[1],num_of_rank))\n",
    "        p_u.fill(0.1)\n",
    "    else:\n",
    "        p_u=p_u_init\n",
    "    \n",
    "    #item latent factors\n",
    "    if len(q_i_init)==0:\n",
    "        q_i=np.zeros((arr.shape[0],num_of_rank))\n",
    "        q_i.fill(0.1)\n",
    "    else:\n",
    "        q_i=q_i_init\n",
    "        \n",
    "    #user implicit feedback\n",
    "    if len(y_j_init)==0:\n",
    "        y_j=np.zeros((arr.shape[0],num_of_rank))\n",
    "        y_j.fill(0.1)\n",
    "    else:\n",
    "        y_j=y_j_init\n",
    "    \n",
    "    #gradient descent\n",
    "    while not stop:\n",
    "        \n",
    "        error,b_u,b_i,p_u,q_i,y_j=svd_plus_plus_epoch(\n",
    "            arr,miu,b_u,b_i,p_u,q_i,y_j,alpha,lambda_)\n",
    "\n",
    "        counter+=1\n",
    "        \n",
    "        #maximum number of epoch\n",
    "        if counter>=max_iter:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print('Not converged.',\n",
    "                      'Consider increase number of iterations or tolerance')\n",
    "                \n",
    "        #use sum of squared error to determine if converged\n",
    "        sse_prev=sse\n",
    "        sse=error\n",
    "        if sse_prev and abs(sse/sse_prev-1)<=tau:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print(f'{counter} iterations to reach convergence\\n')\n",
    "\n",
    "    return b_u,b_i,p_u,q_i,y_j"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "17 iterations to reach convergence\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#svdpp\n",
    "b_u,b_i,p_u,q_i,y_j=svd_plus_plus(arr_train,\n",
    "         num_of_rank=num_of_latent_factors,\n",
    "         alpha=learning_rate,\n",
    "         lambda_=lagrange_multiplier,tau=tolerance,\n",
    "         max_iter=max_num_of_epoch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute implicit feedback\n",
    "feedback=[]\n",
    "for u in range(arr_train.shape[1]):\n",
    "    N_u=np.where(~np.isnan(arr_train[:,u]))[0]\n",
    "    N_u_norm_sqrt=arr_train[:,u][~np.isnan(arr_train[:,u])].shape[0]**0.5\n",
    "    feedback.append((y_j[N_u,:]/N_u_norm_sqrt).sum(axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix completion\n",
    "output=miu+np.repeat(\n",
    "            b_u.reshape(1,-1),\n",
    "            arr_train.shape[0],axis=0)+np.repeat(\n",
    "            b_i.reshape(-1,1),\n",
    "    arr_train.shape[1],axis=1)+q_i@(p_u+np.mat(feedback)).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SVD++ Mean Squared Error: 0.921\n"
     ]
    }
   ],
   "source": [
    "#use mse as benchmark for comparison\n",
    "mse_svdpp=np.square((\n",
    "    output-arr)[testing_idx]).sum()/len(arr[testing_idx])\n",
    "print('SVD++ Mean Squared Error:',round(mse_svdpp,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Non-negative Matrix Factorization\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Funk SVD is cool but not without caveats. Sometimes nonnegativity is required to better interpret the causality of the latent variables. This is where Non-negative Matrix Factorization comes to the theatre. For simplicity, the official optimization problem is formulated with two latent variables only, $p_u$ and $q_i$.\n",
    "\n",
    "$$ \\min_{p_*,q_*}\\,\\sum_{r_{ui}\\,\\in\\,\\mathcal{K}} \\left(r_{ui} - q_i^Tp_u \\right)^2 + \\lambda( ||p_u||^2 + ||q_i||^2 )$$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Intuitively, we can derive the gradient descent update rule via $\\theta:=\\theta-\\alpha\\frac{\\partial J(\\theta)}{\\partial \\theta}$.\n",
    "\n",
    "$$p_u := p_u + \\alpha ( \\left( r_{ui} - \\hat{r}_{ui} \\right) \\cdot q_i - \\lambda p_u)$$\n",
    "\n",
    "$$q_i := q_i + \\alpha ( \\left( r_{ui} - \\hat{r}_{ui} \\right) \\cdot p_u - \\lambda q_i)$$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "The issue with a conventional gradient descent is that it cannot guarantee the non-negativity. Hence, we introduce a modified version called Multiplicative Update. MU can be derived from the learning rate $\\alpha$. Take user preference $p_u$ for example, if we set learning rate $\\alpha = \\frac {p_u}{\\hat{r}_{ui} \\cdot q_i + \\lambda p_u}$. We get rid of the subtraction and obtain an update rule based upon multiplication.\n",
    "\n",
    "$$p_u := p_u \\frac {r_{ui} \\cdot q_i}{\\hat{r}_{ui} \\cdot q_i + \\lambda p_u}$$\n",
    "\n",
    "Similarly, if we set learning rate $\\alpha = \\frac {q_i}{\\hat{r}_{ui} \\cdot p_u + \\lambda q_i}$ for item attribute $q_i$, we can obtain another set of MU.\n",
    "\n",
    "$$q_i := q_i \\frac {r_{ui} \\cdot p_u}{\\hat{r}_{ui} \\cdot p_u + \\lambda q_i}$$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Quite a few existing packages provide NMF. Apart from `surprise.prediction_algorithms.matrix_factorization.NMF`, you can actually find NMF in `sklearn.decomposition.NMF`.\n",
    "\n",
    "A good tutorial on NMF\n",
    "\n",
    "https://perso.telecom-paristech.fr/essid/teach/NMF_tutorial_ICME-2014.pdf\n",
    "\n",
    "Reference to the original paper of MU on NMF\n",
    "\n",
    "https://proceedings.neurips.cc/paper/2000/file/f9d1152547c0bde01830b7e8bd60024c-Paper.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use numba to dramatically boost the speed of linear algebra\n",
    "#this will be a lot faster than surprise library\n",
    "@numba.njit\n",
    "def nmf_epoch(arr,p_u,q_i,lambda_):\n",
    "    \n",
    "    #initialize\n",
    "    error=0\n",
    "    \n",
    "    #only iterate known ratings\n",
    "    for i in range(arr.shape[0]):\n",
    "        for u in range(arr.shape[1]):\n",
    "            r_ui=arr[i,u]\n",
    "            num_of_users=len(np.where(~np.isnan(arr[i]))[0])\n",
    "            num_of_items=len(np.where(~np.isnan(arr[:,u]))[0])\n",
    "            \n",
    "            #another way to identify nan\n",
    "            #r_ui!=r_ui\n",
    "            if np.isnan(r_ui):\n",
    "                continue\n",
    "\n",
    "            #compute error\n",
    "            estimated=q_i[i].T@p_u[u]\n",
    "            epsilon_ui=r_ui-estimated\n",
    "            error+=epsilon_ui**2\n",
    "\n",
    "            #update\n",
    "            p_u[u]=np.multiply(p_u[u],np.divide(q_i[i]*r_ui,\n",
    "                             q_i[i]*estimated+lambda_*p_u[u]*num_of_items))\n",
    "            q_i[i]=np.multiply(q_i[i],np.divide(p_u[u]*r_ui,\n",
    "                             p_u[u]*estimated+lambda_*q_i[i]*num_of_users))\n",
    "    \n",
    "    return error,p_u,q_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#nmf latent factor model\n",
    "def nmf(arr,p_u_init=[],q_i_init=[],num_of_rank=40,\n",
    "             lambda_=0.02,tau=0.0001,\n",
    "             max_iter=20,diagnosis=True\n",
    "             ):\n",
    "\n",
    "    #initialize\n",
    "    stop=False\n",
    "    counter=0\n",
    "    sse=None\n",
    "    \n",
    "    #user latent factors\n",
    "    if len(p_u_init)==0:\n",
    "        p_u=np.zeros((arr.shape[1],num_of_rank))\n",
    "        p_u.fill(0.1)\n",
    "    else:\n",
    "        p_u=p_u_init\n",
    "    \n",
    "    #item latent factors\n",
    "    if len(q_i_init)==0:\n",
    "        q_i=np.zeros((arr.shape[0],num_of_rank))\n",
    "        q_i.fill(0.1)\n",
    "    else:\n",
    "        q_i=q_i_init\n",
    "        \n",
    "    #multiplicative update\n",
    "    while not stop:\n",
    "        \n",
    "        error,p_u,q_i=nmf_epoch(arr,p_u,q_i,lambda_)\n",
    "\n",
    "        counter+=1\n",
    "        \n",
    "        #maximum number of epoch\n",
    "        if counter>=max_iter:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print('Not converged.',\n",
    "                      'Consider increase number of iterations or tolerance')\n",
    "                \n",
    "        #use sum of squared error to determine if converged\n",
    "        sse_prev=sse\n",
    "        sse=error\n",
    "        if sse_prev and abs(sse/sse_prev-1)<=tau:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print(f'{counter} iterations to reach convergence\\n')\n",
    "\n",
    "    return p_u,q_i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#nmf is highly sensitive to initial values\n",
    "p_u_nmf_init=np.random.uniform(size=(arr.shape[1],\n",
    "                            num_of_latent_factors))\n",
    "q_i_nmf_init=np.random.uniform(size=(arr.shape[0],\n",
    "                            num_of_latent_factors))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6 iterations to reach convergence\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#nmf\n",
    "p_u,q_i=nmf(arr_train,p_u_init=p_u_nmf_init,q_i_init=q_i_nmf_init,\n",
    "            num_of_rank=num_of_latent_factors,\n",
    "            lambda_=lagrange_multiplier,\n",
    "            tau=tolerance,\n",
    "            max_iter=max_num_of_epoch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix completion\n",
    "output=q_i@p_u.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NMF Mean Squared Error: 6.987\n"
     ]
    }
   ],
   "source": [
    "#use mse as benchmark for comparison\n",
    "mse_nmf=np.square((\n",
    "    output-arr)[testing_idx]).sum()/len(arr[testing_idx])\n",
    "print('NMF Mean Squared Error:',round(mse_nmf,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Co-clustering\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Co-clustering is an unsupervised multi-label classification in a 2D matrix. Each element in the matrix depends on two labels – row cluster and column cluster. The training strategy is simply one-versus-rest, meaning the user cluster label and item cluster label of each element are trained separately by K Means. Initially \n",
    "$K_u$ user centroids and $K_i$ item centroids are randomly assigned to each rating $r_{ui}$. Repeat the following steps until convergence. The convergence can be defined as maximum number of iterations, the fixed position of each cluster or the relative error change.\n",
    "\n",
    "* E-step: compute the distance from each data point $r_{ui}$ to $K_u$ and $K_i$ centroids. Each label of the data point is determined by the centroid with the shortest distance. \n",
    "\n",
    "* M-step: Each centroid $C_*$ is recalibrated to the mean of all the data points $r_{ui}$ which fall under its jurisdiction.\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "In the context of recommender system, the distance metric is defined as the squared error between estimated rating and actual rating. The estimated rating is influenced by the user mean $\\mu_u$ (think of it as user baseline), the item mean $\\mu_i$ (think of it as item baseline) and the co-clustering effect $\\overline{C_{ui}}$. One thing to bear in mind is that the user clustering effect $\\overline{C_{u}}$ contributes to both user mean and co-clustering effect. It has to be deducted from the estimated rating to avoid double counting. The same applies to the item clustering effect $\\overline{C_{i}}$.\n",
    "\n",
    "$$distance(r_{ui},C_*) = \\left( r_{ui}- \\hat{r}_{ui} \\right)^2 = \\left( r_{ui} - \\overline{C_{ui}} - (\\mu_u - \\overline{C_u}) - (\\mu_i - \\overline{C_i}) \\right)^2$$\n",
    "\n",
    "where \n",
    "\n",
    "$\\overline{C_{ui}}$ denotes the average rating of co-cluster $C_{ui}$\n",
    "\n",
    "$\\overline{C_u}$ denotes the average rating of user cluster $C_u$\n",
    "\n",
    "$\\overline{C_i}$ denotes the average rating of item  cluster $C_i$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Reference to K Means\n",
    "\n",
    "https://github.com/je-suis-tm/machine-learning/blob/master/k%20means.ipynb\n",
    "\n",
    "Reference to the original paper by Thomas George and Srujana Merugu\n",
    "\n",
    "https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.113.6458&rep=rep1&type=pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute mean of user,item and co-cluster\n",
    "@numba.njit\n",
    "def compute_avg(arr,num_of_cluster_item,\n",
    "           num_of_cluster_user,cluster_item,cluster_user):\n",
    "    \n",
    "    #compute co-cluster mean\n",
    "    bar_c_ui=np.zeros((num_of_cluster_item,\n",
    "                       num_of_cluster_user))\n",
    "    for c_i in range(num_of_cluster_item):\n",
    "        for c_u in range(num_of_cluster_user):\n",
    "            row_id=np.where(cluster_item==c_i)[0]\n",
    "            col_id=np.where(cluster_user==c_u)[0]\n",
    "            cocluster_sub=arr[row_id][:,col_id].copy().flatten()\n",
    "            bar_c_ui[c_i][c_u]=(cocluster_sub[~np.isnan(cocluster_sub)]).mean()\n",
    "            \n",
    "    #compute user and item mean\n",
    "    miu_i=[]\n",
    "    for i in range(arr.shape[0]):\n",
    "        subset=arr[i][~np.isnan(arr[i])]\n",
    "        if len(subset)==0:\n",
    "            miu_i.append(np.nan)\n",
    "        else:\n",
    "            miu_i.append(subset.mean())\n",
    "    miu_u=[]\n",
    "    for u in range(arr.shape[1]):\n",
    "        subset=arr[:,u][~np.isnan(arr[:,u])]\n",
    "        if len(subset)==0:\n",
    "            miu_u.append(np.nan)\n",
    "        else:\n",
    "            miu_u.append(subset.mean())\n",
    "    \n",
    "    return bar_c_ui,miu_i,miu_u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute mean of user and item cluster\n",
    "@numba.njit\n",
    "def m_step(arr,num_of_cluster_item,\n",
    "           num_of_cluster_user,cluster_item,\n",
    "           cluster_user):\n",
    "\n",
    "    #compute mean of item cluster\n",
    "    bar_cluster_item=[]\n",
    "    for c_i in range(num_of_cluster_item):\n",
    "        cluster_sub=arr[np.where(cluster_item==c_i)[0]].copy().flatten()\n",
    "        bar_cluster_item.append((cluster_sub[~np.isnan(cluster_sub)]).mean())\n",
    "\n",
    "    #compute mean of user cluster\n",
    "    bar_cluster_user=[]\n",
    "    for c_u in range(num_of_cluster_user):\n",
    "        cluster_sub=arr[:,np.where(cluster_user==c_u)[0]].copy().flatten()\n",
    "        bar_cluster_user.append((cluster_sub[~np.isnan(cluster_sub)]).mean())\n",
    "    \n",
    "    #compute mean of user,item and co-cluster\n",
    "    bar_c_ui,miu_i,miu_u=compute_avg(arr,num_of_cluster_item,\n",
    "                       num_of_cluster_user,cluster_item,cluster_user)\n",
    "    \n",
    "    return bar_c_ui,miu_i,miu_u,bar_cluster_item,bar_cluster_user"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "#reassign user cluster\n",
    "@numba.njit\n",
    "def get_cluster_user(arr,bar_c_ui,miu_i,miu_u,num_of_cluster_item,\n",
    "           num_of_cluster_user,cluster_item,cluster_user,\n",
    "           bar_cluster_item,bar_cluster_user):\n",
    "    \n",
    "    #iterate through users\n",
    "    for u in range(arr.shape[1]):\n",
    "          \n",
    "        #get user cluster\n",
    "        sorting=[]\n",
    "        for c_u in range(num_of_cluster_user):\n",
    "            sorting.append(0)\n",
    "            for i in np.where(~np.isnan(arr[:,u]))[0]:            \n",
    "                c_i=cluster_item[i]\n",
    "                r_ui=arr[i,u]\n",
    "    \n",
    "                #compute estimated\n",
    "                if np.isnan(miu_u[u]):\n",
    "                    user_part=0\n",
    "                else:\n",
    "                    user_part=miu_u[u]-bar_cluster_user[c_u]\n",
    "                if np.isnan(miu_i[i]):\n",
    "                    item_part=0\n",
    "                else:\n",
    "                    item_part=miu_i[i]-bar_cluster_item[c_i]\n",
    "                if np.isnan(bar_c_ui[c_i,c_u]):\n",
    "                    co_part=arr.flatten()[~np.isnan(arr.flatten())].mean()\n",
    "                else:\n",
    "                    co_part=bar_c_ui[c_i,c_u]\n",
    "                est=co_part+user_part+item_part\n",
    "    \n",
    "                #compute sum of squared error\n",
    "                sorting[-1]+=(r_ui-est)**2\n",
    "        \n",
    "        #assign to the user cluster that minimizes sse\n",
    "        cluster_user[u]=sorting.index(list(set(sorting))[0])\n",
    "    \n",
    "    return cluster_user"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "#reassign item cluster\n",
    "@numba.njit\n",
    "def get_cluster_item(arr,bar_c_ui,miu_i,miu_u,num_of_cluster_item,\n",
    "           num_of_cluster_user,cluster_item,cluster_user,\n",
    "           bar_cluster_item,bar_cluster_user):\n",
    "    \n",
    "    #iterate through items\n",
    "    for i in range(arr.shape[0]):\n",
    "            \n",
    "        #get item cluster\n",
    "        sorting=[]\n",
    "        for c_i in range(num_of_cluster_item):\n",
    "            sorting.append(0)\n",
    "            for u in np.where(~np.isnan(arr[i]))[0]: \n",
    "                c_u=cluster_user[u]\n",
    "                r_ui=arr[i,u]\n",
    "    \n",
    "                #compute estimated\n",
    "                if np.isnan(miu_u[u]):\n",
    "                    user_part=0\n",
    "                else:\n",
    "                    user_part=miu_u[u]-bar_cluster_user[c_u]\n",
    "                if np.isnan(miu_i[i]):\n",
    "                    item_part=0\n",
    "                else:\n",
    "                    item_part=miu_i[i]-bar_cluster_item[c_i]\n",
    "                if np.isnan(bar_c_ui[c_i,c_u]):\n",
    "                    co_part=arr.flatten()[~np.isnan(arr.flatten())].mean()\n",
    "                else:\n",
    "                    co_part=bar_c_ui[c_i,c_u]\n",
    "                est=co_part+user_part+item_part\n",
    "                \n",
    "                #compute sum of squared error\n",
    "                sorting[-1]+=(r_ui-est)**2\n",
    "        \n",
    "        #assign to the item cluster that minimizes sse\n",
    "        cluster_item[i]=sorting.index(list(set(sorting))[0])\n",
    "    \n",
    "    return cluster_item"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "#reassign user and item clusters\n",
    "@numba.njit\n",
    "def e_step(arr,bar_c_ui,miu_i,miu_u,\n",
    "           num_of_cluster_item,num_of_cluster_user,\n",
    "           cluster_item,cluster_user,\n",
    "           bar_cluster_item,bar_cluster_user):\n",
    "\n",
    "    cluster_user=get_cluster_user(arr,bar_c_ui,miu_i,miu_u,\n",
    "                                  num_of_cluster_item,num_of_cluster_user,\n",
    "                                  cluster_item,cluster_user,\n",
    "                                   bar_cluster_item,bar_cluster_user)\n",
    "\n",
    "    cluster_item=get_cluster_item(arr,bar_c_ui,miu_i,miu_u,\n",
    "                                  num_of_cluster_item,num_of_cluster_user,\n",
    "                                  cluster_item,cluster_user,\n",
    "                                   bar_cluster_item,bar_cluster_user)\n",
    "    \n",
    "    return cluster_item,cluster_user"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix completion\n",
    "@numba.njit\n",
    "def cocluster_predict(arr,bar_c_ui,miu_i,miu_u,bar_cluster_item,\n",
    "                      bar_cluster_user,cluster_item,cluster_user):\n",
    "\n",
    "    #map cocluster mean into array shape\n",
    "    cluster_mapping=[bar_c_ui[\n",
    "        cluster_item[i],cluster_user[u]] for i in range(\n",
    "        arr.shape[0]) for u in range(arr.shape[1])]\n",
    "    cocluster_mean=np.array(cluster_mapping)\n",
    "\n",
    "    #fill na with global mean\n",
    "    cocluster_mean[np.isnan(cocluster_mean)]=arr.flatten()[\n",
    "        ~np.isnan(arr.flatten())].mean()\n",
    "    copart=cocluster_mean.reshape(arr.shape)\n",
    "    \n",
    "    #map item mean into array shape\n",
    "    item_cluster_mean=np.array(\n",
    "        [bar_cluster_item[c_i] for c_i in cluster_item for _ in range(\n",
    "            arr.shape[1])])    \n",
    "    item_mean=np.array(\n",
    "        [bar_i for bar_i in miu_i for _ in range(\n",
    "            arr.shape[1])])\n",
    "\n",
    "    #compute\n",
    "    item_part=item_mean-item_cluster_mean\n",
    "\n",
    "    #fill na with zero\n",
    "    item_part[np.isnan(item_part)]=0\n",
    "    item_part=item_part.reshape(arr.shape)\n",
    "\n",
    "    #map user mean into array shape\n",
    "    user_cluster_mean=np.array([bar_cluster_user[\n",
    "                c_u] for c_u in cluster_user]*arr.shape[0])\n",
    "    user_mean=np.array(miu_u*arr.shape[0])           \n",
    "\n",
    "    #compute\n",
    "    user_part=user_mean-user_cluster_mean\n",
    "\n",
    "    #fill na with zero\n",
    "    user_part[np.isnan(user_part)]=0\n",
    "    user_part=user_part.reshape(arr.shape)\n",
    "    \n",
    "    #forecast\n",
    "    est=copart+user_part+item_part\n",
    "    \n",
    "    return est"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "#co-clustering model\n",
    "def coclustering(arr,cluster_user_init=[],cluster_item_init=[],\n",
    "                 num_of_cluster_user=3,\n",
    "                 num_of_cluster_item=3,tau=0.0001,\n",
    "                 max_iter=200,diagnosis=True\n",
    "                 ):\n",
    "\n",
    "    #initialize\n",
    "    stop=False\n",
    "    counter=0\n",
    "    sse=None\n",
    "    \n",
    "    #initialize item cluster\n",
    "    if len(cluster_item_init)==0:\n",
    "        cluster_item=np.random.choice(\n",
    "            range(num_of_cluster_item),arr.shape[0])\n",
    "    else:\n",
    "        cluster_item=cluster_item_init\n",
    "    \n",
    "    #initialize user cluster\n",
    "    if len(cluster_user_init)==0:\n",
    "        cluster_user=np.random.choice(\n",
    "            range(num_of_cluster_user),arr.shape[1])\n",
    "    else:\n",
    "        cluster_user=cluster_user_init\n",
    "            \n",
    "    #compute mean of user and item cluster\n",
    "    bar_c_ui,miu_i,miu_u,bar_cluster_item,bar_cluster_user=m_step(\n",
    "            arr,num_of_cluster_item,num_of_cluster_user,cluster_item,\n",
    "                   cluster_user)\n",
    "    \n",
    "    #keep track of cluster assignment\n",
    "    cluster_item_prev=cluster_item.copy()\n",
    "    cluster_user_prev=cluster_user.copy()\n",
    "    \n",
    "    #k means\n",
    "    while not stop:\n",
    "        \n",
    "        #e step\n",
    "        cluster_item,cluster_user=e_step(arr,bar_c_ui,miu_i,miu_u,\n",
    "                           num_of_cluster_item,num_of_cluster_user,\n",
    "                           cluster_item,cluster_user,\n",
    "                           bar_cluster_item,bar_cluster_user)\n",
    "\n",
    "        #m step\n",
    "        bar_c_ui,miu_i,miu_u,bar_cluster_item,bar_cluster_user=m_step(\n",
    "            arr,num_of_cluster_item,num_of_cluster_user,cluster_item,\n",
    "                   cluster_user)\n",
    "\n",
    "        counter+=1\n",
    "        \n",
    "        #maximum number of epoch\n",
    "        if counter>=max_iter:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print('Not converged.',\n",
    "                      'Consider increase number of iterations',\n",
    "                      'or change number of clusters')\n",
    "                \n",
    "        #check cluster assignment to determine if converged\n",
    "        if (cluster_item_prev==cluster_item).all() and \\\n",
    "        (cluster_user_prev==cluster_user).all():\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print(f'{counter} iterations to reach convergence\\n')\n",
    "        \n",
    "        #use sum of squared error to determine if converged\n",
    "        estimated=cocluster_predict(arr,bar_c_ui,\n",
    "                                 miu_i,miu_u,bar_cluster_item,\n",
    "                                  bar_cluster_user,\n",
    "                                 cluster_item,cluster_user)\n",
    "        sse_prev=sse\n",
    "        sse=np.square(arr-estimated).sum()\n",
    "        if sse_prev and abs(sse/sse_prev-1)<=tau:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print(f'{counter} iterations to reach convergence\\n')\n",
    "\n",
    "    return bar_c_ui,miu_i,miu_u,bar_cluster_item, \\\n",
    "bar_cluster_user,cluster_item,cluster_user"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize the number of clusters\n",
    "num_c_u=3\n",
    "num_c_i=3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Not converged. Consider increase number of iterations or change number of clusters\n"
     ]
    }
   ],
   "source": [
    "#co-clustering\n",
    "bar_c_ui,miu_i,miu_u,bar_cluster_item,bar_cluster_user, \\\n",
    "cluster_item,cluster_user=coclustering(\n",
    "             arr_train,num_of_cluster_user=num_c_u,\n",
    "             num_of_cluster_item=num_c_i,\n",
    "             tau=tolerance,\n",
    "             max_iter=max_num_of_epoch\n",
    "                 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix completion\n",
    "output=cocluster_predict(arr_train,bar_c_ui,miu_i,miu_u,bar_cluster_item,\n",
    "                      bar_cluster_user,cluster_item,cluster_user)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Co-clustering Mean Squared Error: 0.983\n"
     ]
    }
   ],
   "source": [
    "#use mse as benchmark for comparison\n",
    "mse_cocluster=np.square((\n",
    "    output-arr)[testing_idx]).sum()/len(arr[testing_idx])\n",
    "print('Co-clustering Mean Squared Error:',round(mse_cocluster,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Memory-based"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Slope One\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Slope one is probably the easiest memory-based method. The name comes from the naïve form of $f(x) = kx + b$. In this particular version, slope $k$ is set at $1$, hence the name `Slope One`. The intercept $b$ is defined as the user mean $\\mu_u$ and the variable $x$ is defined as the average deviation of relevant items. By using the previous notations, we are able to obtain the following equation of slope one.\n",
    "\n",
    "$$\\hat{r}_{ui} = \\frac{1}{|R_i(u)|}\n",
    "        \\sum\\limits_{j \\in R_i(u)} \\text{dev}(i, j) + \\mu_u$$\n",
    "\n",
    "where\n",
    "\n",
    "$\\mu_u$ denotes the average rating by user $u$\n",
    "\n",
    "$R_i(u)$ denotes the set of relevant items which both user $u$ and user $v$, who has rated item $i$, have rated\n",
    "\n",
    "$\\text{dev}_(i, j)$ denotes the average rating deviation between item $i$ and relevant item $j$ rated by users who have rated item $i$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "By the definition of $\\text{dev}_(i, j)$, the average deviation of relevant items can be thought as the relative item baseline $b_i$ which explains why item $i$ attracts higher/lower score than item $j$ across all users. The equation is formulated below.\n",
    "\n",
    "$$\\text{dev}(i, j) = \\frac{1}{        |U_{ij}|}\\sum\\limits_{u \\in U_{ij}} r_{ui} - r_{uj}$$\n",
    "\n",
    "where\n",
    "\n",
    "$U_{ij}$ denotes the set of users who have rated both item $i$ and relevant item $j$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Reference to the original paper by Daniel Lemire and Anna Maclachlan\n",
    "\n",
    "https://www.researchgate.net/publication/1960789_Slope_One_Predictors_for_Online_Rating-Based_Collaborative_Filtering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "#basic slope one algorithm\n",
    "@numba.njit\n",
    "def slope_one(arr,miu):\n",
    "                \n",
    "    #create a copy\n",
    "    estimated=arr.copy()\n",
    "        \n",
    "    for i in range(arr.shape[0]):\n",
    "        for u in range(arr.shape[1]):\n",
    "            if np.isnan(arr[i,u]):\n",
    "                \n",
    "                #get user mean\n",
    "                miu_u=(arr[:,u][~np.isnan(arr[:,u])]).mean()\n",
    "                \n",
    "                #take global mean for cold start\n",
    "                if np.isnan(miu_u):\n",
    "                    miu_u=miu\n",
    "                \n",
    "                #items rated by user u\n",
    "                items_rated_u=set(np.where(~np.isnan(arr[:,u]))[0])\n",
    "\n",
    "                #users who have rated item i\n",
    "                users_rated_i=np.where(~np.isnan(arr[i]))[0]\n",
    "                \n",
    "                #items rated by users who have rated item i\n",
    "                items_rated_others=set(np.where(~np.isnan(\n",
    "                    arr_train[:,users_rated_i]))[0])\n",
    "                    \n",
    "                #find relevant items\n",
    "                R_i_u=items_rated_u.intersection(items_rated_others)\n",
    "\n",
    "                #take user mean if user u share nothing in common with others\n",
    "                if len(R_i_u)==0:                    \n",
    "                    estimated[i,u]=miu_u\n",
    "                else:\n",
    "\n",
    "                    #compute average deviation of relevant items\n",
    "                    relevant_items=[]\n",
    "                    for j in R_i_u:\n",
    "                                                   \n",
    "                        #find users who have rated both i and j\n",
    "                        U_i_j=set(users_rated_i).intersection(\n",
    "                            set(np.where(~np.isnan(arr[j]))[0]))\n",
    "\n",
    "                        #compute average deviation \n",
    "                        #for the users who have rated both i and j\n",
    "                        dev_i_j=np.mean(arr[i][\n",
    "                            np.array(list(U_i_j))]-arr[j][\n",
    "                            np.array(list(U_i_j))])\n",
    "                        relevant_items.append(dev_i_j)\n",
    "                    \n",
    "                    #predict\n",
    "                    estimated[i,u]=miu_u+np.mean(np.array(relevant_items))\n",
    "                    \n",
    "    return estimated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix completion\n",
    "output=slope_one(arr_train,miu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Slope One Mean Squared Error: 0.985\n"
     ]
    }
   ],
   "source": [
    "#use mse as benchmark for comparison\n",
    "mse_slop_one=np.square((\n",
    "    output-arr)[testing_idx]).sum()/len(arr[testing_idx])\n",
    "print('Slope One Mean Squared Error:',round(mse_slop_one,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### K Nearest Neighbors\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "KNN is undoubtedly an easy supervised learning model. In the iris dataset, we have shown how KNN works in classification problem. In the case of recommender system, it is effectively a regression problem. We have to slightly change the methodology. \n",
    "\n",
    "* KNN in recommender system takes the weighted average of neighbors to predict the rating, albeit KNN in the iris dataset takes the majority label of neighbors to predict the classification.\n",
    "\n",
    "* KNN in recommender system uses cosine similarity, mean squared distance similarity or Pearson correlation similarity to define neighbors, whereas KNN in the iris dataset uses Euclidean distance, Chebyshev distance or Manhattan distance to define neighbors.\n",
    "\n",
    "In KNN, you also have the choice of picking neighbors based upon the users (how similar users rate the same item) or the items (how similar items get rated by the same user). In this script, every neighborhood model will be based upon the users. For item approach, just switch user and item in the equation below.\n",
    "\n",
    "$$\\hat{r}_{ui} = \\frac{  \\sum_{v \\in N^k_i(u)} \\mathcal{similarity}(u, v) \\cdot r_{vi}}  {\\sum_{v \\in N^k_i(u)} \\mathcal{similarity}(u, v)}$$\n",
    "\n",
    "where \n",
    "\n",
    "$v$ denotes the neighbors of user $u$\n",
    "\n",
    "$N^k_i(u)$ denotes the set of top $k$ neighbors of user $u$ that have rated item $i$\n",
    "\n",
    "$\\mathcal{similarity}(u, v)$ denotes the similarity metric between user $u$ and its neighbor $v$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Each similarity metric has its own pros and cons. In this case, we will pick the MSD similarity which is formulated below. The biggest selling point of MSD is its resemblance to Euclidean distance. One of its malaise is it doesn't penalize a small number of $|I_{uv}|$. A sum of squared error of 0 with $|I_{uv}|$ at 1 is inevitably inferior to a sum of squared error of 1 with $|I_{uv}|$ at 10 but the latter yields a smaller MSD even the latter pair has rated more items in common.\n",
    "\n",
    "$$\\mathcal{similarity}(u, v) = \\frac{1}{\\frac{\\sum_{i \\in I_{uv}} (r_{ui} - r_{vi})^2}{|I_{uv}|}\n",
    "         + 1}$$\n",
    "\n",
    "where \n",
    "\n",
    "$I_{uv}$ denotes the set of items where both user $u$ and its neighbor $v$ have rated\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Reference to KNN in iris dataset\n",
    "\n",
    "https://github.com/je-suis-tm/machine-learning/blob/master/k%20nearest%20neighbors.ipynb\n",
    "\n",
    "Reference to cosine similarity\n",
    "\n",
    "https://github.com/je-suis-tm/machine-learning/blob/master/latent%20semantic%20indexing.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "#obtain msd similarity matrix\n",
    "@numba.njit\n",
    "def get_msd_similarity_matrix(arr):\n",
    "\n",
    "    similarity_matrix=np.zeros((arr.shape[1],arr.shape[1]))\n",
    "    for u in range(arr.shape[1]):\n",
    "        for v in range(u+1,arr.shape[1]):\n",
    "            \n",
    "            #self correlation distorts knn selection\n",
    "            if u==v:\n",
    "                continue\n",
    "\n",
    "            #compute msd first then eliminate nan\n",
    "            I_uv=np.square(arr[:,u]-arr[:,v])\n",
    "            valid_I_uv=I_uv[~np.isnan(I_uv)]\n",
    "\n",
    "            #avoid the case where two users havent rated any items in common\n",
    "            if len(valid_I_uv)>0:\n",
    "                msd=1/(valid_I_uv.sum()/len(valid_I_uv)+1)\n",
    "            else:\n",
    "                msd=0\n",
    "\n",
    "            #symmetric matrix\n",
    "            similarity_matrix[u,v]=msd\n",
    "            similarity_matrix[v,u]=msd\n",
    "            \n",
    "    return similarity_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "#knn for matrix completion\n",
    "@numba.njit\n",
    "def knn(arr,similarity_matrix,top_k=10):\n",
    "\n",
    "    estimated=arr.copy()\n",
    "    for i in range(arr.shape[0]):\n",
    "        for u in range(arr.shape[1]):\n",
    "            if np.isnan(arr[i,u]):\n",
    "\n",
    "                #find top k neighbor based upon similarity matrix\n",
    "                rated_users=np.where(~np.isnan(arr[i]))[0]\n",
    "                similarities=similarity_matrix[u][rated_users]\n",
    "                top_k_neighbors=np.argsort(similarities)[-top_k:]\n",
    "                N_k_i_u=np.array([\n",
    "                    rated_users[neighbor] for neighbor in top_k_neighbors])\n",
    "\n",
    "                #compute weighted average                \n",
    "                if len(N_k_i_u)!=0:\n",
    "                    numerator=arr[i][N_k_i_u].T@similarity_matrix[u][N_k_i_u]\n",
    "                    denominator=similarity_matrix[u][N_k_i_u].sum()\n",
    "                    if denominator!=0:\n",
    "                        estimated[i,u]=numerator/denominator\n",
    "                        continue\n",
    "                        \n",
    "                #when the users who rated item i\n",
    "                #have nothing in common with user u\n",
    "                #take the average of user mean and item mean\n",
    "                #if both user mean and item mean are empty\n",
    "                #take global mean\n",
    "                miu_i=arr[i][~np.isnan(arr[i])]\n",
    "                miu_u=arr[:,u][~np.isnan(arr[:,u])]                \n",
    "                if len(miu_i)==0 and len(miu_u)==0:\n",
    "                    estimated[i,u]=arr[np.where(~np.isnan(arr))[0]][\n",
    "                        np.where(~np.isnan(arr))[1]].mean()\n",
    "                else:\n",
    "                    estimated[i,u]=np.array(list(miu_i)+list(miu_u)).mean()\n",
    "                    \n",
    "    return estimated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_k=10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute similarity matrix\n",
    "similarity_matrix=get_msd_similarity_matrix(arr_train)\n",
    "\n",
    "#matrix completion\n",
    "output=knn(arr_train,similarity_matrix,top_k=top_k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KNN Mean Squared Error: 0.986\n"
     ]
    }
   ],
   "source": [
    "#use mse as benchmark for comparison\n",
    "mse_knn=np.square((\n",
    "    output-arr)[testing_idx]).sum()/len(arr[testing_idx])\n",
    "print('KNN Mean Squared Error:',round(mse_knn,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### K Nearest Neighbors with mean\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Now that we are familiar with KNN, it's about time to introduce many of its variations. The first one is KNN with mean. In the user approach, it takes consideration of user mean and neighbor mean into the equation.\n",
    "\n",
    "$$\\hat{r}_{ui} = \\mu_u + \\frac{  \\sum_{v \\in N^k_i(u)} \\mathcal{similarity}(u, v) \\cdot (r_{vi}-\\mu_v)}  {\\sum_{v \\in N^k_i(u)} \\mathcal{similarity}(u, v)}$$\n",
    "\n",
    "where \n",
    "\n",
    "$\\mu_u$ denotes the average rating of user $u$\n",
    "\n",
    "$\\mu_v$ denotes the average rating of neighbor $v$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "#knn with mean for matrix completion\n",
    "@numba.njit\n",
    "def knn_with_mean(arr,similarity_matrix,top_k=10):\n",
    "    \n",
    "    #compute user mean\n",
    "    user_mean=np.array([arr[:,u][\n",
    "        ~np.isnan(arr[:,u])].mean() for u in range(arr.shape[1])])\n",
    "\n",
    "    estimated=arr.copy()\n",
    "    for i in range(arr.shape[0]):\n",
    "        for u in range(arr.shape[1]):\n",
    "            if np.isnan(arr[i,u]):\n",
    "\n",
    "                #find top k neighbor based upon similarity matrix\n",
    "                rated_users=np.where(~np.isnan(arr[i]))[0]\n",
    "                similarities=similarity_matrix[u][rated_users]\n",
    "                top_k_neighbors=np.argsort(similarities)[-top_k:]\n",
    "                N_k_i_u=np.array([\n",
    "                    rated_users[neighbor] for neighbor in top_k_neighbors])\n",
    "\n",
    "                #compute weighted average                \n",
    "                if len(N_k_i_u)!=0:\n",
    "                    numerator=(arr[i][\n",
    "                        N_k_i_u]-user_mean[\n",
    "                        N_k_i_u]).T@similarity_matrix[u][N_k_i_u]\n",
    "                    denominator=similarity_matrix[u][N_k_i_u].sum()\n",
    "                    if denominator!=0:\n",
    "                        estimated[i,u]=user_mean[u]+numerator/denominator\n",
    "                        continue\n",
    "                        \n",
    "                #when the users who rated item i\n",
    "                #have nothing in common with user u\n",
    "                #take the average of user mean and item mean\n",
    "                #if both user mean and item mean are empty\n",
    "                #take global mean\n",
    "                miu_i=arr[i][~np.isnan(arr[i])]\n",
    "                miu_u=arr[:,u][~np.isnan(arr[:,u])]                \n",
    "                if len(miu_i)==0 and len(miu_u)==0:\n",
    "                    estimated[i,u]=arr[np.where(~np.isnan(arr))[0]][\n",
    "                        np.where(~np.isnan(arr))[1]].mean()\n",
    "                else:\n",
    "                    estimated[i,u]=np.array(list(miu_i)+list(miu_u)).mean()\n",
    "                    \n",
    "    return estimated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix completion\n",
    "output=knn_with_mean(arr_train,similarity_matrix,top_k=top_k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KNN with μ Mean Squared Error: 0.971\n"
     ]
    }
   ],
   "source": [
    "#use mse as benchmark for comparison\n",
    "mse_knn_with_miu=np.square((\n",
    "    output-arr)[testing_idx]).sum()/len(arr[testing_idx])\n",
    "print('KNN with μ Mean Squared Error:',round(mse_knn_with_miu,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### K Nearest Neighbors with z score\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "On top of KNN with mean, we can include the standard deviation of user into the equation. After all, some users may be moody so their ratings tend to have a larger fluctuation.\n",
    "\n",
    "$$\\hat{r}_{ui} = \\mu_u + \\sigma_u \\cdot \\frac{  \\sum_{v \\in N^k_i(u)} \\mathcal{similarity}(u, v) \\cdot \\frac {(r_{vi}-\\mu_v)} {\\sigma_v}  }  {\\sum_{v \\in N^k_i(u)} \\mathcal{similarity}(u, v) } $$\n",
    "\n",
    "where \n",
    "\n",
    "$\\sigma_u$ denotes the standard deviation of the rating from user $u$\n",
    "\n",
    "$\\sigma_v$ denotes the standard deviation of the rating from neighbor $v$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "#knn with z score for matrix completion\n",
    "@numba.njit\n",
    "def knn_with_zscore(arr,similarity_matrix,top_k=10):\n",
    "    \n",
    "    #compute user mean and std\n",
    "    user_mean=np.array([arr[:,u][\n",
    "        ~np.isnan(arr[:,u])].mean() for u in range(arr.shape[1])])\n",
    "    user_std=np.array([arr[:,u][\n",
    "        ~np.isnan(arr[:,u])].std() for u in range(arr.shape[1])])\n",
    "\n",
    "    estimated=arr.copy()\n",
    "    for i in range(arr.shape[0]):\n",
    "        for u in range(arr.shape[1]):\n",
    "            if np.isnan(arr[i,u]):\n",
    "\n",
    "                #find top k neighbor based upon similarity matrix\n",
    "                rated_users=np.where(~np.isnan(arr[i]))[0]\n",
    "                similarities=similarity_matrix[u][rated_users]\n",
    "                top_k_neighbors=np.argsort(similarities)[-top_k:]\n",
    "                N_k_i_u=np.array([\n",
    "                    rated_users[neighbor] for neighbor in top_k_neighbors])\n",
    "\n",
    "                #compute weighted average                \n",
    "                if len(N_k_i_u)!=0:\n",
    "                    z_score=np.divide((arr[i][N_k_i_u]-user_mean[N_k_i_u]),\n",
    "                                      user_std[N_k_i_u])\n",
    "                    numerator=z_score.T@similarity_matrix[u][N_k_i_u]\n",
    "                    denominator=similarity_matrix[u][N_k_i_u].sum()\n",
    "                    if denominator!=0:\n",
    "                        estimated[i,u]=user_mean[u]+user_std[\n",
    "                            u]*numerator/denominator\n",
    "                        continue\n",
    "                        \n",
    "                #when the users who rated item i\n",
    "                #have nothing in common with user u\n",
    "                #take the average of user mean and item mean\n",
    "                #if both user mean and item mean are empty\n",
    "                #take global mean\n",
    "                miu_i=arr[i][~np.isnan(arr[i])]\n",
    "                miu_u=arr[:,u][~np.isnan(arr[:,u])]                \n",
    "                if len(miu_i)==0 and len(miu_u)==0:\n",
    "                    estimated[i,u]=arr[np.where(~np.isnan(arr))[0]][\n",
    "                        np.where(~np.isnan(arr))[1]].mean()\n",
    "                else:\n",
    "                    estimated[i,u]=np.array(list(miu_i)+list(miu_u)).mean()\n",
    "                    \n",
    "    return estimated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix completion\n",
    "output=knn_with_zscore(arr_train,similarity_matrix,top_k=top_k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KNN with Z Score Mean Squared Error: 0.997\n"
     ]
    }
   ],
   "source": [
    "#use mse as benchmark for comparison\n",
    "mse_knn_with_zscore=np.square((\n",
    "    output-arr)[testing_idx]).sum()/len(arr[testing_idx])\n",
    "print('KNN with Z Score Mean Squared Error:',round(mse_knn_with_zscore,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hybird\n",
    "\n",
    "#### K Nearest Neighbors with baseline\n",
    "\n",
    "##### Baseline Model via ALS\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Before we talk about KNN with baseline, it is necessarily to talk about the baseline model. Even though we have given extensive coverage of baseline in Funk SVD, KNN does not really do any gradient descent. Thus, the baseline involved in KNN has to be estimated separately. The baseline model is formulated below. It is merely Funk SVD without latent factors.\n",
    "\n",
    "$$ \\min_{b_*}\\,\\sum_{r_{ui}\\,\\in\\,\\mathcal{K}} \\left(r_{ui} - \\mu - b_u - b_i \\right)^2 + \\lambda( b_u^2 + b_i^2 )$$\n",
    "\n",
    "In gradient descent, the equation can be solved via\n",
    "\n",
    "$$b_u := b_u + \\alpha (\\epsilon_{ui} - \\lambda b_u)$$\n",
    "$$b_i := b_i + \\alpha (\\epsilon_{ui} - \\lambda b_i)$$\n",
    "\n",
    "However, we intend to introduce another method called Alternating Least Square. The way it works is similar to coordinate descent. At each epoch, we estimate one parameter and fix the rest. The haute couture of ALS is its ability to estimate each parameter in parallel. It can execute a lot quicker than GD in high dimensional data. Another benefit is we do not need to tune the learning rate of GD. The derivation of ALS is formulated below.\n",
    "\n",
    "Let's start to estimate $b_u$. The term $r_{ui}-\\mu-b_i$ can be considered as a constant $\\phi$. Yet, in a standard least square problem, the equation is in the form of $y=x \\cdot \\beta^T$. We have to assume $b_u$ follows the form of $b_u \\cdot I^T$.\n",
    "\n",
    "$$ J(b_u)= \\left(\\phi - b_u \\cdot I^T \\right)^{T} \\left(\\phi - b_u \\cdot I^T \\right) + \\lambda( b_u^2 + b_i^2 )$$\n",
    "\n",
    "$$\\frac{\\partial J(b_u)}{\\partial\\,b_u}=2 I^T \\left(\\phi - b_u \\cdot I^T \\right)-2\\lambda b_u=0$$\n",
    "\n",
    "By using $I \\cdot b_u^T=b_u \\cdot I^T$, we obtain\n",
    "\n",
    "$$\\frac{\\partial J(b_u)}{\\partial\\,b_u}=I^T \\phi -I^T I \\cdot b_u^T - \\lambda b_u=0$$\n",
    "\n",
    "$$ b_u := \\left( I^T I + \\lambda \\right) ^ {-1} I^T \\phi $$\n",
    "\n",
    "Now $b_u$ can be estimated from a Ridge Regression estimator. In our case, $I$ is a vector of ones. $b_u \\cdot I^T$ is an operation which repeats $b_u$ into the dimension of a full matrix $r$ (or $\\phi$). $I^T I$ is normally the length of vector $b_u$ (a scalar) but $r$ is an incomplete matrix so $I^T I$ denotes a vector that contains the number of users which rated item $i$. $I^T \\phi$ is merely the summation of $r_{ui}-\\mu-b_i$. Similarly, we define $\\psi=r_{ui}-\\mu-b_u$ and $H$ is a vector of ones. We can easily derive the update rule for $b_i$ as well.\n",
    "\n",
    "$$ b_i := \\left( H^T H + \\lambda \\right) ^ {-1} H^T \\psi $$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Reference to ALS\n",
    "\n",
    "https://solomonik.cs.illinois.edu/teaching/cs554/slides/slides_15.pdf\n",
    "\n",
    "Reference to Coordinate Descent\n",
    "\n",
    "https://github.com/je-suis-tm/machine-learning/blob/master/coordinate%20descent%20for%20elastic%20net.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use numba to dramatically boost the speed of linear algebra\n",
    "#this will be a lot faster than surprise library\n",
    "@numba.njit\n",
    "def als_epoch(arr,miu,b_u,b_i,lambda_):\n",
    "    \n",
    "    #initialize\n",
    "    error=0\n",
    "    \n",
    "    #avoid np repeat to use numba\n",
    "    #equivalent to np.repeat(b_u.reshape(1,-1),arr.shape[0],axis=0)\n",
    "    #np.repeat(b_i.reshape(-1,1),arr.shape[1],axis=1)\n",
    "    b_u_map=np.array(list(b_u)*arr.shape[0]).reshape(arr.shape)\n",
    "    b_i_map=np.array(list(b_i)*arr.shape[1]).reshape(\n",
    "        (arr.shape[1],arr.shape[0])).T\n",
    "    \n",
    "    #estimate b_i while fix b_u\n",
    "    phi=arr-miu-b_u_map\n",
    "    H_HT=np.array(\n",
    "        [arr[i][~np.isnan(arr[i])].shape[0] for i in range(arr.shape[0])])\n",
    "    numerator=[phi[i][~np.isnan(phi[i])].sum() for i in range(phi.shape[0])]\n",
    "    b_i=np.divide(np.array(numerator),(lambda_+H_HT))\n",
    "    \n",
    "    #estimate b_u while fix b_i\n",
    "    psi=arr-miu-b_i_map\n",
    "    I_IT=np.array(\n",
    "        [arr[:,u][~np.isnan(arr[:,u])].shape[0] for u in range(arr.shape[1])])\n",
    "    numerator=[psi[:,u][~np.isnan(psi[:,u])].sum() for u in range(psi.shape[1])]\n",
    "    b_u=np.divide(np.array(numerator),(lambda_+I_IT))\n",
    "    \n",
    "    #update repeat\n",
    "    b_u_map=np.array(list(b_u)*arr.shape[0]).reshape(arr.shape)\n",
    "    b_i_map=np.array(list(b_i)*arr.shape[1]).reshape(\n",
    "        (arr.shape[1],arr.shape[0])).T\n",
    "    \n",
    "    #compute error\n",
    "    epsilon=(arr-miu-b_u_map-b_i_map).ravel()\n",
    "    error=np.square(epsilon[~np.isnan(epsilon)]).sum()\n",
    "    \n",
    "    return error,b_u,b_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "#als to estimate baseline model\n",
    "def baseline_als(arr,miu_init=None,b_u_init=[],b_i_init=[],\n",
    "             lambda_=0.02,tau=0.0001,\n",
    "             max_iter=20,diagnosis=True\n",
    "             ):\n",
    "\n",
    "    #initialize\n",
    "    stop=False\n",
    "    counter=0\n",
    "    sse=None\n",
    "    \n",
    "    #global mean\n",
    "    if not miu_init:       \n",
    "        miu=arr[~np.isnan(arr)].mean()\n",
    "    else:\n",
    "        miu=miu_init\n",
    "        \n",
    "    #user baseline\n",
    "    if len(b_u_init)==0:\n",
    "        b_u=np.zeros(arr.shape[1])\n",
    "    else:\n",
    "        b_u=b_u_init\n",
    "    \n",
    "    #item baseline\n",
    "    if len(b_i_init)==0:\n",
    "        b_i=np.zeros(arr.shape[0])\n",
    "    else:\n",
    "        b_i=b_i_init\n",
    "       \n",
    "    #als\n",
    "    while not stop:\n",
    "        \n",
    "        error,b_u,b_i=als_epoch(arr,miu,b_u,b_i,\n",
    "                                lambda_)\n",
    "\n",
    "        counter+=1\n",
    "\n",
    "        #maximum number of epoch\n",
    "        if counter>=max_iter:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print('Not converged.',\n",
    "                      'Consider increase number of iterations or tolerance')\n",
    "                \n",
    "        #use sum of squared error to determine if converged\n",
    "        sse_prev=sse\n",
    "        sse=error\n",
    "        \n",
    "        if sse_prev and abs(sse/sse_prev-1)<=tau:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print(f'{counter} iterations to reach convergence\\n')\n",
    "\n",
    "    return b_u,b_i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 iterations to reach convergence\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#als to estimate baseline\n",
    "b_u,b_i=baseline_als(arr_train,\n",
    "         lambda_=lagrange_multiplier,tau=tolerance,\n",
    "         max_iter=max_num_of_epoch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute global mean\n",
    "miu=arr_train[~np.isnan(arr_train)].mean()\n",
    "\n",
    "#matrix completion\n",
    "output=miu+np.repeat(\n",
    "            b_u.reshape(1,-1),\n",
    "            arr_train.shape[0],axis=0)+np.repeat(\n",
    "            b_i.reshape(-1,1),arr_train.shape[1],axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Baseline Model Mean Squared Error: 0.94\n"
     ]
    }
   ],
   "source": [
    "#use mse as benchmark for comparison\n",
    "mse_baseline=np.square((\n",
    "    output-arr)[testing_idx]).sum()/len(arr[testing_idx])\n",
    "print('Baseline Model Mean Squared Error:',round(mse_baseline,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Pearson Correlation with baseline\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Since we are using baseline in KNN, it's also worth a while to switch the distance metric from MSD to Pearson Correlation with Baseline. This metric is a variation of Pearson Correlation in statistics.\n",
    "\n",
    "$$\\rho_{uv} = \\frac{cov(u,v)}{\\sigma_u \\cdot \\sigma_v} = \\frac{ \\sum\\limits_{i \\in I_{uv}}\n",
    "        (r_{ui} -  \\mu_u) \\cdot (r_{vi} - \\mu_{v})} {\\sqrt{\\sum\\limits_{i\n",
    "        \\in I_{uv}} (r_{ui} -  \\mu_u)^2} \\cdot \\sqrt{\\sum\\limits_{i \\in\n",
    "        I_{uv}} (r_{vi} -  \\mu_{v})^2} }$$\n",
    "        \n",
    "To replace $\\mu_u$ and $\\mu_{v}$ with $b_u$ and $b_{v}$ in the equation, we obtain the designated metric – Pearson Correlation with Baseline.\n",
    "\n",
    "$$\\text{similarity}(u, v) = \\frac{\n",
    "            \\sum\\limits_{i \\in I_{uv}} (r_{ui} -  b_{u}) \\cdot (r_{vi} -\n",
    "            b_{v})} {\\sqrt{\\sum\\limits_{i \\in I_{uv}} (r_{ui} -  b_{u})^2}\n",
    "            \\cdot \\sqrt{\\sum\\limits_{i \\in I_{uv}} (r_{vi} -  b_{v})^2}}$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "#obtain pearson correlation with baseline\n",
    "#if b_u is substituted with arr.mean(axis=1) and arr is a complete matrix\n",
    "#the function will be equivalent to np.corrcoef(arr.T)\n",
    "#except the diagonal line would be filled with zeros\n",
    "@numba.njit\n",
    "def get_pearson_corr_baseline_matrix(arr,b_u):\n",
    "\n",
    "    similarity_matrix=np.zeros((arr.shape[1],arr.shape[1]))\n",
    "    for u in range(arr.shape[1]):\n",
    "        for v in range(u+1,arr.shape[1]):\n",
    "            \n",
    "            #self correlation distorts knn selection\n",
    "            if u==v:\n",
    "                continue\n",
    "\n",
    "            #find items rated by both user u and v\n",
    "            arr_sub1=arr[:,u]\n",
    "            arr_sub2=arr[:,v]\n",
    "            set_u=np.where(~np.isnan(arr_sub1))\n",
    "            set_v=np.where(~np.isnan(arr_sub2))\n",
    "            I_uv=set(set_u[0]).intersection(set(set_v[0]))\n",
    "          \n",
    "            #avoid the case where two users havent rated any items in common\n",
    "            if len(I_uv)>0:\n",
    "                \n",
    "                #extract ratings of inner join items from user u and v\n",
    "                arr_u=arr_sub1[np.array(list(I_uv))]\n",
    "                arr_v=arr_sub2[np.array(list(I_uv))]\n",
    "                numerator=(arr_u-b_u[u]).T@(arr_v-b_u[v])\n",
    "                denominator=np.sqrt(\n",
    "                    (arr_u-b_u[u]).T@(arr_u-b_u[u]))*np.sqrt(\n",
    "                    (arr_v-b_u[v]).T@(arr_v-b_u[v]))\n",
    "                pearson_baseline=numerator/denominator\n",
    "            else:\n",
    "                pearson_baseline=0\n",
    "\n",
    "            #symmetric matrix\n",
    "            similarity_matrix[u,v]=pearson_baseline\n",
    "            similarity_matrix[v,u]=pearson_baseline\n",
    "            \n",
    "    return similarity_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute similarity matrix\n",
    "similarity_matrix=get_pearson_corr_baseline_matrix(arr,b_u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### KNN with baseline\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "There is nothing trivial about KNN with baseline. Like Pearson Correlation with baseline, we replace $\\mu_u$ and $\\mu_v$ with $b_u$ and $b_v$ in KNN with mean. Ça y est!\n",
    "\n",
    "$$\\hat{r}_{ui} = b_{u} + \\frac{  \\sum_{v \\in N^k_i(u)} \\mathcal{similarity}(u, v) \\cdot (r_{vi}-b_{v})}  {\\sum_{v \\in N^k_i(u)} \\mathcal{similarity}(u, v)}$$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Reference to the paper with KNN with baseline (equation 3)\n",
    "\n",
    "https://people.engr.tamu.edu/huangrh/Spring16/papers_course/matrix_factorization.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "#knn with baseline for matrix completion\n",
    "@numba.njit\n",
    "def knn_with_baseline(arr,similarity_matrix,\n",
    "                      b_u,top_k=10):\n",
    "    \n",
    "    estimated=arr.copy()\n",
    "    for i in range(arr.shape[0]):\n",
    "        for u in range(arr.shape[1]):\n",
    "            if np.isnan(arr[i,u]):\n",
    "\n",
    "                #find top k neighbor based upon similarity matrix\n",
    "                rated_users=np.where(~np.isnan(arr[i]))[0]\n",
    "                similarities=similarity_matrix[u][rated_users]\n",
    "                top_k_neighbors=np.argsort(similarities)[-top_k:]\n",
    "                N_k_i_u=np.array([\n",
    "                    rated_users[neighbor] for neighbor in top_k_neighbors])\n",
    "\n",
    "                #compute weighted average                \n",
    "                if len(N_k_i_u)!=0:\n",
    "                    numerator=(arr[i][\n",
    "                        N_k_i_u]-b_u[\n",
    "                        N_k_i_u]).T@similarity_matrix[u][N_k_i_u]\n",
    "                    denominator=similarity_matrix[u][N_k_i_u].sum()\n",
    "                    if denominator!=0:\n",
    "                        estimated[i,u]=b_u[u]+numerator/denominator\n",
    "                        continue\n",
    "                        \n",
    "                #when the users who rated item i\n",
    "                #have nothing in common with user u\n",
    "                #take the average of user baseline and item baseline\n",
    "                #if both user baseline and item baseline are empty\n",
    "                #take global baseline\n",
    "                miu_i=arr[i][~np.isnan(arr[i])]\n",
    "                miu_u=arr[:,u][~np.isnan(arr[:,u])]                \n",
    "                if len(miu_i)==0 and len(miu_u)==0:\n",
    "                    estimated[i,u]=arr[np.where(~np.isnan(arr))[0]][\n",
    "                        np.where(~np.isnan(arr))[1]].mean()\n",
    "                else:\n",
    "                    estimated[i,u]=np.array(list(miu_i)+list(miu_u)).mean()\n",
    "                    \n",
    "    return estimated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix completion\n",
    "output=knn_with_baseline(arr_train,similarity_matrix,b_u,top_k=top_k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KNN with baseline Mean Squared Error: 0.85\n"
     ]
    }
   ],
   "source": [
    "#use mse as benchmark for comparison\n",
    "mse_knn_with_baseline=np.square((\n",
    "    output-arr)[testing_idx]).sum()/len(arr[testing_idx])\n",
    "print('KNN with baseline Mean Squared Error:',round(mse_knn_with_baseline,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Koren Neighborhood Model\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "To facilitate global optimization, Koren proposed a neighborhood model with global weights independent of a specific user. The similarity weight is solved via optimization rather than correlation matrix. However, the model below is a bit different from the original version. As the only implicit feedback we get from the data is the items rated by the users, there is no point to create an extra variable which is merely a linear combination of neighbor-weighted deviation from the baseline. The official optimization problem is formulated as below. The notations are almost the same with everything above.\n",
    "\n",
    "$$ \\min_{w_*,b_*}\\,\\sum_{r_{ui}\\,\\in\\,\\mathcal{K}} \\left(r_{ui} - \\mu - b_u - b_i - |N_i^k(u)|^{-\\frac{1}{2}} \\sum_{j \\in N_i^k(u)} (r_{uj} - \\mu - b_u - b_j)w_{ij} \\right)^2 + \\lambda\\left( \\sum_{j \\in N_i^k(u)}w_{ij}^2 + b_u^2 + b_i^2 \\right)$$\n",
    "\n",
    "where \n",
    "\n",
    "$w_{ij}$ denotes the similarity weight between item $i$ and item $j$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "The actual optimization solver is similar to Funk SVD with one more parameter $w_{ij}$. \n",
    "\n",
    "$$b_u := b_u + \\alpha (\\epsilon_{ui} - \\lambda b_u)$$\n",
    "$$b_i := b_i + \\alpha (\\epsilon_{ui} - \\lambda b_i)$$\n",
    "$$w_{ij} := w_{ij} + \\alpha (|N(u)|^{-\\frac{1}{2}} \\cdot \\epsilon_{ui} \\cdot (r_{uj} - \\mu - b_u - b_j) - \\lambda w_{ij})$$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Reference to the original paper of Koren neighborhood model (equation 11)\n",
    "\n",
    "https://people.engr.tamu.edu/huangrh/Spring16/papers_course/matrix_factorization.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use numba to dramatically boost the speed of linear algebra\n",
    "@numba.njit\n",
    "def koren_neighbor_epoch(arr,miu,b_u,b_i,\n",
    "                         similarity_matrix,\n",
    "                         w_ij,alpha,lambda_,top_k):\n",
    "    \n",
    "    #initialize\n",
    "    error=0\n",
    "    \n",
    "    #only iterate known ratings\n",
    "    for i in range(arr.shape[0]):\n",
    "        for u in range(arr.shape[1]):\n",
    "            r_ui=arr[i,u]\n",
    "            \n",
    "            #another way to identify nan\n",
    "            #r_ui!=r_ui\n",
    "            if np.isnan(r_ui):\n",
    "                continue\n",
    "                \n",
    "            #find top k neighbor based upon similarity matrix\n",
    "            rated_items=np.where(~np.isnan(arr[:,u]))[0]\n",
    "            similarities=similarity_matrix[i][rated_items]\n",
    "            top_k_neighbors=np.argsort(similarities)[-top_k:]\n",
    "            N_k_i_u=np.array([\n",
    "                    rated_items[neighbor] for neighbor in top_k_neighbors])\n",
    "            \n",
    "            #compute error\n",
    "            if len(N_k_i_u)!=0:\n",
    "                deviation=arr[:,u][N_k_i_u]-miu-b_u[u]-b_i[N_k_i_u]\n",
    "                weighted_sum=(w_ij[i][N_k_i_u]).T@deviation\n",
    "                epsilon_ui=(r_ui-miu-b_u[u]-b_i[i]-weighted_sum).item()\n",
    "                error+=epsilon_ui**2\n",
    "            else:\n",
    "                epsilon_ui=(r_ui-miu-b_u[u]-b_i[i]).item()\n",
    "                error+=epsilon_ui**2\n",
    "\n",
    "            #update baseline\n",
    "            b_u[u]+=alpha*(epsilon_ui-lambda_*b_u[u])\n",
    "            b_i[i]+=alpha*(epsilon_ui-lambda_*b_i[i])\n",
    "            \n",
    "            #only update weights when there are similar items                \n",
    "            if len(N_k_i_u)!=0:\n",
    "                N_k_i_u_norm_sqrt=N_k_i_u.shape[0]**0.5                \n",
    "                w_ij[i][N_k_i_u]=w_ij[i][N_k_i_u]+alpha*(\n",
    "                    epsilon_ui*deviation/N_k_i_u_norm_sqrt-lambda_*w_ij[\n",
    "                        i][N_k_i_u])\n",
    "                w_ij[N_k_i_u][:,i]=w_ij[i][N_k_i_u]\n",
    "                                        \n",
    "    return error,b_u,b_i,w_ij"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "#koren weighted neighborhood model\n",
    "def koren_neighbor(arr,similarity_matrix,\n",
    "                   miu_init=None,b_u_init=[],\n",
    "                   b_i_init=[],w_ij_init=[],\n",
    "                   alpha=0.005,lambda_=0.02,\n",
    "                   tau=0.0001,top_k=10,\n",
    "                   max_iter=20,diagnosis=True):\n",
    "\n",
    "    #initialize\n",
    "    stop=False\n",
    "    counter=0\n",
    "    sse=None\n",
    "    \n",
    "    #global mean\n",
    "    if not miu_init:       \n",
    "        miu=arr[~np.isnan(arr)].mean()\n",
    "    else:\n",
    "        miu=miu_init\n",
    "        \n",
    "    #user baseline\n",
    "    if len(b_u_init)==0:\n",
    "        b_u=np.zeros(arr.shape[1])\n",
    "    else:\n",
    "        b_u=b_u_init\n",
    "    \n",
    "    #item baseline\n",
    "    if len(b_i_init)==0:\n",
    "        b_i=np.zeros(arr.shape[0])\n",
    "    else:\n",
    "        b_i=b_i_init\n",
    "        \n",
    "    #weighted neighbor\n",
    "    if len(w_ij_init)==0:\n",
    "        w_ij=np.zeros((arr.shape[0],arr.shape[0]))\n",
    "        w_ij.fill(0.001)\n",
    "    else:\n",
    "        w_ij=w_ij_init\n",
    "    \n",
    "    #gradient descent\n",
    "    while not stop:\n",
    "        \n",
    "        error,b_u,b_i,w_ij=koren_neighbor_epoch(\n",
    "                         arr,miu,b_u,b_i,\n",
    "                         similarity_matrix,\n",
    "                         w_ij,alpha,lambda_,top_k)\n",
    "\n",
    "        counter+=1\n",
    "        \n",
    "        #maximum number of epoch\n",
    "        if counter>=max_iter:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print('Not converged.',\n",
    "                      'Consider increase number of iterations or tolerance')\n",
    "                \n",
    "        #use sum of squared error to determine if converged\n",
    "        sse_prev=sse\n",
    "        sse=error\n",
    "        if sse_prev and abs(sse/sse_prev-1)<=tau:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print(f'{counter} iterations to reach convergence\\n')\n",
    "\n",
    "    return b_u,b_i,w_ij"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute item similarity matrix\n",
    "similarity_matrix=get_msd_similarity_matrix(\n",
    "    arr_train.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "310 iterations to reach convergence\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#koren neighbor\n",
    "b_u,b_i,w_ij=koren_neighbor(arr_train,\n",
    "                   similarity_matrix,\n",
    "                   alpha=learning_rate,\n",
    "                   lambda_=lagrange_multiplier,\n",
    "                   tau=tolerance,top_k=top_k,\n",
    "                   max_iter=max_num_of_epoch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute deviation\n",
    "deviation=arr_train-np.repeat(\n",
    "            b_u.reshape(1,-1),\n",
    "            arr_train.shape[0],axis=0)+np.repeat(\n",
    "            b_i.reshape(-1,1),\n",
    "    arr_train.shape[1],axis=1)-miu\n",
    "\n",
    "#set nan to zero for dot product\n",
    "deviation[np.isnan(deviation)]=0\n",
    "\n",
    "#find top k neighbors for each item\n",
    "top_k_neighbors=np.argsort(similarity_matrix,axis=0)[-top_k:]\n",
    "\n",
    "#identify the col index of top k neighbors\n",
    "col=top_k_neighbors.flatten()\n",
    "\n",
    "#identify the row index of each item\n",
    "row=np.array([j for i in range(arr_train.shape[0]) for j in [i]*top_k])\n",
    "\n",
    "#create weights matrix\n",
    "weights=np.zeros(w_ij.shape)\n",
    "\n",
    "#only keep weights of top k neighbors\n",
    "weights[(row,col)]=w_ij[(row,col)]\n",
    "\n",
    "#compute influence from weighted neighbors\n",
    "weighted_neighbors=weights@deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix completion\n",
    "output=miu+np.repeat(\n",
    "            b_u.reshape(1,-1),\n",
    "            arr_train.shape[0],axis=0)+np.repeat(\n",
    "            b_i.reshape(-1,1),\n",
    "    arr_train.shape[1],axis=1)+weighted_neighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Koren Neighborhood Model Mean Squared Error: 1.141\n"
     ]
    }
   ],
   "source": [
    "#use mse as benchmark for comparison\n",
    "mse_koren_ngbr=np.square((\n",
    "    output-arr)[testing_idx]).sum()/len(arr[testing_idx])\n",
    "print('Koren Neighborhood Model Mean Squared Error:',round(mse_koren_ngbr,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Koren Integrated Model\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "There is nothing special about Koren Integrated Model. It is a combination of SVD++ and Koren Neighborhood Model which intends to get the best from both worlds. The official optimization problem is the aggregation of both.\n",
    "\n",
    "$$ \\min_{w_*,p_*,q_*,b_*,y_*}\\,\\sum_{r_{ui}\\,\\in\\,\\mathcal{K}} \\left(r_{ui} - \\mu - b_u - b_i - |N_i^k(u)|^{-\\frac{1}{2}} \\sum_{j \\in N_i^k(u)} (r_{uj} - \\mu - b_u - b_j)w_{ij} - q_i^T\\left(p_u+|N(u)|^{-\\frac{1}{2}} \\sum_{j \\in N(u)}y_j\\right) \\right)^2 + \\lambda\\left( ||p_u||^2 + ||q_i||^2 + \\sum_{j \\in N(u)}||y_j||^2 + \\sum_{j \\in N_i^k(u)}w_{ij}^2 + b_u^2 + b_i^2 \\right)$$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "The actual optimization solver is the aggregation of both as well. \n",
    "\n",
    "$$b_u := b_u + \\alpha (\\epsilon_{ui} - \\lambda b_u)$$\n",
    "$$b_i := b_i + \\alpha (\\epsilon_{ui} - \\lambda b_i)$$\n",
    "$$w_{ij} := w_{ij} + \\alpha (|N(u)|^{-\\frac{1}{2}} \\cdot \\epsilon_{ui} \\cdot (r_{uj} - \\mu - b_u - b_j) - \\lambda w_{ij})$$\n",
    "$$p_u := p_u + \\alpha (\\epsilon_{ui} \\cdot q_i - \\lambda p_u)$$\n",
    "$$q_i := q_i + \\alpha (\\epsilon_{ui} \\cdot (p_u+|N(u)|^{-\\frac{1}{2}} \\sum_{j \\in N(u)}y_j) - \\lambda q_i)$$\n",
    "$$y_j := y_j + \\alpha (|N(u)|^{-\\frac{1}{2}} \\cdot \\epsilon_{ui} \\cdot q_i - \\lambda y_j)$$\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Reference to the original paper of Koren integrated model (equation 16)\n",
    "\n",
    "https://people.engr.tamu.edu/huangrh/Spring16/papers_course/matrix_factorization.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use numba to dramatically boost the speed of linear algebra\n",
    "@numba.njit\n",
    "def koren_integrated_epoch(arr,similarity_matrix,\n",
    "                           miu,b_u,b_i,p_u,q_i,y_j,\n",
    "                           w_ij,alpha,lambda_,top_k):\n",
    "    \n",
    "    #initialize\n",
    "    error=0\n",
    "    \n",
    "    #only iterate known ratings\n",
    "    for i in range(arr.shape[0]):\n",
    "        for u in range(arr.shape[1]):\n",
    "            r_ui=arr[i,u]\n",
    "            \n",
    "            #another way to identify nan\n",
    "            #r_ui!=r_ui\n",
    "            if np.isnan(r_ui):\n",
    "                continue\n",
    "                \n",
    "            #compute implicit feedback\n",
    "            N_u=np.where(~np.isnan(arr[:,u]))[0]\n",
    "            N_u_norm_sqrt=arr[:,u][~np.isnan(arr[:,u])].shape[0]**0.5\n",
    "            feedback=(y_j[N_u,:]/N_u_norm_sqrt).sum(axis=0)\n",
    "\n",
    "            #find top k neighbor based upon similarity matrix\n",
    "            rated_items=np.where(~np.isnan(arr[:,u]))[0]\n",
    "            similarities=similarity_matrix[i][rated_items]\n",
    "            top_k_neighbors=np.argsort(similarities)[-top_k:]\n",
    "            N_k_i_u=np.array([\n",
    "                    rated_items[neighbor] for neighbor in top_k_neighbors])\n",
    "            \n",
    "            #compute error\n",
    "            if len(N_k_i_u)!=0:\n",
    "                deviation=arr[:,u][N_k_i_u]-miu-b_u[u]-b_i[N_k_i_u]\n",
    "                weighted_sum=(w_ij[i][N_k_i_u]).T@deviation\n",
    "                epsilon_ui=(r_ui-miu-b_u[u]-b_i[\n",
    "                    i]-weighted_sum-q_i[i].T@(\n",
    "                feedback+p_u[u]).reshape(-1,1)).item()\n",
    "                error+=epsilon_ui**2\n",
    "            else:\n",
    "                epsilon_ui=(r_ui-miu-b_u[u]-b_i[\n",
    "                    i]-q_i[i].T@(\n",
    "                feedback+p_u[u]).reshape(-1,1)).item()\n",
    "                error+=epsilon_ui**2\n",
    "\n",
    "            #update\n",
    "            b_u[u]+=alpha*(epsilon_ui-lambda_*b_u[u])\n",
    "            b_i[i]+=alpha*(epsilon_ui-lambda_*b_i[i])\n",
    "            p_u[u]+=alpha*(epsilon_ui*q_i[i]-lambda_*p_u[u])\n",
    "            q_i[i]+=alpha*(epsilon_ui*(p_u[u]+feedback)-lambda_*q_i[i])\n",
    "            y_j[N_u]+=alpha*(epsilon_ui*q_i[N_u]/N_u_norm_sqrt-lambda_*y_j[N_u])\n",
    "            \n",
    "            #only update weights when there are similar items                \n",
    "            if len(N_k_i_u)!=0:\n",
    "                N_k_i_u_norm_sqrt=N_k_i_u.shape[0]**0.5                \n",
    "                w_ij[i][N_k_i_u]=w_ij[i][N_k_i_u]+alpha*(\n",
    "                    epsilon_ui*deviation/N_k_i_u_norm_sqrt-lambda_*w_ij[\n",
    "                        i][N_k_i_u])\n",
    "                w_ij[N_k_i_u][:,i]=w_ij[i][N_k_i_u]\n",
    "                  \n",
    "    return error,b_u,b_i,p_u,q_i,y_j,w_ij"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "#koren integrated model\n",
    "def koren_integrated(arr,similarity_matrix,\n",
    "                     miu_init=None,b_u_init=[],\n",
    "                     b_i_init=[],p_u_init=[],\n",
    "                     q_i_init=[],y_j_init=[],\n",
    "                     w_ij_init=[],\n",
    "                     num_of_rank=40,alpha=0.005,\n",
    "                     lambda_=0.02,tau=0.0001,top_k=10,\n",
    "                     max_iter=20,diagnosis=True\n",
    "                     ):\n",
    "\n",
    "    #initialize\n",
    "    stop=False\n",
    "    counter=0\n",
    "    sse=None\n",
    "    \n",
    "    #global mean\n",
    "    if not miu_init:       \n",
    "        miu=arr[~np.isnan(arr)].mean()\n",
    "    else:\n",
    "        miu=miu_init\n",
    "        \n",
    "    #user baseline\n",
    "    if len(b_u_init)==0:\n",
    "        b_u=np.zeros(arr.shape[1])\n",
    "    else:\n",
    "        b_u=b_u_init\n",
    "    \n",
    "    #item baseline\n",
    "    if len(b_i_init)==0:\n",
    "        b_i=np.zeros(arr.shape[0])\n",
    "    else:\n",
    "        b_i=b_i_init\n",
    "        \n",
    "    #user latent factors\n",
    "    if len(p_u_init)==0:\n",
    "        p_u=np.zeros((arr.shape[1],num_of_rank))\n",
    "        p_u.fill(0.1)\n",
    "    else:\n",
    "        p_u=p_u_init\n",
    "    \n",
    "    #item latent factors\n",
    "    if len(q_i_init)==0:\n",
    "        q_i=np.zeros((arr.shape[0],num_of_rank))\n",
    "        q_i.fill(0.1)\n",
    "    else:\n",
    "        q_i=q_i_init\n",
    "        \n",
    "    #user implicit feedback\n",
    "    if len(y_j_init)==0:\n",
    "        y_j=np.zeros((arr.shape[0],num_of_rank))\n",
    "        y_j.fill(0.1)\n",
    "    else:\n",
    "        y_j=y_j_init\n",
    "        \n",
    "    #weighted neighbor\n",
    "    if len(w_ij_init)==0:\n",
    "        w_ij=np.zeros((arr.shape[0],arr.shape[0]))\n",
    "        w_ij.fill(0.001)\n",
    "    else:\n",
    "        w_ij=w_ij_init\n",
    "    \n",
    "    #gradient descent\n",
    "    while not stop:\n",
    "        \n",
    "        error,b_u,b_i,p_u,q_i,y_j,w_ij=koren_integrated_epoch(\n",
    "                           arr,similarity_matrix,\n",
    "                           miu,b_u,b_i,p_u,q_i,y_j,\n",
    "                           w_ij,alpha,lambda_,top_k)\n",
    "\n",
    "        counter+=1\n",
    "        \n",
    "        #maximum number of epoch\n",
    "        if counter>=max_iter:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print('Not converged.',\n",
    "                      'Consider increase number of iterations or tolerance')\n",
    "                \n",
    "        #use sum of squared error to determine if converged\n",
    "        sse_prev=sse\n",
    "        sse=error\n",
    "        if sse_prev and abs(sse/sse_prev-1)<=tau:\n",
    "            stop=True\n",
    "            if diagnosis:\n",
    "                print(f'{counter} iterations to reach convergence\\n')\n",
    "\n",
    "    return b_u,b_i,p_u,q_i,y_j,w_ij"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "292 iterations to reach convergence\n",
      "\n"
     ]
    }
   ],
   "source": [
    "b_u,b_i,p_u,q_i,y_j,w_ij=koren_integrated(\n",
    "                 arr_train,similarity_matrix,\n",
    "                 num_of_rank=num_of_latent_factors,\n",
    "                 alpha=learning_rate,\n",
    "                 lambda_=lagrange_multiplier,tau=tolerance,\n",
    "                 max_iter=max_num_of_epoch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute implicit feedback\n",
    "feedback=[]\n",
    "for u in range(arr_train.shape[1]):\n",
    "    N_u=np.where(~np.isnan(arr_train[:,u]))[0]\n",
    "    N_u_norm_sqrt=arr_train[:,u][~np.isnan(arr_train[:,u])].shape[0]**0.5\n",
    "    feedback.append((y_j[N_u,:]/N_u_norm_sqrt).sum(axis=0))\n",
    "\n",
    "#compute deviation\n",
    "deviation=arr_train-np.repeat(\n",
    "            b_u.reshape(1,-1),\n",
    "            arr_train.shape[0],axis=0)+np.repeat(\n",
    "            b_i.reshape(-1,1),\n",
    "    arr_train.shape[1],axis=1)-miu\n",
    "\n",
    "#set nan to zero for dot product\n",
    "deviation[np.isnan(deviation)]=0\n",
    "\n",
    "#find top k neighbors for each item\n",
    "top_k_neighbors=np.argsort(similarity_matrix,axis=0)[-top_k:]\n",
    "\n",
    "#identify the col index of top k neighbors\n",
    "col=top_k_neighbors.flatten()\n",
    "\n",
    "#identify the row index of each item\n",
    "row=np.array([j for i in range(arr_train.shape[0]) for j in [i]*top_k])\n",
    "\n",
    "#create weights matrix\n",
    "weights=np.zeros(w_ij.shape)\n",
    "\n",
    "#only keep weights of top k neighbors\n",
    "weights[(row,col)]=w_ij[(row,col)]\n",
    "\n",
    "#compute influence from weighted neighbors\n",
    "weighted_neighbors=weights@deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix completion\n",
    "output=miu+np.repeat(\n",
    "            b_u.reshape(1,-1),\n",
    "            arr_train.shape[0],axis=0)+np.repeat(\n",
    "            b_i.reshape(-1,1),\n",
    "    arr_train.shape[1],axis=1)+q_i@(\n",
    "    p_u+np.mat(feedback)).T+weighted_neighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Koren Integrated Model Mean Squared Error: 1.143\n"
     ]
    }
   ],
   "source": [
    "#use mse as benchmark for comparison\n",
    "mse_koren_integrated=np.square((\n",
    "    output-arr)[testing_idx]).sum()/len(arr[testing_idx])\n",
    "print('Koren Integrated Model Mean Squared Error:',\n",
    "      round(mse_koren_integrated,3))"
   ]
  },
  {
   "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.6.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "154px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}