<jupyter_start><jupyter_text>CF Part 6 - Non-negative Matrix factorization method> Collaborative Filtering on MovieLens Latest-small Part 6 - Finding recommendations using model based non-negative matrix factorization method- toc: true- badges: true- comments: true- categories: [movie, collaborative]- image: Jusl like Matrix Factorization (MF) [(<NAME> et al., 2009)](https://ieeexplore.ieee.org/document/5197422), Non-negative Matrix Factorization (NMF in short) factors the rating matrix $R$ in two matrices in such a way that $R = PQ^{\top}$. One limitation of Matrix Factorization$P$ and $Q$ values in MF are non interpretable since their components can take arbitrary (positive and negative) values. Particulariy of Non-negative Matrix FactorizationNMF [(Lee and Seung, 1999)](http://www.dm.unibo.it/~simoncin/nmfconverge.pdf) allows the reconstruction of $P$ and $Q$ in such a way that $P,Q \ge 0$. Constraining $P$ and $Q$ values to be taken from $[0,1]$ allows  a probabilistic interpretation- Latent factors represent groups of users who share the same tastes,- The value  $P_{u,l}$  represents the probability that user $u$ belongs to the group $l$ of users and - The value $Q_{l,i}$ represents the probability that users in the group $l$  likes item $i$. Objective functionWith the Euclidian distance, the NMF objective function is defined by\begin{equation}J = \frac{1}{2}\sum_{(u,i) \in \kappa}||R_{u,i} - P_uQ_i^{\top}||^2 + \lambda_P||P_u||^2 + \lambda_Q||Q_i||^2\end{equation}The goal is to minimize the cost function $J$ by optimizing parameters $P$ and $Q$, with $\lambda_P$ and $\lambda_Q$ the regularizer parameters.  Multiplicative update ruleAccording [(Lee and Seung, 1999)](https://www.nature.com/articles/44565.pdf), to the multiplicative update rule for $P$ and $Q$ are as follows :\begin{equation}P \leftarrow P \cdot \frac{RQ}{PQ^{\top}Q}\end{equation}\begin{equation}Q \leftarrow Q \cdot \frac{R^{\top}P}{QP^{\top}P}\end{equation}However, since $R$ is a sparse matrix, we need to update each $P_u$ according to existing ratings of user $u$. Similarily, we need to update $Q_i$ according to existing ratings on item $i$. Hence :\begin{equation}P_{u,k} \leftarrow P_{u,k} \cdot \frac{\sum_{i \in I_u}Q_{i,k}\cdot r_{u,i}}{\sum_{i \in I_u}Q_{i,k}\cdot \hat{r}_{u,i} + \lambda_P|I_u|P_{u,k}}\end{equation}\begin{equation}Q_{i,k} \leftarrow Q_{i,k} \cdot \frac{\sum_{u \in U_i}P_{u,k}\cdot r_{u,i}}{\sum_{u \in U_i}P_{u,k}\cdot \hat{r}_{u,i} + \lambda_Q|U_i|Q_{i,k}}\end{equation}Where- $P_{u,k}$ is the $k^{th}$ latent factor of $P_u$- $Q_{i,k}$ is the $k^{th}$ latent factor of $Q_i$- $I_u$ the of items rated by user $u$- $U_i$ the set of users who rated item $i$<jupyter_code>import os

if not (os.path.exists("recsys.zip") or os.path.exists("recsys")):
    !wget https://github.com/nzhinusoftcm/review-on-collaborative-filtering/raw/master/recsys.zip    
    !unzip recsys.zip<jupyter_output><empty_output><jupyter_text>Install and import useful packages  requirements```matplotlib==3.2.2numpy==1.19.2pandas==1.0.5python==3.7scikit-learn==0.24.1scikit-surprise==1.1.1scipy==1.6.2```<jupyter_code>from recsys.preprocessing import mean_ratings
from recsys.preprocessing import normalized_ratings
from recsys.preprocessing import ids_encoder
from recsys.preprocessing import train_test_split
from recsys.preprocessing import rating_matrix
from recsys.preprocessing import get_examples
from recsys.preprocessing import scale_ratings

from recsys.datasets import ml1m
from recsys.datasets import ml100k

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

import os<jupyter_output><empty_output><jupyter_text>Load and preprocess rating<jupyter_code># load data
ratings, movies = ml100k.load()

# prepare data
ratings, uencoder, iencoder = ids_encoder(ratings)

# convert ratings from dataframe to numpy array
np_ratings = ratings.to_numpy()

# get examples as tuples of userids and itemids and labels from normalize ratings
raw_examples, raw_labels = get_examples(ratings, labels_column="rating")

# train test split
(x_train, x_test), (y_train, y_test) = train_test_split(examples=raw_examples, labels=raw_labels)<jupyter_output><empty_output><jupyter_text>Non-negative Matrix Factorization Model<jupyter_code>class NMF:
    
    def __init__(self, ratings, m, n, uencoder, iencoder, K=10, lambda_P=0.01, lambda_Q=0.01):
        
        np.random.seed(32)
        
        # initialize the latent factor matrices P and Q (of shapes (m,k) and (n,k) respectively) that will be learnt
        self.ratings = ratings
        self.np_ratings = ratings.to_numpy()
        self.K = K
        self.P = np.random.rand(m, K)
        self.Q = np.random.rand(n, K)
        
        # hyper parameter initialization
        self.lambda_P = lambda_P
        self.lambda_Q = lambda_Q

        # initialize encoders
        self.uencoder = uencoder
        self.iencoder = iencoder
        
        # training history
        self.history = {
            "epochs": [],
            "loss": [],
            "val_loss": [],
        }
    
    def print_training_parameters(self):
        print('Training NMF ...')
        print(f'k={self.K}')
        
    def mae(self, x_train, y_train):
        """
        returns the Mean Absolute Error
        """
        # number of training examples
        m = x_train.shape[0]
        error = 0
        for pair, r in zip(x_train, y_train):
            u, i = pair
            error += abs(r - np.dot(self.P[u], self.Q[i]))
        return error / m
    
    def update_rule(self, u, i, error):
        I = self.np_ratings[self.np_ratings[:, 0] == u][:, [1, 2]]
        U = self.np_ratings[self.np_ratings[:, 1] == i][:, [0, 2]]    
                    
        num = self.P[u] * np.dot(self.Q[I[:, 0]].T, I[:, 1])
        dem = np.dot(self.Q[I[:, 0]].T, np.dot(self.P[u], self.Q[I[:, 0]].T)) + self.lambda_P * len(I) * self.P[u]
        self.P[u] = num / dem

        num = self.Q[i] * np.dot(self.P[U[:, 0]].T, U[:, 1])
        dem = np.dot(self.P[U[:, 0]].T, np.dot(self.P[U[:, 0]], self.Q[i].T)) + self.lambda_Q * len(U) * self.Q[i]
        self.Q[i] = num / dem
    
    @staticmethod
    def print_training_progress(epoch, epochs, error, val_error, steps=5):
        if epoch == 1 or epoch % steps == 0:
            print(f"epoch {epoch}/{epochs} - loss : {round(error, 3)} - val_loss : {round(val_error, 3)}")
                
    def fit(self, x_train, y_train, validation_data, epochs=10):

        self.print_training_parameters()
        x_test, y_test = validation_data
        for epoch in range(1, epochs+1):
            for pair, r in zip(x_train, y_train):
                u, i = pair
                r_hat = np.dot(self.P[u], self.Q[i])
                e = abs(r - r_hat)
                self.update_rule(u, i, e)                
            # training and validation error  after this epochs
            error = self.mae(x_train, y_train)
            val_error = self.mae(x_test, y_test)
            self.update_history(epoch, error, val_error)
            self.print_training_progress(epoch, epochs, error, val_error, steps=1)
        
        return self.history
    
    def update_history(self, epoch, error, val_error):
        self.history['epochs'].append(epoch)
        self.history['loss'].append(error)
        self.history['val_loss'].append(val_error)
    
    def evaluate(self, x_test, y_test):        
        error = self.mae(x_test, y_test)
        print(f"validation error : {round(error,3)}")
        print('MAE : ', error)        
        return error
      
    def predict(self, userid, itemid):
        u = self.uencoder.transform([userid])[0]
        i = self.iencoder.transform([itemid])[0]
        r = np.dot(self.P[u], self.Q[i])
        return r

    def recommend(self, userid, N=30):
        # encode the userid
        u = self.uencoder.transform([userid])[0]

        # predictions for users userid on all product
        predictions = np.dot(self.P[u], self.Q.T)

        # get the indices of the top N predictions
        top_idx = np.flip(np.argsort(predictions))[:N]

        # decode indices to get their corresponding itemids
        top_items = self.iencoder.inverse_transform(top_idx)

        # take corresponding predictions for top N indices
        preds = predictions[top_idx]

        return top_items, preds<jupyter_output><empty_output><jupyter_text>Train the NMF model with ML-100K datasetmodel parameters :- $k = 10$ : (number of factors)- $\lambda_P = 0.6$- $\lambda_Q = 0.6$- epochs = 10Note that it may take some time to complete the training on 10 epochs (around 7 minutes).<jupyter_code>m = ratings['userid'].nunique()   # total number of users
n = ratings['itemid'].nunique()   # total number of items

# create and train the model
nmf = NMF(ratings, m, n, uencoder, iencoder, K=10, lambda_P=0.6, lambda_Q=0.6)
history = nmf.fit(x_train, y_train, epochs=10, validation_data=(x_test, y_test))
nmf.evaluate(x_test, y_test)<jupyter_output>validation error : 0.917
MAE :  0.916504134301954<jupyter_text>Evaluation of NMF with Scikit-suprise We can use the [scikt-suprise](https://surprise.readthedocs.io/en/stable/) package to train the NMF model. It  is an easy-to-use Python scikit for recommender systems.1. Import the NMF class from the suprise scikit.2. Load the data with the built-in function3. Instanciate NMF with ```k=10``` (```n_factors```) and we use 10 epochs (```n_epochs```)4. Evaluate the model using cross-validation with 5 folds.<jupyter_code>from surprise import NMF
from surprise import Dataset
from surprise.model_selection import cross_validate

# Load the movielens-100k dataset (download it if needed).
data = Dataset.load_builtin('ml-100k')

# Use the NMF algorithm.
nmf = NMF(n_factors=10, n_epochs=10)

# Run 5-fold cross-validation and print results.
history = cross_validate(nmf, data, measures=['MAE'], cv=5, verbose=True)<jupyter_output>Evaluating MAE of algorithm NMF on 5 split(s).

                  Fold 1  Fold 2  Fold 3  Fold 4  Fold 5  Mean    Std     
MAE (testset)     0.9626  0.9435  0.9777  0.9607  0.9510  0.9591  0.0116  
Fit time          0.63    0.69    0.74    0.65    0.71    0.68    0.04    
Test time         0.15    0.11    0.14    0.10    0.15    0.13    0.02<jupyter_text>As result, the mean MAE on the test set is **mae = 0.9510** which is equivalent to the result we have obtained on *ml-100k* with our own implementation **mae = 0.9165**  ML-1MThe mean MAE on a 5-fold cross-validation is **mae = 0.9567**This may take arount 2 minutes<jupyter_code>data = Dataset.load_builtin('ml-1m')
nmf = NMF(n_factors=10, n_epochs=10)
history = cross_validate(nmf, data, measures=['MAE'], cv=5, verbose=True)<jupyter_output>Evaluating MAE of algorithm NMF on 5 split(s).

                  Fold 1  Fold 2  Fold 3  Fold 4  Fold 5  Mean    Std     
MAE (testset)     0.9559  0.9580  0.9566  0.9639  0.9490  0.9567  0.0047  
Fit time          6.37    6.55    6.49    6.45    6.51    6.47    0.06    
Test time         1.68    1.31    1.64    1.46    1.75    1.57    0.16
{'fraction_non_alphanumeric': 0.12613101701049584, 'fraction_numerical': 0.0326637712631198, 'mean_word_length': 3.096738324684952, 'pattern_counts': {'":': 3, '<': 27, '<?xml version=': 0, '>': 28, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '575248', 'license': 'Apache-2.0', 'n_tokens_mistral': 4020, 'n_tokens_neox': 3768, 'n_words': 1068, 'path': '_notebooks/2021-06-23-collaborative-filtering-movielens-latest-small-06.ipynb', 'repo_name': 'recohut/notebook', 'size': '11058'}