Collaborative filtering and matrix factorization

The big idea

In the last lesson, Ana's rating of Interstellar was predicted by finding users who rate like her. That works, but it compares raw rows of a matrix that is more than 99 percent empty in real life, so two users often share almost no rated items at all.

Matrix factorization takes a different route. It assumes that a handful of hidden traits explain most of the ratings. Every film gets a short vector saying how much of each trait it has, and every user gets a vector saying how much they care about each trait. A rating is then roughly the dot product of the two. A film critic might call the traits "serious vs escapist" or "sci-fi vs romance"; the model just calls them factor 1, factor 2, and learns them from the data.

Think of a music shop assistant who has never seen your purchase history, only two numbers: how much you like loud guitars and how much you like slow songs. With those two numbers and the same two numbers for every album, they can guess your opinion of any album in the shop. Matrix factorization learns those numbers for everyone at once.

How it works

Step 1: Neighbourhoods, a quick recap

The collaborative filtering in Recommender fundamentals is a neighbourhood methodCollaborative filtering that predicts from similar users or similar items, using a similarity measure such as cosine on co-rated items.. User-based: find users similar to you and average their ratings. Item-based: find items similar to the ones you liked, where two items are similar if the same users rated them alike. Both are easy to explain ("because you watched Alien"), but similarity is computed only from the few items two rows share, so it is noisy when the matrix is sparse, and the full item-item table grows with the square of the number of items.

Step 2: Latent factors

The model, in its plainest form, gives user $u$ a vector $p_u$ and item $i$ a vector $q_i$, both of length $k$, the number of latent factorsHidden dimensions learned by a factorization model; each user and item gets a coordinate on each, and their dot product predicts the interaction., and predicts

The extra terms are biasesPer-user and per-item offsets in a factorization model that capture generous raters and well-liked items, independent of taste.. Koren, Bell and Volinsky's overview of the Netflix Prize work gives the example: the average rating $\mu$ is 3.7 stars, Titanic tends to be rated 0.5 above average and Joe is a critical user who rates 0.3 below, so the baseline for Joe and Titanic is 3.7 + 0.5 - 0.3 = 3.9 stars. The factors then only have to explain the part of the rating that is about taste.

Training minimizes squared error on the known ratings only, plus a penalty on the size of every vector so that a user with three ratings does not get extreme factors:

This is not the SVD of Eigenvectors, eigenvalues and SVD . Classic SVD needs every cell filled; here the loss simply skips the empty ones, which is exactly what makes it work on sparse data.

Step 3: Training with SGD or ALS

Two ways to minimize that loss:

• Stochastic gradient descent. Loop over the known ratings. For each, compute the error $e_{ui}=r_{ui}-{\hat{r}}_{ui}$ and nudge both vectors: $q_i\leftarrow q_i+\gamma (e_{ui}{p}_u-\lambda q_i)$ and the same for $p_u$. Simon Funk popularized this during the Netflix Prize; it is simple and fast.
• Alternating least squaresAlternating least squares: fix all item vectors and solve each user vector exactly as a small ridge regression, then fix users and solve items, and repeat. (ALS). With both $p$ and $q$ unknown the problem is not convex. But fix all the item vectors and each user's vector becomes an ordinary ridge regression with $k$ unknowns, which has an exact solution. So ALS alternates: solve every user with items fixed, then every item with users fixed. Each half-step can only lower the loss.

Koren and colleagues note that SGD is generally easier and faster, but ALS wins in two cases. First, every user's solve is independent of every other user's, so it parallelizes across a cluster. Second, implicit data, where the "training set" is effectively every cell, which is too many to loop over one by one.

Step 4: Implicit feedback needs weights

With play counts instead of ratings ( Recommender fundamentals ), you cannot just factorize the known cells: they are all positive, so a model trained on them alone learns that everyone likes everything. Hu, Koren and Volinsky's fix, from a study of TV viewing on about 300,000 set-top boxes, splits each cell into two numbers:

• a preference $p_{ui}$: 1 if the user played the item at all, 0 otherwise;
• a confidence $c_{ui}=1+\alpha r_{ui}$, where $r_{ui}$ is the count. A 0 still gets confidence 1, a little, because not watching is weak evidence; a show watched five times with $\alpha =40$ gets confidence 201.

The loss is $\sum c_{ui}(p_{ui}-x_u^{\top }{y}_i{)}^2$ plus regularization, summed over all user-item pairs, missing ones included. The paper reports that $\alpha =40$ worked well on its data. Summing over billions of pairs sounds hopeless, but ALS has a trick: the product $Y^{\top }{C}^{u}Y$ needed for user $u$ equals $Y^{\top }Y+Y^{\top }(C^u-I)Y$. The first part is the same for every user and is computed once; the second only touches the items user $u$ actually played. The cost per sweep becomes linear in the number of non-zero cells, and the paper uses about 10 sweeps.

Try it

1. At k = 2 and lambda 0.10, press 10 sweeps. Training RMSE falls to 0.15 and held-out RMSE settles at 0.59. The empty cells fill in: Ana gets 5.0 for Alien and Dee 2.4 for Interstellar. In the film map, the three sci-fi films (blue) sit apart from the three romances (pink): the factors found the genres on their own.
2. Drag lambda to 0. Training RMSE drops to 0.00, a perfect fit, while held-out RMSE jumps to 1.56: Ben's held-out Heat, which he rated 4, is now predicted at 0.5. With 33 training ratings and 48 free numbers (two factors and a bias for each of 16 users and films), the model memorizes.
3. Set lambda back to 0.10 and k to 0, which leaves only the average and the biases: held-out RMSE is 0.81. Then try k = 1 with lambda 0.30: 0.47, the best here. One taste factor plus biases beats biases alone, and a bit more regularization helps a tiny data set.

Going deeper

Code: explicit ratings with ALS

import numpy as np

# Synthetic ratings: 600 users x 300 films, 3 hidden taste dimensions plus biases, 1-5 stars.
rng = np.random.default_rng(0)
m, n, k_true = 600, 300, 3
U, V = rng.normal(0, 0.8, (m, k_true)), rng.normal(0, 0.8, (n, k_true))
full = np.clip(np.round(3.4 + rng.normal(0, .4, (m, 1)) + rng.normal(0, .5, n) + U @ V.T + rng.normal(0, .5, (m, n))), 1, 5)
seen = rng.random((m, n)) < 0.08                       # each user rated about 8% of films
test = seen & (rng.random((m, n)) < 0.2)               # hold out 20% of the known ratings
train = seen & ~test
mu = full[train].mean()
print(f"{train.sum()} training ratings, {test.sum()} held out, {train.mean():.1%} of cells known")

def als(k, lam, sweeps=15):
    """Explicit-feedback ALS with user/item biases; ridge scaled by each row's count."""
    r = np.random.default_rng(1)
    P, Q = r.normal(0, .1, (m, k)), r.normal(0, .1, (n, k))
    bu, bi = np.zeros(m), np.zeros(n)
    for _ in range(sweeps):
        for rows, B, F, cols, cB, mask in ((range(m), bu, P, Q, bi, train), (range(n), bi, Q, P, bu, train.T)):
            for a in rows:
                obs = mask[a]
                X = np.c_[np.ones(obs.sum()), cols[obs]]          # [1, other side's factors]
                y = (full[a, obs] if mask is train else full[obs, a]) - mu - cB[obs]
                w = np.linalg.solve(X.T @ X + lam * obs.sum() * np.eye(k + 1), X.T @ y)
                B[a], F[a] = w[0], w[1:]
    return mu + bu[:, None] + bi[None, :] + P @ Q.T

rmse = lambda pred, mask: np.sqrt(np.mean((pred[mask] - full[mask]) ** 2))
print(f"global mean only:      test RMSE {rmse(np.full((m, n), mu), test):.3f}")
for k, lam in [(0, 0.1), (3, 0.01), (3, 0.1), (10, 0.01), (10, 0.1)]:
    pred = als(k, lam)
    label = "biases only" if k == 0 else f"k={k:<2d} lambda={lam}"
    print(f"{label:22s} train RMSE {rmse(pred, train):.3f}  test RMSE {rmse(pred, test):.3f}")

11499 training ratings, 2950 held out, 6.4% of cells known
global mean only:      test RMSE 1.193
biases only            train RMSE 1.020  test RMSE 1.092
k=3  lambda=0.01       train RMSE 0.460  test RMSE 0.707
k=3  lambda=0.1        train RMSE 0.503  test RMSE 0.693
k=10 lambda=0.01       train RMSE 0.201  test RMSE 1.387
k=10 lambda=0.1        train RMSE 0.380  test RMSE 0.760

600 synthetic users and 300 films whose ratings come from 3 hidden taste dimensions, biases and noise, with 6.4 percent of cells known. Each configuration runs 15 ALS sweeps.

• Biases alone help a little: test RMSE 1.092 against 1.193 for the global mean.
• Factors help a lot: k = 3 reaches 0.693 with lambda 0.1.
• Too many factors without regularization overfit: k = 10 with lambda 0.01 has the lowest training error, 0.201, and the worst test error, 1.387. Raising lambda to 0.1 brings it back to 0.760. Choose k and lambda on held-out data.

Code: implicit plays with confidence weights

import numpy as np

# Synthetic play counts: 800 listeners x 400 songs. Liked songs get replayed; some plays are accidents.
rng = np.random.default_rng(0)
m, n = 800, 400
taste = rng.normal(0, 1, (m, 4)) @ rng.normal(0, 1, (4, n)) + rng.normal(0, 1, n)
liked = rng.random((m, n)) < 1 / (1 + np.exp(-(taste - 5)))
accident = (rng.random((m, n)) < 0.03) & ~liked                    # one stray play each
plays = np.where(liked, 2 + rng.poisson(4, (m, n)), 0) + accident
test = {u: rng.choice(np.flatnonzero(liked[u])) for u in range(m) if liked[u].sum() >= 3}
R = plays.copy()
for u, i in test.items(): R[u, i] = 0                              # hide one liked song per user
print(f"{(R > 0).mean():.1%} of cells played; one liked song hidden for {len(test)} listeners")

def als(R, k=20, lam=10.0, alpha=1.0, use_missing=True, sweeps=10):
    """Hu, Koren, Volinsky: p = [plays > 0], confidence c = 1 + alpha * plays."""
    r = np.random.default_rng(1)
    X, Y = r.normal(0, .1, (R.shape[0], k)), r.normal(0, .1, (R.shape[1], k))
    for _ in range(sweeps):
        for A, B, M in ((X, Y, R), (Y, X, R.T)):
            BtB = B.T @ B if use_missing else 0                     # every cell counts, at weight 1
            for a in range(len(A)):
                nz = np.flatnonzero(M[a])
                c = 1 + alpha * M[a, nz]
                lhs = BtB + (B[nz].T * (c - use_missing)) @ B[nz] + lam * np.eye(k)
                A[a] = np.linalg.solve(lhs, (B[nz].T * c).sum(axis=1))
    return X @ Y.T

def hit_rate(scores, K=10):                                        # is the hidden song in the top 10?
    s = np.where(R > 0, -np.inf, scores)                           # never re-recommend played songs
    return np.mean([i in np.argsort(-s[u])[:K] for u, i in test.items()])

print(f"most popular songs            hit rate@10 {hit_rate(np.tile((R > 0).sum(0), (m, 1)) * 1.0):.3f}")
print(f"observed cells only           hit rate@10 {hit_rate(als(R, use_missing=False)):.3f}")
for alpha in [0.0, 1.0, 10.0, 40.0]:
    print(f"all cells, alpha={alpha:<4}         hit rate@10 {hit_rate(als(R, alpha=alpha)):.3f}")

6.7% of cells played; one liked song hidden for 793 listeners
most popular songs            hit rate@10 0.115
observed cells only           hit rate@10 0.124
all cells, alpha=0.0          hit rate@10 0.288
all cells, alpha=1.0          hit rate@10 0.252
all cells, alpha=10.0         hit rate@10 0.155
all cells, alpha=40.0         hit rate@10 0.096

800 synthetic listeners and 400 songs. Liked songs get several plays; about 3 percent of the other cells get one accidental play. One liked song per listener is hidden, and the score is the share of listeners whose hidden song appears in their top 10.

• Ignoring the missing cells fails: fitting the played cells only scores 0.124, no better than recommending the most popular songs (0.115).
• Counting every cell works: the same factorization over all cells, with missing ones at preference 0, scores 0.288 at alpha 0.
• Alpha must be tuned. Here the number of replays of a liked song is pure noise, so a large alpha amplifies noise: 0.096 at alpha 40. On its TV data, where repeat viewing did carry signal, the paper reports that alpha = 40 gave good results. Tune alpha together with lambda on held-out data.

At scale

Spark MLlib's ALS implements both variants. Its parameters mirror this lesson: rank (the number of factors, default 10), regParam (lambda), implicitPrefs to switch to the implicit-feedback variant, and alpha for its confidence. Spark scales lambda by the number of ratings each user or item has, the "ALS-WR" approach, which the lab above also uses. For a user or item missing from training it predicts NaN by default, or drops those rows with coldStartStrategy="drop", so evaluation is not broken by cold-start cases.

Trade-offs

Common mistakes

• Factorizing implicit data on the observed cells only. Every cell is a positive, so every prediction comes out as "likes it".
• Forgetting biases, so the factors waste capacity modelling generous raters and blockbuster films.
• Picking k and lambda by training error. Training error always falls as k grows; use a held-out set of ratings.
• Filling missing ratings with zeros or averages and running ordinary SVD. The fill-in values swamp the real ones.
• Expecting factors to have names. Sometimes one lines up with a genre, as in the lab, but usually they are mixtures.

Where you will see this

The Netflix Prize made matrix factorization famous. In 2006 Netflix released more than 100 million ratings from about 500,000 customers on more than 17,000 films and offered a million dollars to the first team that cut the RMSE of its own system by 10 percent. That system had an RMSE of 0.9514, so the target was about 0.856. Koren, Bell and Volinsky, whose team was in first place when they wrote their overview, report that most of their winning predictors were factorization models, and that accuracy kept improving as they added biases, implicit feedback and time effects. Today the implicit-feedback ALS of Hu, Koren and Volinsky ships in Spark MLlib, and the same idea of scoring a user vector against an item vector with a dot product is the core of the two-tower retrieval models in the next lesson.

Check yourself

Recap