<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Background2&nbsp;&nbsp;Setup3&nbsp;&nbsp;Data3.1&nbsp;&nbsp;About the data3.2&nbsp;&nbsp;Read data3.3&nbsp;&nbsp;Matrix form4&nbsp;&nbsp;Aux functions5&nbsp;&nbsp;User recommendations5.1&nbsp;&nbsp;Product similarity based recommendations5.1.1&nbsp;&nbsp;Distances5.1.2&nbsp;&nbsp;Rules for similarity based recommendations5.1.3&nbsp;&nbsp;Example 15.2&nbsp;&nbsp;User similarity based recommendations6&nbsp;&nbsp;Collaborative Filtering6.1&nbsp;&nbsp;Parameters6.2&nbsp;&nbsp;Simple Gradient Descent6.2.1&nbsp;&nbsp;Initialize6.2.2&nbsp;&nbsp;Gradient Descent6.3&nbsp;&nbsp;Regularisation6.4&nbsp;&nbsp;Stochastic GD6.5&nbsp;&nbsp;User/product clusters6.6&nbsp;&nbsp;Calculating recommendations7&nbsp;&nbsp;Finally, a note on scalability  BackgroundThis is based on ideas from this [post](http://www.albertauyeung.com/post/python-matrix-factorization/), but takes a more simple, script-like approach. Our motivation is to try and recommend _new_ products to our users, based on their rating of products they already engaged with.  Setup Data & Algebra libs<jupyter_code>import pandas as pd
import numpy as np
from sklearn.metrics import pairwise_distances<jupyter_output><empty_output><jupyter_text>Graphic libs<jupyter_code>from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
from scipy.cluster import hierarchy<jupyter_output><empty_output><jupyter_text>Data  About the dataYou can find everything you need to know about the MovieLens dataset [here](http://files.grouplens.org/datasets/movielens/ml-100k/README)  Read data<jupyter_code>ml = \
pd.read_csv(
    filepath_or_buffer = 'http://files.grouplens.org/datasets/movielens/ml-100k/u.data',
    header = None,
    sep = '\t'
)
ml.columns = 'user id | item id | rating | timestamp'.split(' | ')<jupyter_output><empty_output><jupyter_text>For off-line usage<jupyter_code># ml.to_feather('movie_rec.feather')
# ml = pd.read_feather('movie_rec.feather')<jupyter_output><empty_output><jupyter_text>Our dataset at a glance<jupyter_code>ml.shape
ml.head()
ml.describe()<jupyter_output><empty_output><jupyter_text>Matrix form<jupyter_code>data_mat = \
ml.drop(columns=['timestamp']).\
    pivot(
        index = 'user id', 
        columns = 'item id', 
        values = 'rating'
    ).\
    fillna(0)<jupyter_output><empty_output><jupyter_text>Aux functions<jupyter_code>def vis_dist_mat(dist, bins = 30):
    '''Visualize different aspects of a distance matrix
    
    :param dist: a distance matrix calculated by `pairwise_distances`
    
    :Example:
    >>> vis_dist_mat(dist)
    '''
    
    dist_vec = \
    pd.Series(
        dist[
            np.tril_indices(dist.shape[0], k=1)
        ]
    )
    
    print('Distance Matrix Details\n=======================\n')
    
    print(
        '* have have',
        (dist_vec == 0).sum(),
        'identical items (with distance = 0), and',
        (dist_vec == 1).sum(),
        'with no common ratings (distance = 1)\n'
    )
    print('* Percentiles:')
    print(dist_vec.describe(), '\n')
    
    sns.distplot(
        dist_vec,
        bins = bins,
        kde = False,
        norm_hist = False,
    ).\
    set_title('Distance distribution')
    
    plt.figure()

    dn = hierarchy.dendrogram(hierarchy.ward(dist))
    plt.xticks([])
    plt.title('Cluster dendogram')

    plt.show()
    
    return None<jupyter_output><empty_output><jupyter_text>User recommendations  Product similarity based recommendationsThis approach looks at all the products recommended by the users and try to find "similar" products based on some simple rules.  Distances<jupyter_code>product_dist = pairwise_distances(X = data_mat.T, metric = 'cosine')
vis_dist_mat(product_dist)<jupyter_output>Distance Matrix Details
=======================

* have have 2392 identical items (with distance = 0), and 430641 with no common ratings (distance = 1)

* Percentiles:
count    1.417084e+06
mean     9.129765e-01
std      1.089063e-01
min      0.000000e+00
25%      8.696132e-01
50%      9.451297e-01
75%      1.000000e+00
max      1.000000e+00
dtype: float64<jupyter_text>It's interesting that a lot of products have close to no similarity at all (distance is close to 1 on the histogram). This is probably a reflection of the almost trivial 2 clusters we can see in the dendogram. It's possible that more data on the movies in each cluster will reveal a straight forward way to split our recommendation problem into 2 separate, smaller problems:  if, for example, the clusters here turn out to be something like "drama" vs. "comedy" movies clusters this would suggest that there are separate groups of users that typically rate one genre (either "comedy" or "drama") and that only a few users rate both. We can then think of building separate recommendation systems based on this fact. For the purpose of this demonstration we keep a single system, but such insights reflect the importance of looking at your data...  Rules for similarity based recommendations There are several strategies for picking the products we want to recommend based on product similarity, but essentially we need to make several choices1. What products we look at: all the products the user rated? Only products with $n$ stars or more? Or do we weights the distances by the number of stars? And even with weights - do 1 star reviews have a positive weight or should we give a low number of stars a negative weight? 2. How do we aggregate the distances? Say the user only rated 2 products, and both of them are similar to product $j$. Do we rank product $j$ compared to all the other products that have some similarity to the 2 rated products? Do we sum up distances for each product and rank? Do we take $min$/$max$/... distance for each product?3. What would be the threshold for a product to be recommended? This is dependent on practical limitations, e.g. what's the minimum/maximum number of products we have to/can recommend, and what is the distribution of distances and other considerations (e.g. do we want an absolute threshold of similarity and not just the top 5) Looking at the second choice we have to make: there's one huge advantage to summing up distances over picking min/max or other aggregation methods - __we can use matrix multiplication!__  Example 1We decide to look at all the movies a user rated (no matter what the rating was) and sum up the similarities to other products, and recommend the top 5 most similar products (regardless of the absolute level of similarity - in this example we _have_ to recommend exactly 5 products, no matter what) We define 1. The symmetric matrix $P$ containing product to product distances as we calculated above (our `product_dist` matrix based on `cosine` distance)2. A matrix with binary values $C$ \begin{equation*}    C_{I \times J} = \left[         \begin{array}{lll}            c_{1,1} & \ldots & c_{1,J} \\            \vdots  & \ddots & \vdots  \\            c_{J,1} & \ldots & c_{J,J}        \end{array}    \right]\end{equation*}where\begin{equation*}    c_{i,j} =     \left\{        \begin{array}{ll}          1 & R_{i,j} > 0 \\          0 & \text{otherwise}        \end{array}    \right.\end{equation*} For every user $i$ we want to look at all rated products, and find the distances from other products. This can be written as a sparse matrix, where only rows matching rated products will have non-zero values: \begin{equation*}   M^{(i)}_{J \times J} =  \left[        \begin{array}{lll}            c_{i,1} p_{1,1} & \ldots & c_{i,1} p_{1,J} \\            \vdots  & \ddots & \vdots  \\            c_{i,J} p_{J,1} & \ldots & c_{i,J} p_{J,J}        \end{array}    \right]\end{equation*} And if our way to rank products is to sum up distances, we can sum up across columns, to get $\left( \sum_{k =1}^{J} c_{i,k} p_{k,1} , \ldots, \sum_{k=1}^{J}  c_{i,k} p_{k,J} \right)$, which is simply the $i$'th rows of the matrix $CP$.  Naturally, some of the products we see in this list with non-zero weight have already been encountered by the user (rated, viewed or purchased, depending on the context), so we would not want to recommend them at all. We can eliminate these products [element-by-element product](https://en.wikipedia.org/wiki/Hadamard_product) with the matrix $(1 - C)$, so our final matrix is: \begin{equation*}   CP \circ \left( 1 - C \right)\end{equation*} which we calculate as:<jupyter_code>CP = (data_mat > 0).dot(product_dist) 
CP.columns = data_mat.columns
user_prod_weights = CP * (data_mat == 0)<jupyter_output><empty_output><jupyter_text>For practical reasons we need to work with the long format, and we can filter all zero weights:<jupyter_code>user_rec = pd.DataFrame({'dist': user_prod_weights.T.unstack()})
user_rec = user_rec.loc[user_rec['dist'] > 0, ]
print(
    'long format has', 
    user_rec.shape[0], 
    'rows, which is', 
    round(100.0 * user_rec.shape[0] / (data_mat.shape[0] * data_mat.shape[1]), 2), 
    '% of all possible',
    data_mat.shape[0] * data_mat.shape[1],
    'combinations'
)<jupyter_output>long format has 1486126 rows, which is 93.7 % of all possible 1586126 combinations<jupyter_text>This did not save much space (our distance matrix is not very sparse and there's a lot of overlap in user preference) but at least we easily can pull out only the top 5 most similar products (or apply any other selection rule). These can be exported back to our production system (as CSV, to the production DB or wrapped with a micro-service)<jupyter_code>user_rec.\
groupby(user_rec.index.get_level_values(0)).\
apply(
    lambda x: x.sort_values(by = 'dist').head(5).reset_index()
).head(15)<jupyter_output><empty_output><jupyter_text>User similarity based recommendations<jupyter_code>user_dist = pairwise_distances(X = data_mat, metric = 'cosine')
vis_dist_mat(user_dist)<jupyter_output>Distance Matrix Details
=======================

* have have 943 identical items (with distance = 0), and 15075 with no common ratings (distance = 1)

* Percentiles:
count    446038.000000
mean          0.826180
std           0.124427
min           0.000000
25%           0.753425
50%           0.848349
75%           0.920774
max           1.000000
dtype: float64<jupyter_text>As a reflection of our finding on products we see again 2 major clusters, and again we can consider splitting the problem in 2. In any case the next steps are to find the the top $n$ most similar users (in general or within their cluster) and select the products these users liked most.  Collaborative FilteringSince the above solution does not scale up well with the number of users $I$ or the number of products $J$, we may want to explore dimensionality reduction methods. In this case we try to find 2 lower dimension matrices $P_{I,K}, Q_{K,J}$ such that our original user X product matrix $R_{I,J}$ can be approximated by $PQ$. By setting $K$ to be relatively low we can tightly control the dimensionality of the problem.  Parameters<jupyter_code>K = 2
iterations = 1000
tolerance = 0.05
alpha = 0.0001<jupyter_output><empty_output><jupyter_text>Simple Gradient Descent  Initialize<jupyter_code>R = np.array(data_mat)
P = np.random.normal(scale=1./K, size=(R.shape[0], K))
Q = np.random.normal(scale=1./K, size=(K, R.shape[1]))
training_process = pd.DataFrame(columns=['iteration', 'mse'])
t = 1
mse_prev = 0
mse = 1e6<jupyter_output><empty_output><jupyter_text>Gradient Descent<jupyter_code>while t <= iterations and abs(mse - mse_prev) > tolerance:
    mse_prev = mse
    
    ## Error matrix
    error = R - P.dot(Q)

    ## Update down the gradient
    P = P + 2 * alpha * error.dot(Q.T)
    Q = Q + 2 * alpha * P.T.dot(error)

    mse = np.power(error, 2).sum()
    training_process = training_process.append(pd.DataFrame({'iteration': [t], 'mse': [mse]}))

    if (t+1) % 50 == 0:
        print("Iteration: %d ; error = %.4f" % (t+1, mse))
    t += 1
training_process.set_index('iteration', inplace=True)
plt.figure()
plt.plot(training_process['mse'])
plt.title('Algorithm progress')<jupyter_output><empty_output><jupyter_text>Regularisation TBD* Add a penalty on the values of $P,Q$ to avoid extreme solutions and over-fitting. We can use the $L2$ norm $\| Q \|_{2} +\| P \|_{2} = \sqrt{\sum q_{i,j}^2} + \sqrt{\sum p_{i,j}^2}$) or the $L1$ norm $\| Q \|_{1} +\| P \|_{1} = \sum \left| q_{i,j} \right| + \sum \left| p_{i,j} \right| $ as a penalty on our MSE goal function by adding a penalty coeff $\gamma$ on the MSE goal function and optimizing $MSE - \gamma (\|P\| + \|Q\|)$ * This probably requires some form of normalization of $P,Q$ prior to optimization   Stochastic GD TBD  User/product clusters Since we picked $K=2$ we can actually visualize our results<jupyter_code>plt.figure()
plt.scatter(x = Q[0:1], y = Q[1:])
plt.title('Product dimensions')
plt.figure()
plt.scatter(x = P.T[0:1], y = P.T[1:])
plt.title('User dimensions')<jupyter_output><empty_output><jupyter_text>And of course the fact we do not see a good separation as we've seen before suggests we need to increase $K$. For higher dimensions we can use both clustering methods or more dimension reduction techniques (e.g. t-SNE) to try and find clusters of users.   Calculating recommendations We can now use the embedding matrices $P,Q$ to generate recommendations in a similar way to what we did before, except similarity of product or user is now defined by the Euclidian distances defined by $Q$ or $P$ respectively:<jupyter_code>product_dist = pairwise_distances(Q.T, metric='euclidean')
user_dist = pairwise_distances(P, metric='euclidean')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11036715473475094, 'fraction_numerical': 0.02567875800161872, 'mean_word_length': 4.123256690539012, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 73, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '434815', 'license': 'Apache-2.0', 'n_tokens_mistral': 4413, 'n_tokens_neox': 4125, 'n_words': 1522, 'path': 'matrix-factorization/recommendation_strategies.ipynb', 'repo_name': 'ytoren/reproduce', 'size': '13591'}