Today, we look at a topic that has become enormously important: recommender systems.
In this lecture, we will:
The optional module Recommender Systems II carries the latent factors into deep models (DLRM).
This section draws heavily on
The concept of recommender systems emerged in the late 1990s / early 2000s as social life moved online:
In these systems the amount of content was exploding and users were having a hard time finding things they were interested in.
Users wanted recommendations.
Over time, the problem has only gotten worse:
An enormous need has emerged for systems to help sort through products, services, and content items.
This often goes by the term personalization.
Some examples (rough, recent catalogue sizes):
A more formal view:
Unfortunately, users generally have a hard time explaining what types of content they prefer.
Some early systems worked by interviewing users to ask what they liked. Those systems did not work very well.
A very interesting article about the earliest personalization systems is User Modeling via Stereotypes by Elaine Rich, dating from 1979.
Instead, modern systems work by capturing user’s opinions about specific items.
This can be done actively:
or it can be done passively:
This data is a subset of the data presented in: “From amateurs to connoisseurs: modeling the evolution of user expertise through online reviews,” by J. McAuley and J. Leskovec. WWW, 2013
Let’s look at a dataset for testing recommender systems consisting of Amazon movie reviews.
The full set is 1.7 million reviews with their text (647 MB). We only need two small files derived from it, both in data/amazon-movies/:
review-counts.csv — how many reviews each user wrote and each movie received, stored as a frequency table (kind, reviews, count: “count users wrote exactly reviews reviews”);dense-block.csv — the (user, movie, star rating) triples we will factor later.import os
DATA = 'data/amazon-movies'
if not os.path.exists(DATA): # e.g. on Colab: read the public copy instead
DATA = 'https://raw.githubusercontent.com/tools4ds/DS701-Materials-FA26/main/data/amazon-movies'
counts = pd.read_csv(f'{DATA}/review-counts.csv')
counts.head()| kind | reviews | count | |
|---|---|---|---|
| 0 | item | 5 | 5751 |
| 1 | item | 6 | 4545 |
| 2 | item | 7 | 3546 |
| 3 | item | 8 | 2915 |
| 4 | item | 9 | 2520 |
Expanding the frequency table gives one number per user and one per movie:
def expand(kind):
rows = counts[counts['kind'] == kind]
return np.repeat(rows['reviews'].values, rows['count'].values)
reviews_per_user = expand('user')
reviews_per_movie = expand('item')
n_users = len(reviews_per_user)
n_movies = len(reviews_per_movie)
n_reviews = reviews_per_user.sum()
display(Markdown(f'There are:\n\n* {n_reviews:,} reviews\n* {n_movies:,} movies\n* {n_users:,} users'))
display(Markdown(f'There are {n_users * n_movies:,} potential reviews, meaning sparsity of {(n_reviews/(n_users * n_movies)):0.4%}'))There are:
There are 6,204,445,920 potential reviews, meaning sparsity of 0.0274%
where
\[ \text{sparsity} = \frac{\text{\# of reviews}}{\text{\# of users} \times \text{\# of movies}} = \frac{\text{\# of reviews}}{\text{\# of potential reviews}} \]
Only about 0.03% of the possible reviews exist – 99.97% of the entries are missing.
display(Markdown(f'There are on average {n_reviews/n_movies:0.1f} reviews per movie' +
f' and {n_reviews/n_users:0.1f} reviews per user'))There are on average 33.9 reviews per movie and 13.7 reviews per user
Although on average a movie receives 34 reviews and a user writes 14, almost all movies and users have fewer. On log–log axes the users’ rank–frequency curve is close to a straight line, the signature of a power law; the movie curve is heavy-tailed too. (The dataset keeps only users and movies with at least 5 reviews.)
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
for ax, counts_k, noun in [(axes[0], reviews_per_movie, 'Movie'),
(axes[1], reviews_per_user, 'User')]:
mean = counts_k.mean()
frac_below = np.mean(counts_k < mean)
ranked = np.sort(counts_k)[::-1]
ax.loglog(np.arange(1, len(ranked) + 1), ranked, color=OKABE_ITO[5], lw=2,
label='Number of reviews')
ax.axhline(mean, color=OKABE_ITO[3], lw=2, ls='--', label=f'Average ({mean:.1f})')
ax.set_xlabel(f'{noun} rank (most reviewed first)', fontsize=12)
ax.set_ylabel('Number of reviews', fontsize=12)
ax.set_title(f'Reviews per {noun.lower()}\n{frac_below:0.0%} of {noun.lower()}s below average',
fontsize=13)
ax.legend(fontsize=11)
plt.tight_layout()
Ultimately, our goal is to predict the rating that a user would give to an item.
For that, we need to define a loss or objective function.
A typical objective function is root mean square error (RMSE)
\[ \text{RMSE} = \sqrt{\frac{1}{|S|} \sum_{(u,i)\in S} (\hat{r}_{ui} - r_{ui})^2}, \]
where
The simplest predictor is the global mean, \(\hat r_{ui} = \mu\) for everyone and everything. Every method below must beat its RMSE; we will measure it on held-out ratings at the end.
OK, now we know the problem and the data available. How can we address the problem?
The earliest method developed is called collaborative filtering.
The central idea of collaborative filtering is that the set of known recommendations can be considered to be a bipartite graph.
The nodes of the bipartite graph are users (\(U\)) and items (\(V\)).
Each edge corresponds to a known rating \(r_{ui}.\)
Then recommendations are formed by traversing or processing the bipartite graph.
Joe’s three war films were also watched by two other users. The films they rated highly, ranked by how many of them agree, are Joe’s recommendations.
There are at least two ways this graph can be used.
Two ways to form a rating for the pair \((u, i)\):
Rule of thumb: compare along the short side of the matrix.
Why the short side wins:
For item-item similarity, we’ll look at item-item Collaborative Filtering (CF).
The questions are:
Here is another view of the ratings graph, this time as a matrix that includes missing entries:
Let’s say we want to predict the value of this unknown rating, user 5’s rating of item 1:
We compare item 1 only with items that user 5 has rated: items 3, 4, 5 and 6.
Similarity \(s_{1j}\) = Pearson correlation of the two item rows, computed over the users who rated both (defined on the next slides):
| item \(j\) | co-raters | \(s_{1j}\) |
|---|---|---|
| 3 | users 1, 9, 11 | 0.84 |
| 4 | users 3, 11 | −1.0 |
| 5 | users 3, 6, 11 | −0.87 |
| 6 | users 1, 3, 11 | 1.0 |
We use \(k\)-nearest neighbors with \(k = 2\): the two most similar, items 3 and 6.
These similarity scores tell us how much weight to put on the rating of the other items.
User 5 rated item 3 a 2 and item 6 a 3, so
\[ \hat{r}_{u=5,\,i=1} = \frac{s_{13}\, r_{5,3} + s_{16}\, r_{5,6}}{s_{13} + s_{16}} = \frac{0.84 \cdot 2 + 1.0 \cdot 3}{0.84 + 1.0} = 2.54. \]
The user–user mirror on the same matrix: the users who rated item 1 and share at least three items with user 5 (with two, a correlation is always \(\pm 1\), as for item 4 above) are users 3 and 11. User 3 correlates positively (\(0.87\) over items 4, 5, 6) and rated item 1 a 3; user 11 correlates negatively (\(-0.95\)) and is dropped. User–user CF predicts 3.
How should we assess similarity of items?
A reasonable approach is to consider items similar if their ratings are correlated, i.e.,
\[\rho(X,Y) = \frac{E\left[(X-\mu_X)(Y-\mu_Y)\right]}{\sigma_X \sigma_Y}.\]
However, note that two items will not have ratings in the same positions.
So we want to compute correlation only over the users who rated both the items.
Let’s put the ratings in python lists:
ratings_item_i = [1, np.nan, np.nan, 5, 5, 3, np.nan, np.nan, np.nan, 4, 2, np.nan, np.nan, np.nan, np.nan, 4, np.nan, 5, 4, 1, np.nan]
ratings_item_j = [np.nan, np.nan, 4, 2, 5, np.nan, np.nan, 1, 2, 5, np.nan, np.nan, 2, np.nan, np.nan, 3, np.nan, np.nan, np.nan, 5, 4]
display(Markdown(f'Ratings for item $i$:\n\n{ratings_item_i}'))
display(Markdown(f'Ratings for item $j$:\n\n{ratings_item_j}'))Ratings for item \(i\):
[1, nan, nan, 5, 5, 3, nan, nan, nan, 4, 2, nan, nan, nan, nan, 4, nan, 5, 4, 1, nan]
Ratings for item \(j\):
[nan, nan, 4, 2, 5, nan, nan, 1, 2, 5, nan, nan, 2, nan, nan, 3, nan, nan, nan, 5, 4]
Let’s drop the non-common ratings:
# Keep only the positions where both items have a rating
filtered_ratings_item_i = [rating_i for rating_i, rating_j in zip(ratings_item_i, ratings_item_j) if not np.isnan(rating_i) and not np.isnan(rating_j)]
filtered_ratings_item_j = [rating_j for rating_i, rating_j in zip(ratings_item_i, ratings_item_j) if not np.isnan(rating_i) and not np.isnan(rating_j)]
display(Markdown(f'Common ratings for item $i$: {filtered_ratings_item_i}'))
display(Markdown(f'Common ratings for item $j$: {filtered_ratings_item_j}'))Common ratings for item \(i\): [5, 5, 4, 4, 1]
Common ratings for item \(j\): [2, 5, 5, 3, 5]
Now we can compute the Pearson correlation coefficient:
rho = np.corrcoef(filtered_ratings_item_i, filtered_ratings_item_j)[0,1]
display(Markdown(f'Pearson correlation coefficient: {rho:0.2f}'))Pearson correlation coefficient: -0.43
Which is a moderate negative correlation, meaning that as item \(i\) gets rated higher, item \(j\) gets rated lower.
The same computation reproduces the similarity table of the item–item example.
nan = np.nan
R_toy = np.array([ # rows = items 1..6, columns = users 1..12
[1, nan, 3, nan, nan, 5, nan, nan, 5, nan, 4, nan],
[nan, nan, 5, 4, nan, nan, 4, nan, nan, 2, 1, 3],
[2, 4, nan, 1, 2, nan, 3, nan, 4, 3, 5, nan],
[nan, 2, 4, nan, 5, nan, nan, 4, nan, nan, 2, nan],
[nan, nan, 4, 3, 4, 2, nan, nan, nan, nan, 2, 5],
[1, nan, 3, nan, 3, nan, nan, 2, nan, nan, 4, nan]])
user = 4 # user 5, zero-indexed
for j in np.where(~np.isnan(R_toy[:, user]))[0]:
both = ~np.isnan(R_toy[0]) & ~np.isnan(R_toy[j])
s = np.corrcoef(R_toy[0, both], R_toy[j, both])[0, 1]
raters = ', '.join(str(k + 1) for k in np.where(both)[0])
print(f'item {j + 1}: co-raters {raters}, s_1{j + 1} = {s:5.2f}')item 3: co-raters 1, 9, 11, s_13 = 0.84
item 4: co-raters 3, 11, s_14 = -1.00
item 5: co-raters 3, 6, 11, s_15 = -0.87
item 6: co-raters 1, 3, 11, s_16 = 1.00The Pearson correlation coefficient, often denoted as \(r\), is a measure of the linear correlation between two variables \(X\) and \(Y\). It quantifies the degree to which a linear relationship exists between the variables. The value of \(r\) ranges from -1 to 1, where:
The formula for the Pearson correlation coefficient is:
\[ r = \frac{\sum (X_i - \bar{X})(Y_i - \bar{Y})}{\sqrt{\sum (X_i - \bar{X})^2 \sum (Y_i - \bar{Y})^2}} \]
Where: - \(X_i\) and \(Y_i\) are the individual sample points, - \(\bar{X}\) and \(\bar{Y}\) are the means of the \(X\) and \(Y\) samples, respectively.
The Pearson correlation coefficient is sensitive to outliers and assumes that the relationship between the variables is linear and that the data is normally distributed.
In some cases we will need to work with binary \(r_{ui}\).
For example, purchase histories on an e-commerce site, or clicks on an ad.
In this case, an appropriate replacement for Pearson \(r\) is the Jaccard similarity coefficient or Intersection over Union.
\[ J_{Sim}(\mathbf{x}, \mathbf{y}) = \frac{|\mathbf{x} \cap \mathbf{y}|}{|\mathbf{x} \cup \mathbf{y}|}. \]
See the lecture on similarity measures.
One problem with the story so far arises due to bias.
These properties interfere with similarity assessment.
Bias correction is crucial for CF recommender systems.
We need to include
Hence we need to form a per-item bias of:
\[ b_{ui} = \mu + \alpha_u + \beta_i \]
where
If we gather all these elements together we can form:
How can we estimate the parameters \(\boldsymbol{\alpha}\), \(\boldsymbol{\beta}\), and \(\mu\)?
Let’s assume for a minute that we had a fully-dense matrix of ratings \(R\).
Recall that each of the \(m\) rows of \(R\) represents an item and each of the \(n\) columns represents a user.
One way to do this is to minimize the squared error between the ratings and the biases:
\[ \min_{\boldsymbol{\alpha},\boldsymbol{\beta},\mu} \Vert R - (\mathbf{1}\boldsymbol{\alpha}^T + \boldsymbol{\beta}\mathbf{1}^T + \mu1)\Vert^2 + \lambda(\Vert\boldsymbol{\alpha}\Vert^2 + \Vert\boldsymbol{\beta}\Vert^2). \]
and include a regularization term to minimize the magnitude of the biases.
Here, bold-faced \(\mathbf{1}\) represents appropriately sized vectors of ones, and non-boldfaced \(1\) is an \(m\times n\) matrix of ones.
So \(\mathbf{1}\boldsymbol{\alpha}^T\) is an \(m\times n\) matrix in which each column is constant, equal to that user’s offset \(\alpha_u\).
Similarly, \(\boldsymbol{\beta}\mathbf{1}^T\) is an \(m\times n\) matrix in which each row is constant, equal to that item’s offset \(\beta_i\).
And \(\mu1\) is an \(m\times n\) matrix where each element is \(\mu\).
There is a simple strategy for solving it.
Assume we hold \(\beta\mathbf{1}^T\) and \(\mu1\) constant.
Then the remaining problem is
\[ \min_{\alpha} \Vert R - \mathbf{1}\alpha^T \Vert^2 + \lambda \Vert\alpha\Vert^2, \]
which (for each column of \(R\)) is a standard regularized least squares problem solved via Ridge regression: \(\hat{\boldsymbol{\beta}} = (X^TX + cI)^{-1}X^T\mathbf{y}\).
This problem is convex — in fact jointly convex in \(\mu\), \(\alpha\) and \(\beta\), because the prediction is linear in all three: it is one big least-squares problem. Solving for one block at a time is just convenient. (The matrix-factorization problem below is not jointly convex.)
The strategy for solving is:
Each of the three steps will reduce the overall error. As a result, we iterate over them until convergence.
The last issue is that the matrix \(R\) is not dense - in reality we only have a small subset of its entries.
We simply need to adapt the least-squares solution to only consider the entries in \(R\) that we know.
As a result, the actual calculation is as follows…
Step 1:
\[ \mu = \frac{\sum_{(u, i) \in R} (r_{ui} - \alpha_u - \beta_i)}{|R|} \]
Step 2:
\[ \alpha_u = \frac{\sum_{i \in R(u)}(r_{ui} - \mu - \beta_i)}{\lambda + |R(u)|} \]
Step 3:
\[ \beta_i = \frac{\sum_{u \in R(i)}(r_{ui} - \mu - \alpha_u)}{\lambda + |R(i)|} \]
Step 4: If not converged, go to Step 1.
Here \(i \in R(u)\) means the set of items rated by user \(u\) and \(u \in R(i)\) means the set of users who have rated item \(i\) and \(|R(u)|\) is the number of ratings.
Now that we have learned the biases, we can do a better job of estimating correlation:
\[ \hat{\rho}_{ij} = \frac{\sum_{u\in U(i,j)}(r_{ui} - b_{ui})(r_{uj}-b_{uj})} {\sqrt{\sum_{u\in U(i,j)}(r_{ui} - b_{ui})^2\sum_{u\in U(i,j)}(r_{uj}-b_{uj})^2}}, \]
where
When \(b_{ui}\) is replaced by user \(u\)’s mean rating, this is the adjusted cosine similarity you will use in the activity.
And using biases we can also do a better job of estimating ratings:
\[ \hat{r}_{ui} = b_{ui} + \frac{\sum_{j \in n_k(i, u)} s_{ij}(r_{uj} - b_{uj})}{\sum_{j \in n_k(i, u)} s_{ij}}, \]
where
When using correlation coefficient as the similarity score, negative values have important implications for the weighted calculation.
Looking at the weighted rating formula:
\[ \hat{r}_{ui} = b_{ui} + \frac{\sum_{j \in n_k(i, u)} s_{ij}(r_{uj} - b_{uj})}{\sum_{j \in n_k(i, u)} s_{ij}}, \]
When \(s_{ij}\) is negative (negative correlation), here’s what happens:
Numerator: The term \(s_{ij}(r_{uj} - b_{uj})\) becomes negative when \(s_{ij} < 0\). This means:
This makes intuitive sense: if items \(i\) and \(j\) are negatively correlated, then liking one means you’re likely to dislike the other.
Denominator: When \(s_{ij}\) is negative, it contributes a negative value to the sum in the denominator. This can potentially lead to issues:
To address these issues, practitioners typically use one of these approaches:
Use only positive similarities: Filter out items with negative correlation (only use \(k\)-NN where \(s_{ij} > 0\))
Use absolute values in denominator: \(\displaystyle\hat{r}_{ui} = b_{ui} + \frac{\sum_{j \in n_k(i, u)} s_{ij}(r_{uj} - b_{uj})}{\sum_{j \in n_k(i, u)} |s_{ij}|}\)
The first approach (filtering to positive similarities) is most common because negatively correlated items are typically less informative for prediction than positively correlated items in most domains.
You will sometimes see a third option, rescaling to \([0, 1]\) with \(s'_{ij} = (s_{ij} + 1)/2\). Avoid it: it gives an uncorrelated neighbour weight \(0.5\) and an anti-correlated one weight near \(0\), so it rewards neighbours that carry no information.
This completes the high level view of CF.
Working with user-user similarities is analogous.
Strengths:
Weaknesses:
Note that standard CF forces us to consider similarity among items, or among users, but does not take into account both.
Can we use both kinds of similarity simultaneously?
We can’t use both the rows and columns of the ratings matrix \(R\) at the same time – the user and item vectors live in different vector spaces.
What we could try to do is find a single vector space in which we represent both users and items, along with a similarity function, such that:
We saw this idea in the SVD lecture: Low Rank Defines Latent Factors. That lecture ended with “the SVD needs every entry” — here is the fix.
This new vector space is called a latent space, and the user and item representations are called latent vectors.
The same shared latent space is the core of deep-learning recommenders (Naumov et al. 2019), covered in the optional module Recommender Systems II.
Now, however, we are working with a matrix which is only partially observed. That is, we only know some of the entries in the ratings matrix.
Nonetheless, we can imagine a situation like this:
where we decompose the ratings matrix \(R\) into two matrices. Each row of the left factor is an item vector \(\mathbf v_i\); each column of the right factor is a user vector \(\mathbf u_u\).
We want the product of the two matrices to be as close as possible to the known values of the ratings matrix.
What this setup implies is that our similarity function is the inner product.
Which means that to predict an unknown rating, we take the inner product of latent vectors:
Taking, for example, item 2 (the 2nd row of the item factor) and user 5 (the 5th column of the user factor)…
We have
\[ \hat r_{u=5,\,i=2} = \mathbf u_5^T \mathbf v_2 = (-0.5 \cdot -2)+(0.6 \cdot 0.3)+(0.5 \cdot 2.4) = 1.0 + 0.18 + 1.2 = 2.38, \]
so (rounded):
Notice that in this case we’ve decided that the factorization should be rank 3, i.e., low-rank.
So we want something like an SVD.
(Recall that SVD gives us the most-accurate-possible low-rank factorization of a matrix).
However, we can’t use the SVD algorithm directly, because we don’t know all the entries in \(R\).
Indeed, the unseen entries in \(R\) are exactly what we want to predict.
Here is what we want to solve:
\[ \min_{U,V} \sum_{(u,i)\in S} \left(r_{ui} - \mathbf u_u^T \mathbf v_i\right)^2 + \lambda(\Vert U\Vert^2 + \Vert V\Vert^2), \]
or in matrix form \(\min_{U,V} \Vert (R - UV^T)_S\Vert^2 + \lambda(\Vert U\Vert^2 + \Vert V\Vert^2)\),
where:
The \((\cdot)_S\) notation means that we are only considering the subset of matrix entries that correspond to known reviews (the set \(S\)).
Note that as usual, we add \(\ell_2\) penalization to avoid overfitting (Ridge regression).
This problem is not convex in \(U\) and \(V\) together, because of the product \(UV^T\). It is biconvex: convex in \(U\) when \(V\) is held fixed, and convex in \(V\) when \(U\) is held fixed.
In particular, if we hold either \(U\) or \(V\) constant, then the result is a simple ridge regression.
So one commonly used algorithm for this problem is called alternating least squares (ALS):
Hold \(V\) fixed. The objective splits into one small problem per user. For user \(u\), let \(J_u\) be the items she rated, \(V_{J_u}\) their rows of \(V\), and \(\mathbf r_{u,J_u}\) her ratings:
\[ \min_{\mathbf u_u} \Vert \mathbf r_{u,J_u} - V_{J_u}\mathbf u_u \Vert^2 + \lambda \Vert \mathbf u_u\Vert^2 \quad\Longrightarrow\quad \mathbf u_u = \left(V_{J_u}^T V_{J_u} + \lambda I\right)^{-1} V_{J_u}^T\, \mathbf r_{u,J_u}. \]
This is the ridge formula \(\hat{\boldsymbol\beta} = (X^TX + cI)^{-1}X^T\mathbf y\) with \(X = V_{J_u}\) (the item vectors are the features), \(\mathbf y = \mathbf r_{u,J_u}\) (her ratings) and \(c = \lambda\).
The other half-step is the mirror image: for item \(i\), regress the ratings it received on the factors of the users who rated it.
solve_U in scripts/recommender_als.py is exactly this line, once per user: np.linalg.solve(V[W[i]].T.dot(V[W[i]]) + I, V[W[i]].T.dot(x_i[W[i]])).
Take \(k = 1\), so every factor is a number. Fixed item factors: \(v_1 = 1\), \(v_2 = 2\), \(v_3 = 1.75\).
User \(u\) rated \(r_{u1} = 2\) and \(r_{u2} = 5\), and has not rated item 3.
\[ J(u_u) = (2 - u_u)^2 + (5 - 2u_u)^2 + \lambda u_u^2 \quad\Longrightarrow\quad u_u = \frac{1\cdot 2 + 2 \cdot 5}{1^2 + 2^2 + \lambda} = \frac{12}{5 + \lambda}. \]
| \(u_u\) | \(\hat r_{u1}\) | \(\hat r_{u2}\) | \(\hat r_{u3}\) | |
|---|---|---|---|---|
| \(\lambda = 0\) | 2.4 | 2.4 | 4.8 | 4.2 |
| \(\lambda = 1\) | 2.0 | 2.0 | 4.0 | 3.5 |
The penalty pulls the factor toward 0, and every prediction shrinks with it: shrinkage. It costs some fit on the two ratings she gave, and protects against trusting a factor estimated from very few ratings.
The same numbers through the ridge formula:
V_rated = np.array([[1.0], [2.0]]) # item factors of the items she rated
r_u = np.array([2.0, 5.0]) # her ratings
for lam in [0, 1]:
u_u = np.linalg.solve(V_rated.T @ V_rated + lam * np.eye(1), V_rated.T @ r_u)
print(f'lambda = {lam}: u_u = {u_u[0]:.2f}, prediction for item 3 = {1.75 * u_u[0]:.2f}')lambda = 0: u_u = 2.40, prediction for item 3 = 4.20
lambda = 1: u_u = 2.00, prediction for item 3 = 3.50The entire Amazon reviews dataset is too large to work with easily, and it is too sparse.
Hence, we will take the densest rows and columns of the matrix: users who wrote more than 50 reviews and movies with more than 50 reviews. Those star ratings are the file dense-block.csv.
dense = pd.read_csv(f'{DATA}/dense-block.csv', dtype={'UserId': str, 'ProductId': str})
dense.head()| UserId | ProductId | Score | |
|---|---|---|---|
| 0 | A02755422E9NI29TCQ5W3 | 0767819438 | 5 |
| 1 | A02755422E9NI29TCQ5W3 | 0767826728 | 5 |
| 2 | A02755422E9NI29TCQ5W3 | 0780651359 | 5 |
| 3 | A02755422E9NI29TCQ5W3 | 0783222734 | 5 |
| 4 | A02755422E9NI29TCQ5W3 | 0783223358 | 5 |
Now we create a \(\textnormal{UserID} \times \textnormal{ProductID}\) matrix from these reviews.
Missing entries will be filled with NaNs.
R = dense.pivot(index='UserId', columns='ProductId', values='Score')
print(f"Review matrix shape: {R.shape}")
print(f"Fraction of entries known: {R.count().sum() / R.size:.2%}")Review matrix shape: (3677, 7244)
Fraction of entries known: 0.99%Even this “dense” block is about 99% empty. A typical \(5 \times 4\) corner is all missing:
R.iloc[900:905, 1000:1004]| ProductId | 1417024917 | 1417030321 | 1417030976 | 1417054069 |
|---|---|---|---|---|
| UserId | ||||
| A1WLZYEOIL1HLT | NaN | NaN | NaN | NaN |
| A1WNJVA59HLMO5 | NaN | NaN | NaN | NaN |
| A1WR12AC35R3K6 | NaN | NaN | NaN | NaN |
| A1WSFHRBY2ZD1R | NaN | NaN | NaN | NaN |
| A1WUMTJOASEL5F | NaN | NaN | NaN | NaN |
We’ll use code from the Antidote Data Framework to do the matrix factorization and ALS. We have local copies recommender_MF.py, recommender_als.py and recommender_lmafit.py in the scripts/ folder of our repository.
To evaluate honestly, we first hide 10% of the known ratings (per user), fit on the rest, and score the hidden ones.
import recommender_MF as MF
np.random.seed(0)
omega_train, omega_val = MF.train_val_split(~R.isnull(), 0.1)
R_train = R.where(omega_train)
print(f'{omega_train.sum().sum():,} training ratings, {omega_val.sum().sum():,} held-out ratings')236,710 training ratings, 26,332 held-out ratingsFit rank 20 with \(\lambda = 1\) on the training ratings:
RS = MF.als_MF(rank=20, lambda_=1)
%time pred, error = RS.fit_model(R_train)
print(f'Training RMSE: {MF.compute_RMSE(pred, R, omega_train):0.3f}')CPU times: user 35.3 s, sys: 165 ms, total: 35.5 s
Wall time: 35.4 s
Training RMSE: 0.311The prediction is a dense matrix — every user gets a score for every movie:
pred.iloc[900:905, 1000:1004].round(2)| ProductId | 1417024917 | 1417030321 | 1417030976 | 1417054069 |
|---|---|---|---|---|
| UserId | ||||
| A1WLZYEOIL1HLT | 3.70 | 3.83 | 4.37 | 4.44 |
| A1WNJVA59HLMO5 | 4.38 | 4.19 | 5.15 | 5.92 |
| A1WR12AC35R3K6 | 4.57 | 3.19 | 4.34 | 4.10 |
| A1WSFHRBY2ZD1R | 5.05 | 4.09 | 4.78 | 4.26 |
| A1WUMTJOASEL5F | 4.75 | 4.37 | 2.88 | 2.16 |
mu = R_train.stack().mean()
baseline = pd.DataFrame(mu, index=R.index, columns=R.columns)
results = [('Global mean $\\mu$', MF.compute_RMSE(baseline, R, omega_train),
MF.compute_RMSE(baseline, R, omega_val))]
results.append(('MF, rank 20, $\\lambda = 1$', MF.compute_RMSE(pred, R, omega_train),
MF.compute_RMSE(pred, R, omega_val)))
for rank, lam in [(5, 1), (5, 10)]:
pred_k, _ = MF.als_MF(rank=rank, lambda_=lam).fit_model(R_train)
results.append((f'MF, rank {rank}, $\\lambda = {lam}$', MF.compute_RMSE(pred_k, R, omega_train),
MF.compute_RMSE(pred_k, R, omega_val)))
table = '| Model | Training RMSE | Held-out RMSE |\n|---|---|---|\n'
table += '\n'.join(f'| {name} | {tr:.3f} | {va:.3f} |' for name, tr, va in results)
display(Markdown(table))| Model | Training RMSE | Held-out RMSE |
|---|---|---|
| Global mean \(\mu\) | 1.168 | 1.168 |
| MF, rank 20, \(\lambda = 1\) | 0.311 | 1.293 |
| MF, rank 5, \(\lambda = 1\) | 0.711 | 1.101 |
| MF, rank 5, \(\lambda = 10\) | 0.884 | 1.006 |
Matrix Factorization per se is a good idea.
However, many of the improvements we’ve discussed for CF apply to MF as well.
To illustrate, we’ll look at some of the successive improvements used by the team that won the Netflix prize (“BellKor’s Pragmatic Chaos”).
When the prize was announced, the Netflix supplied solution achieved an RMSE of 0.951.
By the end of the competition (about 3 years), the winning team’s solution achieved RMSE of 0.856.
Recall our MF objective:
\[ \min_{U, V} \sum_{(u, i)\in S}(r_{ui} - \mathbf u_u^T\mathbf v_i)^2 + \lambda(\Vert U\Vert^2 + \Vert V\Vert^2) \]
If we add biases: \[ \min_{U, V} \sum_{(u, i)\in S}\left(r_{ui} - (\mu + \alpha_u + \beta_i + \mathbf u_u^T\mathbf v_i)\right)^2 + \lambda(\Vert U\Vert^2 + \Vert V\Vert^2 + \Vert \alpha\Vert^2 + \Vert \beta \Vert^2) \]
we see improvements in accuracy:
In reality, ratings are not provided at random.
Take note of which users rated the same movies (ala CF) and use this information.
Older movies tend to get higher ratings!
If we add time-varying biases, and let each user’s taste drift (the item factors stay fixed) (Koren 2009):
\[ \min_{U, V} \sum_{(u, i)\in S}\left(r_{ui} - (\mu + \alpha_u(t) + \beta_i(t) + \mathbf u_u(t)^T\mathbf v_i)\right)^2 + \lambda(\Vert U\Vert^2 + \Vert V\Vert^2 + \Vert \alpha\Vert^2 + \Vert \beta \Vert^2) \]
we see further improvements in accuracy:
With the time and implicit-feedback terms the model is no longer bilinear in two blocks of parameters, so the Netflix teams fit it with stochastic gradient descent (Lecture 16) rather than ALS.
Everything so far needs a history. A new user signs up and has rated nothing.
Mitigations — every real system is a hybrid: