23.3  Matrix Factorization

Matrix factorization represents each user and item by a low-dimensional vector and predicts a rating from their inner product (Koren et al. 2009). The method became prominent through the Netflix Prize, where factor models were important components of the winning ensemble (Töscher et al. 2009). Its lasting value is conceptual as well as practical: a sparse interaction matrix can be modeled through two compact sets of embeddings, with parameters shared across all observed ratings.

23.3.1 The Matrix Factorization Model

Let \(\mathbf R\in\mathbb R^{m\times n}\) contain the explicit ratings of \(m\) users for \(n\) items. We approximate it by \(\hat{\mathbf R}=\mathbf P\mathbf Q^\top\), where \(\mathbf P\in\mathbb R^{m\times k}\) contains user embeddings and \(\mathbf Q\in\mathbb R^{n\times k}\) contains item embeddings, with \(k\ll\min(m,n)\).

The row \(\mathbf p_u\) represents user \(u\), and \(\mathbf q_i\) represents item \(i\). Their inner product is large when the two embeddings align. Individual coordinates should not usually be read as fixed semantic attributes: rotating both embedding spaces by the same orthogonal matrix leaves every prediction unchanged. The model identifies a useful relative geometry, not a unique label for each latent coordinate. The predicted matrix is

\[\hat{\mathbf{R}} = \mathbf{PQ}^\top \tag{23.3.1}\]

The plain inner product omits systematic offsets: some users give higher ratings than others, and some items receive higher ratings across users. User and item biases model these effects,

\[ \hat{\mathbf{R}}_{ui} = \mathbf{p}_u\mathbf{q}^\top_i + b_u + b_i \]

We fit the model only on the observed set \(\mathcal K\) and regularize all learned parameters:

\[ \underset{\mathbf P,\mathbf Q,\mathbf b^{(u)},\mathbf b^{(i)}}{\mathrm{argmin}} \sum_{(u,i)\in\mathcal K}(R_{ui}-\hat R_{ui})^2 +\lambda\left(\|\mathbf P\|_F^2+\|\mathbf Q\|_F^2 +\|\mathbf b^{(u)}\|_2^2+\|\mathbf b^{(i)}\|_2^2\right). \]

The loss does not treat missing entries as zero ratings. The coefficient \(\lambda\) penalizes embedding and bias magnitudes, and the parameters can be fitted with stochastic gradient methods.

Matrix factorization looks up a user embedding \(\mathbf p_u\) and an item embedding \(\mathbf q_i\). Their inner product, together with user and item biases, predicts \(R_{ui}\); the parameters share information across the sparse interaction matrix.

PyTorch

from d2l import torch as d2l
import torch
from torch import nn

MXNet

from d2l import mxnet as d2l
from mxnet import autograd, gluon, np, npx
from mxnet.gluon import nn
import mxnet as mx
npx.set_np()

23.3.2 Model Implementation

First, we implement the matrix factorization model described above. The user and item latent factors can be created with the nn.Embedding. The input_dim is the number of items/users and the output_dim is the dimension of the latent factors \(k\). We can also use nn.Embedding to create the user/item biases by setting the output_dim to one. In the forward function, user and item ids are used to look up the embeddings.

PyTorch

class MF(nn.Module):
    def __init__(self, num_factors, num_users, num_items):
        super().__init__()
        self.P = nn.Embedding(num_users, num_factors)
        self.Q = nn.Embedding(num_items, num_factors)
        self.user_bias = nn.Embedding(num_users, 1)
        self.item_bias = nn.Embedding(num_items, 1)
        nn.init.normal_(self.P.weight, std=0.01)
        nn.init.normal_(self.Q.weight, std=0.01)
        nn.init.zeros_(self.user_bias.weight)
        nn.init.zeros_(self.item_bias.weight)

    def forward(self, user_id, item_id):
        P_u = self.P(user_id)
        Q_i = self.Q(item_id)
        b_u = self.user_bias(user_id)
        b_i = self.item_bias(item_id)
        outputs = (P_u * Q_i).sum(dim=1) + b_u.squeeze() + b_i.squeeze()
        return outputs.flatten()

MXNet

class MF(nn.Block):
    def __init__(self, num_factors, num_users, num_items):
        super().__init__()
        self.P = nn.Embedding(input_dim=num_users, output_dim=num_factors)
        self.Q = nn.Embedding(input_dim=num_items, output_dim=num_factors)
        self.user_bias = nn.Embedding(num_users, 1)
        self.item_bias = nn.Embedding(num_items, 1)

    def forward(self, user_id, item_id):
        P_u = self.P(user_id)
        Q_i = self.Q(item_id)
        b_u = self.user_bias(user_id)
        b_i = self.item_bias(item_id)
        outputs = (P_u * Q_i).sum(axis=1) + np.squeeze(b_u) + np.squeeze(b_i)
        return outputs.flatten()

23.3.3 Evaluation Measures

We then implement the RMSE (root-mean-square error) measure, which is commonly used to measure the differences between rating scores predicted by the model and the actually observed ratings (ground truth) (Gunawardana and Shani 2015). RMSE is defined as:

\[ \textrm{RMSE} = \sqrt{\frac{1}{|\mathcal{T}|}\sum_{(u, i) \in \mathcal{T}}(\mathbf{R}_{ui} -\hat{\mathbf{R}}_{ui})^2} \]

where \(\mathcal{T}\) is the set consisting of pairs of users and items that you want to evaluate on. \(|\mathcal{T}|\) is the size of this set. We can use MXNet’s RMSE function provided by mx.gluon.metric.

PyTorch

def evaluator(net, test_iter, devices):
    net.eval()
    sq_err_total, n = 0.0, 0
    with torch.no_grad():
        for users, items, ratings in test_iter:
            users, items, ratings = (users.to(devices[0]),
                                     items.to(devices[0]),
                                     ratings.to(devices[0]))
            preds = net(users, items)
            sq_err_total += ((preds - ratings) ** 2).sum().item()
            n += ratings.numel()
    return (sq_err_total / n) ** 0.5

MXNet

def evaluator(net, test_iter, devices):
    rmse = mx.gluon.metric.RMSE()  # Get the RMSE
    rmse_list = []
    for idx, (users, items, ratings) in enumerate(test_iter):
        u = gluon.utils.split_and_load(users, devices, even_split=False)
        i = gluon.utils.split_and_load(items, devices, even_split=False)
        r_ui = gluon.utils.split_and_load(ratings, devices, even_split=False)
        r_hat = [net(u, i) for u, i in zip(u, i)]
        rmse.update(labels=r_ui, preds=r_hat)
        rmse_list.append(rmse.get()[1])
    return float(np.mean(np.array(rmse_list)))

23.3.4 Training and Evaluating the Model

In the training function, we adopt the \(\ell_2\) loss with weight decay. The weight decay mechanism has the same effect as the \(\ell_2\) regularization.

PyTorch

def train_recsys_rating(net, train_iter, test_iter, loss, trainer, num_epochs,
                        devices=d2l.try_all_gpus(), evaluator=None,
                        **kwargs):
    net = net.to(devices[0])
    timer = d2l.Timer()
    animator = d2l.Animator(xlabel='epoch', xlim=[1, num_epochs], ylim=[0, 2],
                            legend=['train loss', 'test RMSE'])
    for epoch in range(num_epochs):
        net.train()
        metric, l = d2l.Accumulator(3), 0.
        for i, values in enumerate(train_iter):
            timer.start()
            users, items, ratings = [v.to(devices[0]) for v in values]
            trainer.zero_grad()
            preds = net(users, items)
            ls = loss(preds, ratings)
            ls.backward()
            trainer.step()
            l += ls.item()
            metric.add(ls.item() * users.shape[0], users.shape[0],
                       users.shape[0])
            timer.stop()
        if len(kwargs) > 0:  # It will be used in section AutoRec
            test_rmse = evaluator(net, test_iter, kwargs['inter_mat'],
                                  devices)
        else:
            test_rmse = evaluator(net, test_iter, devices)
        train_l = l / (i + 1)
        animator.add(epoch + 1, (train_l, test_rmse))
    print(f'train loss {metric[0] / metric[1]:.3f}, '
          f'test RMSE {test_rmse:.3f}')
    print(f'{metric[2] * num_epochs / timer.sum():.1f} examples/sec '
          f'on {str(devices)}')

MXNet

def train_recsys_rating(net, train_iter, test_iter, loss, trainer, num_epochs,
                        devices=d2l.try_all_gpus(), evaluator=None,
                        **kwargs):
    timer = d2l.Timer()
    animator = d2l.Animator(xlabel='epoch', xlim=[1, num_epochs], ylim=[0, 2],
                            legend=['train loss', 'test RMSE'])
    for epoch in range(num_epochs):
        metric, l = d2l.Accumulator(3), 0.
        for i, values in enumerate(train_iter):
            timer.start()
            input_data = []
            values = values if isinstance(values, list) else [values]
            for v in values:
                input_data.append(gluon.utils.split_and_load(v, devices))
            train_feat = input_data[:-1] if len(values) > 1 else input_data
            train_label = input_data[-1]
            with autograd.record():
                preds = [net(*t) for t in zip(*train_feat)]
                ls = [loss(p, s) for p, s in zip(preds, train_label)]
            [l.backward() for l in ls]
            l += sum([l.asnumpy() for l in ls]).mean() / len(devices)
            trainer.step(values[0].shape[0])
            metric.add(l, values[0].shape[0], values[0].size)
            timer.stop()
        if len(kwargs) > 0:  # It will be used in section AutoRec
            test_rmse = evaluator(net, test_iter, kwargs['inter_mat'],
                                  devices)
        else:
            test_rmse = evaluator(net, test_iter, devices)
        train_l = l / (i + 1)
        animator.add(epoch + 1, (train_l, test_rmse))
    print(f'train loss {metric[0] / metric[1]:.3f}, '
          f'test RMSE {test_rmse:.3f}')
    print(f'{metric[2] * num_epochs / timer.sum():.1f} examples/sec '
          f'on {str(devices)}')

Finally, let’s put all things together and train the model. Here, we set the latent factor dimension to 30.

PyTorch

devices = d2l.try_all_gpus()
num_users, num_items, train_iter, test_iter = d2l.split_and_load_ml100k(
    test_ratio=0.1, batch_size=512)
net = MF(30, num_users, num_items)
loss = nn.MSELoss()
trainer = torch.optim.Adam(net.parameters(), lr=0.002, weight_decay=1e-5)
train_recsys_rating(net, train_iter, test_iter, loss, trainer, num_epochs=20,
                    devices=devices, evaluator=evaluator)

train loss 0.683, test RMSE 1.048
327975.9 examples/sec on [device(type='cuda', index=0)]

MXNet

devices = d2l.try_all_gpus()
num_users, num_items, train_iter, test_iter = d2l.split_and_load_ml100k(
    test_ratio=0.1, batch_size=512)
net = MF(30, num_users, num_items)
net.initialize(ctx=devices, force_reinit=True, init=mx.init.Normal(0.01))
lr, num_epochs, wd, optimizer = 0.002, 20, 1e-5, 'adam'
# `gluon.loss.L2Loss()` returns 0.5 * MSE, so with identical `lr` MX would
# train at half the effective gradient of PT's `nn.MSELoss()`. Scale by 2
# to match PyTorch's mean-MSE convention.
loss = gluon.loss.L2Loss(weight=2)
trainer = gluon.Trainer(net.collect_params(), optimizer,
                        {"learning_rate": lr, 'wd': wd})
train_recsys_rating(net, train_iter, test_iter, loss, trainer, num_epochs,
                    devices, evaluator)

train loss 0.129, test RMSE 1.057
203730.6 examples/sec on [gpu(0)]

Below, we use the trained model to predict the rating that a user (ID 20) might give to an item (ID 30).

PyTorch

scores = net(torch.tensor([20], device=devices[0]),
             torch.tensor([30], device=devices[0]))
scores

tensor([2.8977], device='cuda:0', grad_fn=<AddBackward0>)

MXNet

scores = net(np.array([20], dtype='int', ctx=devices[0]),
             np.array([30], dtype='int', ctx=devices[0]))
scores

array([2.9645655], device=gpu(0))

23.3.5 Summary

• The matrix factorization model is widely used in recommender systems. It can be used to predict ratings that a user might give to an item.
• We can implement and train matrix factorization for recommender systems.

23.3.6 Exercises

• Vary the size of latent factors. How does the size of latent factors influence the model performance?
• Try different optimizers, learning rates, and weight decay rates.
• Check the predicted rating scores of other users for a specific movie.

PyTorch

Discussions

MXNet

Discussions

Gunawardana, Asela, and Guy Shani. 2015. “Evaluating Recommender Systems.” In Recommender Systems Handbook. Springer. https://doi.org/10.1007/978-1-4899-7637-6_8.

Koren, Yehuda, Robert Bell, and Chris Volinsky. 2009. “Matrix Factorization Techniques for Recommender Systems.” Computer 42 (8): 30–37. https://doi.org/10.1109/mc.2009.263.

Töscher, Andreas, Michael Jahrer, and Robert M Bell. 2009. The Bigchaos Solution to the Netflix Grand Prize. https://www.netflixprize.com/assets/GrandPrize2009_BPC_BigChaos.pdf.