# CF Part 5 - Matrix factorization method
> Collaborative Filtering on MovieLens Latest-small Part 5 - Finding recommendations using model based  matrix factorization method

- toc: true
- badges: true
- comments: true
- categories: [movie, collaborative]
- image:

<b>User-based</b> and <b>Item-based</b> collaborative Filtering recommender systems suffer from <i>data sparsity</i> and <i>scalability</i> for online recommendations. <b>Matrix Factorization</b> helps to address these drawbacks of memory-based collaborative filtering by reducing the dimension of the rating matrix $R$.

The movielen lasted small dataset has 100k ratings of $m=610$ users on $n=9724$ items. The rating matrix in then a $m\times n$ matrix (i.e $R\in \mathbb{R}^{m\times n}$). The fact that users usually interact with less than $1\%$ of items leads the rating matrix $R$ to be highly sparse. For example, the degree of sparsity of the movielen lasted small dataset is 

\begin{equation}
sparsity = 100 - \frac{\text{total # ratings}}{m \times n} = 100 - \frac{100000}{610\times 9724} = 98,3\%
\end{equation}

This means that in this dataset, a user has interacted with less than $2\%$ of items. To reduce the dimension of the rating matrix $R$, Matrix Factorization (MF) mappes both users and items to a joint latent factor space of dimensionality $k$ such that user-item interactions are modeled as inner products in that space <a href='https://ieeexplore.ieee.org/document/5197422'>(Yehuda Koren et al., 2009)</a>. MF then decomposes $R$ in two matrices as follows :

\begin{equation}
R = Q^\top P
\end{equation}

Where $P \in \mathbb{R}^{m\times k}$ represents latent factors of users and $Q \in \mathbb{R}^{n\times k}$ is the latent factors of items. Each line of $P$, say $p_u \in \mathbb{R}^k$ denotes the taste of user $u$ and each $q_i \in \mathbb{R}^k$ the features of item $i$. The dot product between $p_u$ and $q_i$ will be the rating prediction of user $u$ on item $i$ :

\begin{equation}
\hat{r}_{u,i} = q_{i}^{\top} p_u.
\end{equation}

Figure 1 presents an example of decomposition of $R$ into two matrices $P$ and $Q$.




To learn the latent factors $p_u$ and $q_i$, the system minimizes the regularized squared error on the set of known ratings. The cost function $J$ is defined as follows : 

\begin{equation}
J = \frac{1}{2}\sum_{(u,i)\in \kappa} (r_{ui} - q_{i}^{\top} p_u)^2 + \lambda(||p_u||^2 + ||q_i||^2)
\end{equation}

where $\kappa$ is the set of $(u,i)$ pairs for which $r_{u,i}$ is known (the training set), and $\lambda$ is the regularizer parameter.

## Learning Algorithms

As described in <a href='https://ieeexplore.ieee.org/document/5197422'>(Yehuda Koren et al., 2009)</a>, to minimize the cost function $J$, the matrix factorization algorithm predicts $\hat{r}_{u,i}$ for each given training case (existing $r_{u,i}$), and computes the associated error defined by the Mean Absolute Error (MAE) as :

\begin{equation}
e_{u,i} = |r_{ui} - q_{i}^{\top} p_u|.
\end{equation}

<b>Note</b> : The overall error $E$ is defined as :

\begin{equation}
E = \frac{1}{M}\sum_{(u,i)\in\kappa} e_{u,i}
\end{equation}

Where $M$ is the number of example. The update rules for parameters $p_u$ and $q_i$ are defined as follows :

\begin{equation}
q_i \leftarrow q_i - \alpha\frac{\partial}{\partial q_i}J_{u,i}, 
\end{equation}

\begin{equation}
p_u \leftarrow p_u - \alpha\frac{\partial}{\partial p_u}J_{u,i}
\end{equation}

where $\alpha$ is the learning rate and $\frac{\partial}{\partial p_u}J_{u,i}$ is the partial derivative of the cost function $J$ according to $p_u$. It computes the extent to which $p_u$ contributes to the total error.

### How to compute $\frac{\partial}{\partial q_i}J_{u,i}$ ?

\begin{align}
\frac{\partial}{\partial q_i}J_{u,i} & = & \frac{1}{2}\frac{\partial}{\partial q_i} \begin{bmatrix}(r_{ui} - q_{i}^{\top} p_u)^2 + \lambda(||p_u||^2 + ||q_i||^2)\end{bmatrix} \\
& = & -(r_{u,i}-q_{i}^{\top} p_u)\cdot p_u + \lambda \cdot q_i  \\
& = & -e_{u,i}\cdot p_u+\lambda \cdot q_i
\end{align}

The update rules are then given by : 

\begin{equation}
q_i \leftarrow q_i + \alpha\cdot (e_{u,i}\cdot p_u-\lambda \cdot q_i), 
\end{equation}

\begin{equation}
p_u \leftarrow p_u + \alpha\cdot (e_{u,i}\cdot q_i-\lambda \cdot p_u)
\end{equation}

## Matrix Factorization : algorithm
<ol>
    <li>Initialize $P$ and $Q$ with random values
    <li>For each training example $(u,i)\in\kappa$ with the corresponding rating $r_{u,i}$ :
        <ul>
            <li>compute $\hat{r}_{u,i}$ as $\hat{r}_{u,i} = q_{i}^{\top} p_u$
            <li>compute the error : $e_{u,i} = |r_{ui} - \hat{r}_{u,i}|$
            <li>update $p_u$ and $q_i$:
                <ul>
                    <li>$p_u \leftarrow p_u + \alpha\cdot (e_{u,i}\cdot q_i-\lambda \cdot p_u)$
                    <li>$q_i \leftarrow q_i + \alpha\cdot (e_{u,i}\cdot p_u-\lambda \cdot q_i)$
                </ul>
        </ul>
    <li> Repeat step 2 until the optimal parameters are reached.
</ol>


### Download useful files


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

### Import 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 ml100k
from recsys.datasets import ml1m

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

import os
```

## Model definition


```python
class MatrixFactorization:
    
    def __init__(self, m, n, k=10, alpha=0.001, lamb=0.01):
        """
        Initialization of the model        
        : param
            - m : number of users
            - n : number of items
            - k : length of latent factor, both for users and items. 50 by default
            - alpha : learning rate. 0.001 by default
            - lamb : regularizer parameter. 0.02 by default
        """
        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.k = k
        self.P = np.random.normal(size=(m, k))
        self.Q = np.random.normal(size=(n, k))
        
        # hyperparameter initialization
        self.alpha = alpha
        self.lamb = lamb
        
        # training history
        self.history = {
            "epochs":[],
            "loss":[],
            "val_loss":[],
            "lr":[]
        }
    
    def print_training_parameters(self):
        print('Training Matrix Factorization Model ...')
        print(f'k={self.k} \t alpha={self.alpha} \t lambda={self.lamb}')
    
    def update_rule(self, u, i, error):
        self.P[u] = self.P[u] + self.alpha * (error * self.Q[i] - self.lamb * self.P[u])
        self.Q[i] = self.Q[i] + self.alpha * (error * self.P[u] - self.lamb * self.Q[i])
        
    def mae(self,  x_train, y_train):
        """
        returns the Mean Absolute Error
        """
        # number of training exemples
        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 print_training_progress(self, epoch, epochs, error, val_error, steps=5):
        if epoch == 1 or epoch % steps == 0 :
                print("epoch {}/{} - loss : {} - val_loss : {}".format(epoch, epochs, round(error,3), round(val_error,3)))
                
    def learning_rate_schedule(self, epoch, target_epochs = 20):
        if (epoch >= target_epochs) and (epoch % target_epochs == 0):
                factor = epoch // target_epochs
                self.alpha = self.alpha * (1 / (factor * 20))
                print("\nLearning Rate : {}\n".format(self.alpha))
    
    def fit(self, x_train, y_train, validation_data, epochs=1000):
        """
        Train latent factors P and Q according to the training set
        
        :param
            - x_train : training pairs (u,i) for which rating r_ui is known
            - y_train : set of ratings r_ui for all training pairs (u,i)
            - validation_data : tuple (x_test, y_test)
            - epochs : number of time to loop over the entire training set. 
            1000 epochs by default
            
        Note that u and i are encoded values of userid and itemid
        """
        self.print_training_parameters()
        
        # validation data
        x_test, y_test = validation_data
        
        # loop over the number of epochs
        for epoch in range(1, epochs+1):
            
            # for each pair (u,i) and the corresponding rating r
            for pair, r in zip(x_train, y_train):
                
                # get encoded values of userid and itemid from pair
                u,i = pair
                
                # compute the predicted rating r_hat
                r_hat = np.dot(self.P[u], self.Q[i])
                
                # compute the prediction error
                e = abs(r - r_hat)
                
                # update rules
                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)
            
            # update history
            self.history['epochs'].append(epoch)
            self.history['loss'].append(error)
            self.history['val_loss'].append(val_error)
            
            # update history
            self.update_history(epoch, error, val_error)
            
            # print training progress after each steps epochs
            self.print_training_progress(epoch, epochs, error, val_error, steps=1)
              
            # leaning rate scheduler : redure the learning rate as we go deeper in the number of epochs
            # self.learning_rate_schedule(epoch)
        
        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)
        self.history['lr'].append(self.alpha)
    
    def evaluate(self, x_test, y_test):
        """
        compute the global error on the test set        
        :param x_test : test pairs (u,i) for which rating r_ui is known
        :param y_test : set of ratings r_ui for all test pairs (u,i)
        """
        error = self.mae(x_test, y_test)
        print(f"validation error : {round(error,3)}")
        
        return error
      
    def predict(self, userid, itemid):
        """
        Make rating prediction for a user on an item
        :param userid
        :param itemid
        :return r : predicted rating
        """
        # encode user and item ids to be able to access their latent factors in
        # matrices P and Q
        u = uencoder.transform([userid])[0]
        i = iencoder.transform([itemid])[0]

        # rating prediction using encoded ids. Dot product between P_u and Q_i
        r = np.dot(self.P[u], self.Q[i])
        return r

    def recommend(self, userid, N=30):
        """
        make to N recommendations for a given user

        :return(top_items,preds) : top N items with the highest predictions 
        with their corresponding predictions
        """
        # encode the userid
        u = 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 = iencoder.inverse_transform(top_idx)

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

        return top_items, preds        
```


```python
epochs = 10
```

## MovieLens 100k

### Evaluation on raw ratings


```python
# load the ml100k dataset
ratings, movies = ml100k.load()

ratings, uencoder, iencoder = ids_encoder(ratings)

m = ratings.userid.nunique()   # total number of users
n = ratings.itemid.nunique()   # total number of items

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

# train test split
(x_train, x_test), (y_train, y_test) = train_test_split(examples=raw_examples, labels=raw_labels)
```


```python
# create the model
MF = MatrixFactorization(m, n, k=10, alpha=0.01, lamb=1.5)

# fit the model on the training set
history = MF.fit(x_train, y_train, epochs=epochs, validation_data=(x_test, y_test))
```

    Training Matrix Factorization Model ...
    k=10 	 alpha=0.01 	 lambda=1.5
    epoch 1/10 - loss : 2.734 - val_loss : 2.779
    epoch 2/10 - loss : 1.764 - val_loss : 1.794
    epoch 3/10 - loss : 1.592 - val_loss : 1.614
    epoch 4/10 - loss : 1.538 - val_loss : 1.556
    epoch 5/10 - loss : 1.515 - val_loss : 1.531
    epoch 6/10 - loss : 1.503 - val_loss : 1.517
    epoch 7/10 - loss : 1.496 - val_loss : 1.509
    epoch 8/10 - loss : 1.491 - val_loss : 1.504
    epoch 9/10 - loss : 1.488 - val_loss : 1.5
    epoch 10/10 - loss : 1.486 - val_loss : 1.497



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

    validation error : 1.497





    1.4973507972141993



### Evaluation on normalized ratings


```python
# load data
ratings, movies = ml100k.load()

ratings, uencoder, iencoder = ids_encoder(ratings)

m = ratings['userid'].nunique()   # total number of users
n = ratings['itemid'].nunique()   # total number of items

# normalize ratings by substracting means
normalized_column_name = "norm_rating"
ratings = normalized_ratings(ratings, norm_column=normalized_column_name)

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

# train test split
(x_train, x_test), (y_train, y_test) = train_test_split(examples=raw_examples, labels=raw_labels)
```


```python
# create the model
MF = MatrixFactorization(m, n, k=10, alpha=0.01, lamb=1.5)

# fit the model on the training set
history = MF.fit(x_train, y_train, epochs=epochs, validation_data=(x_test, y_test))
```

    Training Matrix Factorization Model ...
    k=10 	 alpha=0.01 	 lambda=1.5
    epoch 1/10 - loss : 0.851 - val_loss : 0.847
    epoch 2/10 - loss : 0.831 - val_loss : 0.831
    epoch 3/10 - loss : 0.828 - val_loss : 0.829
    epoch 4/10 - loss : 0.827 - val_loss : 0.828
    epoch 5/10 - loss : 0.827 - val_loss : 0.828
    epoch 6/10 - loss : 0.826 - val_loss : 0.828
    epoch 7/10 - loss : 0.826 - val_loss : 0.828
    epoch 8/10 - loss : 0.826 - val_loss : 0.828
    epoch 9/10 - loss : 0.826 - val_loss : 0.828
    epoch 10/10 - loss : 0.826 - val_loss : 0.828



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

    validation error : 0.828





    0.8276982643684648



## MovieLens 1M

### Evaluation on raw data


```python
# load the ml1m dataset
ratings, movies = ml1m.load()

ratings, uencoder, iencoder = ids_encoder(ratings)

m = ratings.userid.nunique()   # total number of users
n = ratings.itemid.nunique()   # total number of items

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

# train test split
(x_train, x_test), (y_train, y_test) = train_test_split(examples=raw_examples, labels=raw_labels)
```


```python
# create the model
MF = MatrixFactorization(m, n, k=10, alpha=0.01, lamb=1.5)

# fit the model on the training set
history = MF.fit(x_train, y_train, epochs=epochs, validation_data=(x_test, y_test))
```

    Training Matrix Factorization Model ...
    k=10 	 alpha=0.01 	 lambda=1.5
    epoch 1/10 - loss : 1.713 - val_loss : 1.718
    epoch 2/10 - loss : 1.523 - val_loss : 1.526
    epoch 3/10 - loss : 1.496 - val_loss : 1.498
    epoch 4/10 - loss : 1.489 - val_loss : 1.489
    epoch 5/10 - loss : 1.485 - val_loss : 1.486
    epoch 6/10 - loss : 1.484 - val_loss : 1.484
    epoch 7/10 - loss : 1.483 - val_loss : 1.483
    epoch 8/10 - loss : 1.483 - val_loss : 1.483
    epoch 9/10 - loss : 1.482 - val_loss : 1.482
    epoch 10/10 - loss : 1.482 - val_loss : 1.482



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

    validation error : 1.482





    1.4820034560467208



### Evaluation on normalized ratings


```python
# load data
ratings, movies = ml1m.load()

ratings, uencoder, iencoder = ids_encoder(ratings)

m = ratings['userid'].nunique()   # total number of users
n = ratings['itemid'].nunique()   # total number of items

# normalize ratings by substracting means
normalized_column_name = "norm_rating"
ratings = normalized_ratings(ratings, norm_column=normalized_column_name)

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

# train test split
(x_train, x_test), (y_train, y_test) = train_test_split(examples=raw_examples, labels=raw_labels)
```


```python
# create the model
MF = MatrixFactorization(m, n, k=10, alpha=0.01, lamb=1.5)

# fit the model on the training set
history = MF.fit(x_train, y_train, epochs=epochs, validation_data=(x_test, y_test))
```

    Training Matrix Factorization Model ...
    k=10 	 alpha=0.01 	 lambda=1.5
    epoch 1/10 - loss : 0.826 - val_loss : 0.827
    epoch 2/10 - loss : 0.824 - val_loss : 0.825
    epoch 3/10 - loss : 0.823 - val_loss : 0.825
    epoch 4/10 - loss : 0.823 - val_loss : 0.825
    epoch 5/10 - loss : 0.823 - val_loss : 0.825
    epoch 6/10 - loss : 0.823 - val_loss : 0.825
    epoch 7/10 - loss : 0.823 - val_loss : 0.825
    epoch 8/10 - loss : 0.823 - val_loss : 0.825
    epoch 9/10 - loss : 0.823 - val_loss : 0.825
    epoch 10/10 - loss : 0.823 - val_loss : 0.825



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

    validation error : 0.825





    0.8250208634455388



## Predictions

Now that the latent factors $P$ and $Q$, we can use them to make predictions and recommendations. Let's call the ```predict``` function of the ```Matrix Factorization``` class to make prediction for a given.

rating prediction for user 1 on item 1 for which the truth rating $r=5.0$


```python
ratings.userid = uencoder.inverse_transform(ratings.userid.to_list())
ratings.itemid = uencoder.inverse_transform(ratings.itemid.to_list())
ratings.head(5)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>userid</th>
      <th>itemid</th>
      <th>rating</th>
      <th>rating_mean</th>
      <th>norm_rating</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1</td>
      <td>1</td>
      <td>5</td>
      <td>4.188679</td>
      <td>0.811321</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1</td>
      <td>48</td>
      <td>5</td>
      <td>4.188679</td>
      <td>0.811321</td>
    </tr>
    <tr>
      <th>2</th>
      <td>1</td>
      <td>145</td>
      <td>5</td>
      <td>4.188679</td>
      <td>0.811321</td>
    </tr>
    <tr>
      <th>3</th>
      <td>1</td>
      <td>254</td>
      <td>4</td>
      <td>4.188679</td>
      <td>-0.188679</td>
    </tr>
    <tr>
      <th>4</th>
      <td>1</td>
      <td>514</td>
      <td>5</td>
      <td>4.188679</td>
      <td>0.811321</td>
    </tr>
  </tbody>
</table>
</div>




```python
4.188679 + MF.predict(userid=1, itemid=1) # add the mean because we have used the normalised ratings for training
```




    4.188679163563357



## Summary

This is the link to the MatrixFactorization class : [MatrixFactorization.py](https://github.com/nzhinusoftcm/review-on-collaborative-filtering/blob/master/recsys/models/MatrixFactorization.py)

## Non-negative Matrix Factorization

With the Matrix Factorization model, $P$ and $Q$ latent factors are non interpretable since they contain arbitrary negative or positive values. This make the Matrix Factorization model to be non explainable.

The **Non-negative Matrix Factorization (NMF)** algorithm is a variant of Matrix Factorization which generate explainable latent factors for $P$ and $Q$ by constraining their values in the range [0,1]. This allow a probability interpretation of these latent factors, hense the explainability.

Click [here](https://github.com/nzhinusoftcm/review-on-collaborative-filtering/blob/master/6.NonNegativeMatrixFactorization.ipynb) to go to the NMF notebook.

## Reference

1. Yehuda Koren et al. (2009). <a href='https://ieeexplore.ieee.org/document/5197422'>Matrix Factorization Techniques for Recommender Systems</a>

## 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.

654e40668304241d2c538e2a7c1c1d9d7cd12091
ipynb
Jupyter Notebook
_notebooks/2021-06-23-collaborative-filtering-movielens-latest-small-05.ipynb
recohut/notebook
610670666a1c3d8ef430d42f712ff72ecdbd8f86
Apache-2.0
_notebooks/2021-06-23-collaborative-filtering-movielens-latest-small-05.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-05.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