Evaluating recommender systems

Implicit-feedback data

Historically, many models for recommender systems were designed by approaching the problem as regression or rating prediction, by taking as input a matrix \(\mathbf{X}_{ui}\) denoting user likes and dislikes of items in a scale (e.g. users giving a 1-to-5 star rating to different movies), and evaluating such models by seeing how well they predict these ratings on hold-out data.

In many cases, it is impossible or very expensive to obtain such data, but one has instead so called “implicit-feedback” records: that is, observed logs of user interactions with items (e.g. number of times that a user played each song in a music service), which do not signal dislikes in the same way as a 1-star rating would, but can still be used for building and evaluating recommender systems.

In the latter case, the problem is approached more as ranking or classification instead of regression, with the models being evaluated not by how well they perform at predicting ratings, but by how good they are at scoring the observed interactions higher than the non-observed interactions for each user, using metrics more typical of information retrieval.

Generating a ranked list of items for each user according to their predicted score and comparing such lists against hold-out data can nevertheless be very slow (might even be slower than fitting the model itself), and this is where recometrics comes in: it provides efficient routines for calculating many implicit-feedback recommendation quality metrics, which exploit multi-threading, SIMD instructions, and efficient sorted search procedures.

Matrix factorization models

The perhaps most common approach towards building a recommendation model is by trying to approximate the matrix \(\mathbf{X}_{mn}\) as the product of two lower-dimensional matrices \(\mathbf{A}_{mk}\) and \(\mathbf{B}_{nk}\) (with \(k \ll m\) and \(k \ll n\)), representing latent user and item factors/components, respectively (which are the model parameters to estimate) - i.e. \[ \mathbf{X} \approx \mathbf{A} \mathbf{B}^T \] In the explicit-feedback setting (e.g. movie ratings), this is typically done by trying to minimize squared errors with respect to the observed entries in \(\mathbf{X}\), while in implicit-feedback settings this is typically done by turning the \(\mathbf{X}\) matrix into a binary matrix which has a one if the observation is observed and a zero if not, using the actual values (e.g. number of times that a song was played) instead as weights for the positive entries, thereby looking at all entries rather than just the observed (non-zero) values - e.g.: \[ \min_{\mathbf{A}, \mathbf{B}} \sum_{u=1}^{m} \sum_{i=1}^{n} x_{ui} (I_{x_{ui}>0} - \mathbf{a}_u \cdot \mathbf{b}_i)^2 \]

The recommendations for a given user are then produced by calculating the full products between that user vector \(\mathbf{a}_u\) and the \(\mathbf{B}\) matrix, sorting these predicted scores in descending order.

For a better overview of implicit-feedback matrix factorization, see the paper Hu, Yifan, Yehuda Koren, and Chris Volinsky. “Collaborative filtering for implicit feedback datasets.” 2008 Eighth IEEE International Conference on Data Mining. Ieee, 2008.

Evaluating recommendation models

Such matrix factorization models are commonly evaluated by setting aside a small amount of users as hold-out for evaluation, fitting a model to all the remaining users and items. Then, from the evaluation users, a fraction of their interactions data is set as a hold-out test set, while their latent factors are computed using the rest of the data and the previously fitted model from the other users.

Then, top-K recommendations for each user are produced, discarding the non-hold-out items with which their latent factors were just determined, and these top-K lists are compared against the hold-out test items, seeing how well they do at ranking them near the top vs. how they rank the remainder of the items.

This package can be used to calculate many recommendation quality metrics given the user and item factors and the train-test data split that was used, including:

• P@K (“precision-at-k”): this is the most intuitive metric. It calculates the proportion of the top-K recommendations that include items from the test set for a given user - i.e. \[ P@K = \frac{1}{k} \sum_{i=1}^k \begin{cases} 1, & r_i \in \mathcal{T}\\ 0, & \text{otherwise} \end{cases} \] Where \(r_i\) is the item ranked at position \(i\) by the model (sorting the predicted scores in descending order, after excluding the items in the training data for that user), and \(\mathcal{T}\) is the set of items that are in the test set for that user.

  Note that some papers and libraries define \(P@K\) differently, see the second version below.

• TP@K (truncated \(P@K\)): same calculation as \(P@K\), but will instead divide by the minimum between \(k\) and the number of test items: \[ TP@K = \frac{1}{\min\{k, |\mathcal{T}|\}} \sum_{i=1}^k \begin{cases} 1, & r_i \in \mathcal{T}\\ 0, & \text{otherwise} \end{cases} \]

  The “truncated” prefix is a non-standard nomenclature introduced here to differentiate it from the other \(P@K\) metric.

• R@K (“recall-at-k”): while \(P@K\) offers an intuitive metric that captures what a recommender system aims at being good at, it does not capture the fact that, the more test items there are, the higher the chances that they will be included in the top-K recommendations. Recall instead looks at what proportion of the test items would have been retrieved with the top-K recommended list: \[ R@K = \frac{1}{|\mathcal{T}|} \sum_{i=1}^k \begin{cases} 1, & r_i \in \mathcal{T}\\ 0, & \text{otherwise} \end{cases} \]

• AP@K (“average precision-at-k”): precision and recall look at all the items in the top-K equally, whereas one might want to take into account also the ranking within this top-K list, for which this metric comes in handy. “Average Precision” tries to reflect the precisions that would be obtained at different recalls: \[ AP@K = \frac{1}{|\mathcal{T}|} \sum_{i=1}^k \begin{cases} P@i, & r_i \in \mathcal{T}\\ 0, & \text{otherwise} \end{cases} \] \(AP@K\) is a metric which to some degree considers precision, recall, and rank within top-K. Intuitively, it tries to approximate the are under a precision-recall tradeoff curve. Its average across users is typically called “MAP@K” or “Mean Average Precision”.

  Important: many authors define \(AP@K\) differently, such as dividing by the minimum between \(k\) and \(|\mathcal{T}|\) instead, or as the average for P@1..P@K (either as-is or stopping the calculation after already retrieving all test items). See below for the other version.

• TAP@K (truncated \(AP@K\)): a truncated version of the \(AP@K\) metric, which will instead divide it by the minimum between \(k\) and the number of test items. Just like for \(TP@K\), the “truncated” prefix is a non-standard nomenclature used here to differentiate it from the other more typical \(AP@K\).

• NDCG@K (“normalized discounted cumulative gain at k”): while the earlier metrics look at just the presence of an item in the test set, these items might not all be as good, with some of them having higher observed values than others. NDCG aims at judging these values, but discounted according to the rank in the top-K list. First it calculates the unstandardized discounted cumulative gain: \[ DCG@K = \sum_{i=1}^{k} \frac{C_{r_i}}{log_2 (1+i)} \] Where \(C_{r_i}\) indicates the observed interaction value in the test data for item \(r_i\), and is zero if the item was not in the test data. The DCG@K metric is then standardized by dividing it by the maximum achievable DCG@K for the test data: \[ NDCG@K = \frac{DCG@K}{\max DCG@K} \]

  Unlike the other metrics, NDCG can handle data which contains “dislikes” in the form of negative values. If there are no negative values in the test data, it will be bounded between zero and one.

• Hit@K (from which “Hit Rate” is calculated): this is a simpler yes/no metric that looks at whether any of the top-K recommended items were in the test set for a given user: \[ Hit@K = \max_{i=1..K} \begin{cases} 1, & r_i \in \mathcal{T}\\ 0, & \text{otherwise} \end{cases} \] The average of this metric across users is typically called “Hit Rate”.

• RR@K (“reciprocal rank at k”, from which “MRR” or “mean reciprocal rank” is calculated): this metric only looks at the rank of the first recommended item that is in the test set, and outputs its inverse: \[ RR@K = \max_{i=1..K} \frac{1}{i} \:\:\:\text{s.t.}\:\:\: r_i \in \mathcal{T} \] The average of this metric across users is typically called “Mean Reciprocal Rank”.

• ROC AUC (“area under the receiver-operating characteristic curve”): see the Wikipedia entry for details. While the metrics above only looked at the top-K recommended items, this metric looks at the full ranking of items instead, and produces a standardized number between zero and one in which 0.5 denotes random predictions.

• PR AUC (“area under the precision-recall curve”): while ROC AUC provides an overview of the overall ranking, one is typically only interested in how well it retrieves test items within top ranks, and for this the area under the precision-recall curve can do a better job at judging rankings, albeit the metric itself is not standardized and its minimum does not go as low as zero.

  The metric is calculated using the fast but not-so-precise rectangular method, whose formula corresponds to the AP@K metric with K=N. Some papers and libraries call this the average of this metric the “MAP” or “Mean Average Precision” instead (without the “@K”).

(For more details about the metrics, see the package documentation: ?calc.reco.metrics)

NOT covered by this package:

• Metrics that look at the rareness of the items recommended (to evaluate so-called “serendipity”).

• Metrics that look at “discoverability”.

• Metrics that take into account the diversity of the ranked lists.

Now a practical example using the library cmfrec and the MovieLens100K data, taken from the recommenderlab package.

Note that this is an explicit-feedback dataset about movie ratings. Here it will be converted to implicit-feedback by setting movies with a rating of 4 and 5 stars as the positive (observed) data, while the others will be set as negative (unobserved).

library(Matrix)
library(MatrixExtra)
library(data.table)
library(kableExtra)
library(recommenderlab)
library(cmfrec)
library(recometrics)

data(MovieLense)
X_raw <- MovieLense@data

### Converting it to implicit-feedback
X_implicit <- as.coo.matrix(filterSparse(X_raw, function(x) x >= 4))
str(X_implicit)
#> Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
#>   ..@ i       : int [1:55024] 0 1 4 5 9 15 16 17 20 22 ...
#>   ..@ j       : int [1:55024] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..@ Dim     : int [1:2] 943 1664
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:943] "1" "2" "3" "4" ...
#>   .. ..$ : chr [1:1664] "Toy Story (1995)" "GoldenEye (1995)" "Four Rooms (1995)" "Get Shorty (1995)" ...
#>   ..@ x       : num [1:55024] 5 4 4 4 4 5 4 5 5 5 ...
#>   ..@ factors : list()

reco_split <- create.reco.train.test(
    X_implicit,
    users_test_fraction = NULL,
    max_test_users = 100,
    items_test_fraction = 0.3,
    seed = 123
)
X_train <- reco_split$X_train ## Train data for test users
X_test <- reco_split$X_test ## Test data for test users
X_rem <- reco_split$X_rem ## Data to fit the model
users_test <- reco_split$users_test ## IDs of the test users

### Random recommendations (random latent factors)
set.seed(123)
UserFactors_random <- matrix(rnorm(nrow(X_test) * 5), nrow=5)
ItemFactors_random <- matrix(rnorm(ncol(X_test) * 5), nrow=5)

### Non-personalized recommendations
model_baseline <- cmfrec::MostPopular(as.coo.matrix(X_rem), implicit=TRUE)
item_biases <- model_baseline$matrices$item_bias

### Typical implicit-feedback ALS model
### a.k.a. "WRMF" (weighted regularized matrix factorization)
model_wrmf <- cmfrec::CMF_implicit(as.coo.matrix(X_rem), k=10, verbose=FALSE)
UserFactors_wrmf <- t(cmfrec::factors(model_wrmf, X_train))

### As a comparison, this is the typical explicit-feedback model,
### implemented by software such as Spark,
### and called "Weighted-Lambda-Regularized Matrix Factorization".
### Note that it determines the user factors using the train+test data.
model_wlr <- cmfrec::CMF(as.coo.matrix(X_raw[-users_test, ]),
                         lambda=0.1, scale_lam=TRUE,
                         user_bias=FALSE, item_bias=FALSE,
                         k=10, verbose=FALSE)
UserFactors_wlr <- t(cmfrec::factors(model_wlr, as.csr.matrix(X_raw)[users_test,]))

### This is a different explicit-feedback model which
### uses the same regularization for each user and item
### (as opposed to the "weighted-lambda" model) and which
### adds "implicit features", which are a binarized version
### of the input data, but without weights.
### Note that it determines the user factors using the train+test data.
model_hybrid <- cmfrec::CMF(as.coo.matrix(X_raw[-users_test, ]),
                            lambda=20, scale_lam=FALSE,
                            user_bias=FALSE, item_bias=FALSE,
                            add_implicit_features=TRUE,
                            k=10, verbose=FALSE)
UserFactors_hybrid <- t(cmfrec::factors(model_hybrid, as.csr.matrix(X_raw)[users_test,]))

### Processing user side information
U <- as.data.table(MovieLenseUser)[-users_test, ]
mean_age <- mean(U$age)
sd_age <- sd(U$age)
levels_occ <- levels(U$occupation)
MatrixExtra::restore_old_matrix_behavior()
process.U <- function(U, mean_age,sd_age, levels_occ) {
    U[, `:=`(
        id = NULL,
        age = (age - mean_age) / sd_age,
        sex = as.numeric(sex == "M"),
        occupation =  factor(occupation, levels_occ),
        zipcode = NULL
    )]
    U <- Matrix::sparse.model.matrix(~.-1, data=U)
    U <- as.coo.matrix(U)
    return(U)
}
U <- process.U(U, mean_age,sd_age, levels_occ)
U_train <- as.data.table(MovieLenseUser)[users_test, ]
U_train <- process.U(U_train, mean_age,sd_age, levels_occ)

### Processing item side information
I <- as.data.table(MovieLenseMeta)
mean_year <- mean(I$year, na.rm=TRUE)
sd_year <- sd(I$year, na.rm=TRUE)
I[
    is.na(year), year := mean_year
][, `:=`(
    title = NULL,
    year = (year - mean_year) / sd_year,
    url = NULL
)]
I <- as.coo.matrix(I)

### Manually re-creating a binarized matrix and weights
### that will mimic the WRMF model
X_rem_ones <- as.coo.matrix(mapSparse(X_rem, function(x) rep(1, length(x))))
W_rem <- as.coo.matrix(mapSparse(X_rem, function(x) x+1))
X_train_ones <- as.coo.matrix(mapSparse(X_train, function(x) rep(1, length(x))))
W_train <- as.coo.matrix(mapSparse(X_train, function(x) x+1))

### WRMF model, but with item biases/intercepts
model_bwrmf <- cmfrec::CMF(X_rem_ones, weight=W_rem, NA_as_zero=TRUE,
                           lambda=1, scale_lam=FALSE,
                           center=FALSE, user_bias=FALSE, item_bias=TRUE,
                           k=10, verbose=FALSE)
UserFactors_bwrmf <- t(cmfrec::factors(model_bwrmf, X_train_ones, weight=W_train))


### Collective WRMF model (taking user and item attributes)
model_cwrmf <- cmfrec::CMF_implicit(as.coo.matrix(X_rem), U=U, I=I,
                                    NA_as_zero_user=TRUE, NA_as_zero_item=TRUE,
                                    center_U=TRUE, center_I=TRUE,
                                    lambda=0.1,
                                    k=10, verbose=FALSE)
UserFactors_cwrmf <- t(cmfrec::factors(model_cwrmf, X_train, U=U_train))

### Collective WRMF plus item biases/intercepts
model_bcwrmf <- cmfrec::CMF(X_rem_ones, weight=W_rem, NA_as_zero=TRUE,
                            U=U, I=I, center_U=FALSE, center_I=FALSE,
                            NA_as_zero_user=TRUE, NA_as_zero_item=TRUE,
                            lambda=0.1, scale_lam=FALSE,
                            center=FALSE, user_bias=FALSE, item_bias=TRUE,
                            k=10, verbose=FALSE)
UserFactors_bcwrmf <- t(cmfrec::factors(model_bcwrmf, X_train_ones,
                                        weight=W_train, U=U_train))

k <- 5 ## Top-K recommendations to evaluate

### Baselines
metrics_random <- calc.reco.metrics(X_train, X_test,
                                    A=UserFactors_random,
                                    B=ItemFactors_random,
                                    k=k, all_metrics=TRUE)
metrics_baseline <- calc.reco.metrics(X_train, X_test,
                                      A=NULL, B=NULL,
                                      item_biases=item_biases,
                                      k=k, all_metrics=TRUE)

### Simple models
metrics_wrmf <- calc.reco.metrics(X_train, X_test,
                                  A=UserFactors_wrmf,
                                  B=model_wrmf$matrices$B,
                                  k=k, all_metrics=TRUE)
metrics_wlr <- calc.reco.metrics(X_train, X_test,
                                 A=UserFactors_wlr,
                                 B=model_wlr$matrices$B,
                                 k=k, all_metrics=TRUE)
metrics_hybrid <- calc.reco.metrics(X_train, X_test,
                                    A=UserFactors_hybrid,
                                    B=model_hybrid$matrices$B,
                                    k=k, all_metrics=TRUE)

### More complex models
metrics_bwrmf <- calc.reco.metrics(X_train, X_test,
                                   A=UserFactors_bwrmf,
                                   B=model_bwrmf$matrices$B,
                                   item_biases=model_bwrmf$matrices$item_bias,
                                   k=k, all_metrics=TRUE)
metrics_cwrmf <- calc.reco.metrics(X_train, X_test,
                                   A=UserFactors_cwrmf,
                                   B=model_cwrmf$matrices$B,
                                   k=k, all_metrics=TRUE)
metrics_bcwrmf <- calc.reco.metrics(X_train, X_test,
                                    A=UserFactors_bcwrmf,
                                    B=model_bcwrmf$matrices$B,
                                    item_biases=model_bcwrmf$matrices$item_bias,
                                    k=k, all_metrics=TRUE)

metrics_baseline %>%
    head(5) %>%
    kable() %>%
    kable_styling()

all_metrics <- list(
    `Random` = metrics_random,
    `Non-personalized` = metrics_baseline,
    `Weighted-Lambda` = metrics_wlr,
    `Hybrid-Explicit` = metrics_hybrid,
    `WRMF (a.k.a. iALS)` = metrics_wrmf,
    `bWRMF` = metrics_bwrmf,
    `CWRMF` = metrics_cwrmf,
    `bCWRMF` = metrics_bcwrmf
)
results <- all_metrics %>%
    lapply(function(df) as.data.table(df)[, lapply(.SD, mean)]) %>%
    data.table::rbindlist() %>%
    as.data.frame()
row.names(results) <- names(all_metrics)

results %>%
    kable() %>%
    kable_styling()