# 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) [(Yehuda Koren 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 Factorization

NMF [(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 function

With 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 rule

According [(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$


```python
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
```

### Install and import useful packages

### requirements

```
matplotlib==3.2.2
numpy==1.19.2
pandas==1.0.5
python==3.7
scikit-learn==0.24.1
scikit-surprise==1.1.1
scipy==1.6.2
```


```python
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
```

### Load and preprocess rating


```python
# 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)
```

### Non-negative Matrix Factorization Model


```python
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

```

### Train the NMF model with ML-100K dataset

model parameters :

- $k = 10$ : (number of factors)
- $\lambda_P = 0.6$
- $\lambda_Q = 0.6$
- epochs = 10

Note that it may take some time to complete the training on 10 epochs (around 7 minutes).


```python
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))
```

    Training NMF ...
    k=10
    epoch 1/10 - loss : 0.916 - val_loss : 0.917
    epoch 2/10 - loss : 0.915 - val_loss : 0.917
    epoch 3/10 - loss : 0.915 - val_loss : 0.917
    epoch 4/10 - loss : 0.915 - val_loss : 0.917
    epoch 5/10 - loss : 0.915 - val_loss : 0.917
    epoch 6/10 - loss : 0.915 - val_loss : 0.917
    epoch 7/10 - loss : 0.915 - val_loss : 0.917
    epoch 8/10 - loss : 0.915 - val_loss : 0.917
    epoch 9/10 - loss : 0.915 - val_loss : 0.917
    epoch 10/10 - loss : 0.915 - val_loss : 0.917



```python
nmf.evaluate(x_test, y_test)
```

    validation error : 0.917
    MAE :  0.916504134301954





    0.916504134301954



## 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 function
3. Instanciate NMF with ```k=10``` (```n_factors```) and we use 10 epochs (```n_epochs```)
4. Evaluate the model using cross-validation with 5 folds.


```python
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)
```

    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    


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-1M

The mean MAE on a 5-fold cross-validation is **mae = 0.9567**

This may take arount 2 minutes


```python
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)
```

    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    


## Explainable Matrix Factorization

NMF introduced explainability to MF by constraining $P$ and $Q$ values in $[0,1]$. It can be considered is an inner explainable model. It's possible to inject external informations into the model in order to bring explainability. This is what the **Explainable Matrix Factorization (EMF)** algorithm does. It computes explainable scores from user or item similarities taken from the user-based or item-based CF.

Click [here](https://github.com/nzhinusoftcm/review-on-collaborative-filtering/blob/master/7.Explainable_Matrix_Factorization.ipynb) to open the EMF notebook.

## Reference

1. Daniel D. Lee & H. Sebastian Seung (1999). [Learning the parts of objects by non-negative matrix factorization](https://www.nature.com/articles/44565)
2. Deng Cai et al. (2008). [Non-negative Matrix Factorization on Manifold](https://ieeexplore.ieee.org/document/4781101)
3. Yu-Xiong Wang and Yu-Jin Zhang (2011). [Non-negative Matrix Factorization: a Comprehensive Review](https://ieeexplore.ieee.org/document/6165290)
4. Nicolas Gillis (2014). [The Why and How of Nonnegative Matrix Factorization](https://arxiv.org/pdf/1401.5226.pdf)

## Author

[Carmel WENGA](https://www.linkedin.com/in/carmel-wenga-871876178/), <br>
PhD student at Université de la Polynésie Française, <br> 
Applied Machine Learning Research Engineer, <br>
[ShoppingList](https://shoppinglist.cm), NzhinuSoft.

e07f057583c20903556304e7810657a2fec6c331
ipynb
Jupyter Notebook
_notebooks/2021-06-23-collaborative-filtering-movielens-latest-small-06.ipynb
recohut/notebook
610670666a1c3d8ef430d42f712ff72ecdbd8f86
Apache-2.0
_notebooks/2021-06-23-collaborative-filtering-movielens-latest-small-06.ipynb
recohut/notebook
610670666a1c3d8ef430d42f712ff72ecdbd8f86
Apache-2.0
2022-01-12T05:40:57.000Z
2022-01-12T05:40:57.000Z
_notebooks/2021-06-23-collaborative-filtering-movielens-latest-small-06.ipynb
recohut/notebook
610670666a1c3d8ef430d42f712ff72ecdbd8f86
Apache-2.0
2021-08-13T19:00:26.000Z
2021-08-13T19:00:26.000Z
Qwen/Qwen-72B
1. YES
2. YES