.. _movielensfmrst:

===========================================
Factorisation de matrice et recommandations
===========================================


.. only:: html

    **Links:** :download:`notebook <movielens_fm.ipynb>`, :downloadlink:`html <movielens_fm2html.html>`, :download:`PDF <movielens_fm.pdf>`, :download:`python <movielens_fm.py>`, :downloadlink:`slides <movielens_fm.slides.html>`, :githublink:`GitHub|_doc/notebooks/lectures/movielens_fm.ipynb|*`


Le notebook utilise la factorisation de matrice pour calculer des
recommandations sur la base
`movielens <https://grouplens.org/datasets/movielens/>`__. On utilise le
jeu de données
`ml-latest-small.zip <http://files.grouplens.org/datasets/movielens/ml-latest-small.zip>`__.

.. code:: ipython3

    %matplotlib inline

.. code:: ipython3

    from papierstat.datasets import load_movielens_dataset
    data = load_movielens_dataset(cache='movielens.zip')
    list(sorted(data))




.. parsed-literal::
    ['links', 'movies', 'ratings', 'tags']



.. code:: ipython3

    data['movies'].tail(n=2)






.. raw:: html

    <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>movieId</th>
          <th>title</th>
          <th>genres</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>9123</th>
          <td>164977</td>
          <td>The Gay Desperado (1936)</td>
          <td>Comedy</td>
        </tr>
        <tr>
          <th>9124</th>
          <td>164979</td>
          <td>Women of '69, Unboxed</td>
          <td>Documentary</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    import pandas
    rate = data["ratings"]
    rate["dt"] = pandas.to_datetime(rate["timestamp"], unit='s')
    rate.tail(n=2)






.. raw:: html

    <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>movieId</th>
          <th>rating</th>
          <th>timestamp</th>
          <th>dt</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>100002</th>
          <td>671</td>
          <td>6385</td>
          <td>2.5</td>
          <td>1070979663</td>
          <td>2003-12-09 14:21:03</td>
        </tr>
        <tr>
          <th>100003</th>
          <td>671</td>
          <td>6565</td>
          <td>3.5</td>
          <td>1074784724</td>
          <td>2004-01-22 15:18:44</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    len(set(rate['userId'])), len(set(rate['movieId']))




.. parsed-literal::
    (671, 9066)



.. code:: ipython3

    rate.describe()






.. raw:: html

    <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>movieId</th>
          <th>rating</th>
          <th>timestamp</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>count</th>
          <td>100004.000000</td>
          <td>100004.000000</td>
          <td>100004.000000</td>
          <td>1.000040e+05</td>
        </tr>
        <tr>
          <th>mean</th>
          <td>347.011310</td>
          <td>12548.664363</td>
          <td>3.543608</td>
          <td>1.129639e+09</td>
        </tr>
        <tr>
          <th>std</th>
          <td>195.163838</td>
          <td>26369.198969</td>
          <td>1.058064</td>
          <td>1.916858e+08</td>
        </tr>
        <tr>
          <th>min</th>
          <td>1.000000</td>
          <td>1.000000</td>
          <td>0.500000</td>
          <td>7.896520e+08</td>
        </tr>
        <tr>
          <th>25%</th>
          <td>182.000000</td>
          <td>1028.000000</td>
          <td>3.000000</td>
          <td>9.658478e+08</td>
        </tr>
        <tr>
          <th>50%</th>
          <td>367.000000</td>
          <td>2406.500000</td>
          <td>4.000000</td>
          <td>1.110422e+09</td>
        </tr>
        <tr>
          <th>75%</th>
          <td>520.000000</td>
          <td>5418.000000</td>
          <td>4.000000</td>
          <td>1.296192e+09</td>
        </tr>
        <tr>
          <th>max</th>
          <td>671.000000</td>
          <td>163949.000000</td>
          <td>5.000000</td>
          <td>1.476641e+09</td>
        </tr>
      </tbody>
    </table>
    </div>



671 utilisateurs et 9066 films. C’est petit mais assez pour voir la
factorisation et le temps que cela prend. Quelques idées sur les
données.

.. code:: ipython3

    ax = rate['rating'].hist(bins=10, figsize=(3,3))
    ax.set_title('Distribution des ratings');



.. image:: movielens_fm_8_0.png


Les gens préfèrent les ratings arrondis.

.. code:: ipython3

    ax = rate['dt'].hist(bins=50, figsize=(10,3))
    ax.set_title('Distribution des dates');



.. image:: movielens_fm_10_0.png


.. code:: ipython3

    import matplotlib.pyplot as plt
    fig, ax = plt.subplots(1,2, figsize=(6,3))
    gr = rate[["userId", "movieId"]].groupby('userId').count()
    gr.hist('movieId', bins=50, figsize=(3,3), ax=ax[0])
    ax[0].set_yscale('log')
    ax[0].set_title('Distribution du nombre de\nfilms notés par utilisateur')
    gr = rate[["userId", "movieId"]].groupby('movieId').count()
    gr.hist('userId', bins=50, figsize=(3,3), ax=ax[1])
    ax[1].set_yscale('log')
    ax[1].set_title('Distribution du nombre de\nnotes par film');



.. image:: movielens_fm_11_0.png


Il y a quelques utilisateurs zélés et quelques films suscitant beaucoup
d’intérêt. Ce ne sont pas des valeurs aberrantes mais il faudra songer à
regarder ça de plus près un jour. Noter plus de 2000 films paraît
suspect. Même si les votes s’étalent sur les 20 ans de collecte, cela
fait un film tous les 3-4 jours. Il faut transformer les données sous la
forme d’une matrice
`sparse <https://docs.scipy.org/doc/scipy/reference/sparse.html>`__.

.. code:: ipython3

    from scipy.sparse import csr_matrix
    import numpy
    
    
    def conversion(rating, shape=None, movieId_col=None, userId_row=None):
        rating = rating[['userId', 'movieId', 'rating']].dropna()
        coefs = {}
        posix = {}
        movieId_col = movieId_col.copy() if movieId_col is not None else {}
        userId_row = userId_row.copy() if userId_row is not None else {}
        for ind, uid, mid, note in rating.itertuples():
            if uid not in userId_row:
                userId_row[uid] = len(userId_row)
            if mid not in movieId_col:
                movieId_col[mid] = len(movieId_col)
            row = userId_row[uid]
            col = movieId_col[mid]
            if row not in coefs:
                coefs[row] = []
                posix[row] = []
            coefs[row].append(note)
            posix[row].append(col)
            
        nbcoefs = sum(map(len, coefs.values()))
        indptr = numpy.zeros(len(coefs)+1)
        indices = numpy.zeros(nbcoefs)
        data = numpy.zeros(nbcoefs)    
        nb = 0
        for row in range(len(userId_row)):
            cs = coefs[row]
            ps = posix[row]
            indptr[row] = nb
            for i, (p, c) in enumerate(sorted(zip(ps, cs))):
                indices[nb] = p
                data[nb] = c
                nb += 1    
        
        indptr[-1] = nb
        if shape is None:
            shape = (len(userId_row), len(movieId_col))
        mat = csr_matrix((data, indices, indptr), shape=shape)
        if mat.max() != data.max():
            end = min(10, len(indptr))
            raise RuntimeError("La conversion s'est mal passée.\ndata={0}\nindices={1}\nindptr={2}".format(
                            data[:end], indices[:end], indptr[:end]))
        return mat, userId_row, movieId_col
    
    
    petit = pandas.DataFrame(dict(userId=[0, 1, 1, 5, 5], movieId=[0, 1, 2, 4, 10], 
                                  rating=[1, 2, 3, 4, 5]))
    
    mat, userId_row, movieId_col = conversion(petit)
    numpy.nan_to_num(mat.todense())




.. parsed-literal::
    matrix([[1., 0., 0., 0., 0.],
            [0., 2., 3., 0., 0.],
            [0., 0., 0., 4., 5.]])



.. code:: ipython3

    userId_row, '*', movieId_col




.. parsed-literal::
    ({0: 0, 1: 1, 5: 2}, '*', {0: 0, 1: 1, 2: 2, 4: 3, 10: 4})



.. code:: ipython3

    mat, userId_row, movieId_col = conversion(rate)
    numpy.nan_to_num(mat[:5,:5].todense())




.. parsed-literal::
    matrix([[2.5, 3. , 3. , 2. , 4. ],
            [0. , 0. , 0. , 0. , 0. ],
            [0. , 0. , 0. , 0. , 0. ],
            [0. , 0. , 0. , 0. , 0. ],
            [0. , 0. , 0. , 0. , 0. ]])



On cale une factorisation de matrice.

.. code:: ipython3

    from sklearn.decomposition import NMF
    mf = NMF(n_components=20, shuffle=True, max_iter=400)
    mf.fit(mat)




.. parsed-literal::
    NMF(alpha=0.0, beta_loss='frobenius', init=None, l1_ratio=0.0, max_iter=400,
      n_components=20, random_state=None, shuffle=True, solver='cd',
      tol=0.0001, verbose=0)



.. code:: ipython3

    mf.reconstruction_err_ 




.. parsed-literal::
    897.9791684368183



.. code:: ipython3

    wh = mf.transform(mat)
    wh[:5,:5]




.. parsed-literal::
    array([[0.        , 0.05039834, 0.        , 0.00358679, 0.        ],
           [0.        , 0.01736155, 1.91278876, 0.        , 0.        ],
           [0.20592908, 0.        , 0.29535495, 0.        , 0.        ],
           [0.        , 0.11052953, 0.        , 0.39806458, 0.        ],
           [0.71489676, 0.        , 0.20115088, 0.02163221, 0.        ]])



L’erreur ne dit pas grand chose sur la pertinence de la recommandation.
Le plus simple est d’enlever des notes pour voir si on les retrouve.

.. code:: ipython3

    from sklearn.model_selection import train_test_split
    rate_train, rate_test = train_test_split(rate)

Il faut quand même s’assurer que la matrice à décomposer a les mêmes
dimensions que la précédente avec toutes les données.

.. code:: ipython3

    shape0 = mat.shape
    mat_train, userId_row_train, movieId_col_train = conversion(rate, 
                        shape=shape0, userId_row=userId_row, movieId_col=movieId_col)

.. code:: ipython3

    mf.fit(mat_train)




.. parsed-literal::
    NMF(alpha=0.0, beta_loss='frobenius', init=None, l1_ratio=0.0, max_iter=400,
      n_components=20, random_state=None, shuffle=True, solver='cd',
      tol=0.0001, verbose=0)



.. code:: ipython3

    mf.reconstruction_err_




.. parsed-literal::
    898.1781492558509



On calcule l’erreur sur les bases d’apprentissage et de test.

.. code:: ipython3

    def predict(mf, mat_train, test, userId_row, movieId_col):    
        W = mf.transform(mat_train)
        H = mf.components_
        wh = W @ H
        test = test[['userId', 'movieId', 'rating']]
        predictions = []
        for ind, uid, mid, note in test.itertuples(): 
            row = userId_row[uid]
            col = movieId_col[mid]
            try:
                pred = wh[row, col]        
            except Exception as e:
                raise Exception("Issue with uid={} mid={} row={} col={} shape={}".format(uid, mid, row, col, wh.shape))
            predictions.append((ind, pred))
        dfpred = pandas.DataFrame(data=predictions, columns=['index', 'prediction']).set_index('index')
        dfall = pandas.concat([test, dfpred], axis=1)
        return dfall
    
    pred = predict(mf, mat_train, rate_test, userId_row_train, movieId_col_train)
    pred.head()






.. raw:: html

    <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>movieId</th>
          <th>rating</th>
          <th>prediction</th>
        </tr>
        <tr>
          <th>index</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>13240</th>
          <td>85</td>
          <td>356</td>
          <td>4.0</td>
          <td>4.934174</td>
        </tr>
        <tr>
          <th>75789</th>
          <td>527</td>
          <td>1276</td>
          <td>4.5</td>
          <td>0.345162</td>
        </tr>
        <tr>
          <th>42713</th>
          <td>306</td>
          <td>1513</td>
          <td>4.0</td>
          <td>2.439180</td>
        </tr>
        <tr>
          <th>56562</th>
          <td>407</td>
          <td>3635</td>
          <td>3.0</td>
          <td>0.406101</td>
        </tr>
        <tr>
          <th>63681</th>
          <td>457</td>
          <td>122882</td>
          <td>0.5</td>
          <td>0.903360</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    from sklearn.metrics import r2_score
    r2_score(pred['rating'], pred['prediction'])




.. parsed-literal::
    -4.659895568960519



Pas extraordinaire. Faisons varier *k*.

.. code:: ipython3

    from time import perf_counter as clock
    from sklearn.metrics import mean_squared_error
    values = []
    
    for k in [5, 10, 15, 20, 25, 30, 35]:
        mem = {}
        mf = NMF(n_components=k, shuffle=True, max_iter=400)
        cl = clock()
        mf.fit(mat_train)
        mem['train_time'] = clock() - cl
        pred = predict(mf, mat_train, rate_test, userId_row_train, movieId_col_train)
        mem['k'] = k
        mem['r2'] = r2_score(pred['rating'], pred['prediction'])
        mem['err_test'] = mean_squared_error(pred['rating'], pred['prediction'])
        mem['err_train'] = mf.reconstruction_err_
        values.append(mem)
        print(k, mem)    


.. parsed-literal::
    5 {'train_time': 0.7479112074008754, 'k': 5, 'r2': -6.014138361502092, 'err_test': 7.929004708705666, 'err_train': 981.0205349504541}
    10 {'train_time': 1.910340634477734, 'k': 10, 'r2': -5.388353487045178, 'err_test': 7.221597617417111, 'err_train': 943.8568240851894}
    15 {'train_time': 2.2207374467998306, 'k': 15, 'r2': -4.96900640010828, 'err_test': 6.747554355716667, 'err_train': 918.0198374341521}
    20 {'train_time': 5.637187125555101, 'k': 20, 'r2': -4.694288211755458, 'err_test': 6.437004193066304, 'err_train': 897.9357561628665}
    25 {'train_time': 7.713239363839193, 'k': 25, 'r2': -4.420273628728207, 'err_test': 6.127249408216037, 'err_train': 878.84542031377}
    30 {'train_time': 12.43995074364534, 'k': 30, 'r2': -4.195368644607753, 'err_test': 5.873009673241581, 'err_train': 862.2473126443812}
    35 {'train_time': 15.610665020047463, 'k': 35, 'r2': -3.9997376232229183, 'err_test': 5.6518621552167065, 'err_train': 846.2449517351943}


.. code:: ipython3

    df = pandas.DataFrame(values)
    df






.. raw:: html

    <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>err_test</th>
          <th>err_train</th>
          <th>k</th>
          <th>r2</th>
          <th>train_time</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>7.929005</td>
          <td>981.020535</td>
          <td>5</td>
          <td>-6.014138</td>
          <td>0.747911</td>
        </tr>
        <tr>
          <th>1</th>
          <td>7.221598</td>
          <td>943.856824</td>
          <td>10</td>
          <td>-5.388353</td>
          <td>1.910341</td>
        </tr>
        <tr>
          <th>2</th>
          <td>6.747554</td>
          <td>918.019837</td>
          <td>15</td>
          <td>-4.969006</td>
          <td>2.220737</td>
        </tr>
        <tr>
          <th>3</th>
          <td>6.437004</td>
          <td>897.935756</td>
          <td>20</td>
          <td>-4.694288</td>
          <td>5.637187</td>
        </tr>
        <tr>
          <th>4</th>
          <td>6.127249</td>
          <td>878.845420</td>
          <td>25</td>
          <td>-4.420274</td>
          <td>7.713239</td>
        </tr>
        <tr>
          <th>5</th>
          <td>5.873010</td>
          <td>862.247313</td>
          <td>30</td>
          <td>-4.195369</td>
          <td>12.439951</td>
        </tr>
        <tr>
          <th>6</th>
          <td>5.651862</td>
          <td>846.244952</td>
          <td>35</td>
          <td>-3.999738</td>
          <td>15.610665</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    fig, ax = plt.subplots(1, 4, figsize=(12, 3))
    df.plot(x='k', y="r2", style='o-', ax=ax[0])
    ax[0].set_title("NMF\nr2 base de test\net k")
    df.plot(x='k', y="err_test", style='o-', ax=ax[1])
    ax[1].set_title("NMF\nerreur de test\net k");
    df.plot(x='k', y="err_train", style='o-', ax=ax[2])
    ax[2].set_title("NMF\nerreur d'apprentissage\net k")
    df.plot(y='train_time', x="k", style='o-', ax=ax[3])
    ax[3].set_title("NMF\nk\net temps d'apprentissage");



.. image:: movielens_fm_32_0.png


Il faudrait explorer de plus grandes valeurs de *k*, il faudrait aussi
faire de la cross-validation puis regarder si on peut corrélérer les
plus autres erreurs à certains type d’utilisateurs ou de films, si on
arrive à déterminer s’ils se distingue des autres par un faible ou fort
taux de *ratings*, moyenne, plus ou moins proches des utilisateurs
typiques (~*H*) ou des films typiques (~*W*). Bref, ce n’est pas fini.