Matrix Factorization with Side Info

Matrix Factorization

One of the most popular techniques for building recommender systems is to frame the problem as matrix completion, in which a large sparse matrix is built containing the ratings that users give to products (in this case, movies), with rows representing users, columns representing items, and entries corresponding to the ratings that they’ve given (e.g. “5 stars”). Most of these entries will be missing, as each users is likely to consume only a handful of the available products (thus, the matrix is sparse), and the goal is to construct a model which would be able to predict the value of the known interactions (i.e. predict which rating would each user give to each movie), which is compared against the observed values. The items to recommend to each user are then the ones with highest predicted values among those which the user has not yet consumed.

Typically, the problem is approached by trying to approximate the interactions matrix as the product of two lower-dimension matrices (a.k.a. latent factor matrices), which when multiplied by each other would produce something that resembles the original matrix, having the nice property that it will produce predictions for all user-item combinations - i.e.

\[ \mathbf{X} \approx \mathbf{A} \mathbf{B}^T \] Where:

• \(\mathbf{X}\) is the interactions matrix (users are rows, items are columns).
• \(\mathbf{A}\) and \(\mathbf{B}\) are the matrices estimated by the model (a.k.a. latent factors), which have a low number of columns, typically 30-100.

For a better and more stable model, the \(\mathbf{X}\) matrix is typically centered by substracting its mean, a bias/intercept is added for each user and item, and a regularization penalty is applied to the model matrices and biases (typically on the L2 norm) - i.e.:

\[ \mathbf{X} \approx \mathbf{A} \mathbf{B}^T + \mu + \mathbf{b}_A + \mathbf{b}_B \] Where:

• \(\mu\) is the global mean used to center \(\mathbf{X}\).
• \(\mathbf{b}_A\) are user-specific biases (row vector).
• \(\mathbf{b}_B\) are item-specific biases (column vector).

The matrices are typically fitted by initializing them to random numbers and then iteratively updating them in a way that decreases the reconstruction error with respect to the observed entries in \(\mathbf{X}\), using either gradient-based procedures (e.g. stochastic gradient descent) or the ALS (alternating least-squares) method, which optimizes one matrix at a time while leaving the other fixed, performing a few sweeps until convergence.

This library (cmfrec) will by default use the ALS method with L2 regularization, and will use user/item biases which are model parameters (updated at each iteration) rather than being pre-estimated.

Loading the data

The MovieLens100k data is taken from the recommenderlab package. As the data is sparse, it is represented as sparse matrices from the Matrix package. The data comes in CSC format, whereas cmfrec requires COO/triplets format - the conversion is handled by the MatrixExtra package for convenience, which also provides extra slicing functionality that will be used later.

library(cmfrec)
library(Matrix)
library(MatrixExtra)
library(recommenderlab)

data("MovieLense")
X <- as.coo.matrix(MovieLense@data)
str(X)
#> Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
#>   ..@ i       : int [1:99392] 0 1 4 5 9 12 14 15 16 17 ...
#>   ..@ j       : int [1:99392] 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:99392] 5 4 4 4 4 3 1 5 4 5 ...
#>   ..@ factors : list()

subsample_coo_matrix <- function(X, indices) {
    X@i <- X@i[indices]
    X@j <- X@j[indices]
    X@x <- X@x[indices]
    return(X)
}

n_ratings <- length(X@x)
set.seed(123)
ix_train <- sample(n_ratings, floor(0.75 * n_ratings), replace=FALSE)
X_train <- subsample_coo_matrix(X, ix_train)
X_test <- subsample_coo_matrix(X, -ix_train)

Classical model

Now fitting the classical matrix factorization model, with global mean centering, user/item biases, L2 regularization which scales with the number of ratings for each user/item, and no side information. This is the model explained in the earlier section: \[ \mathbf{X} \approx \mathbf{A} \mathbf{B}^T + \mu + \mathbf{b}_A + \mathbf{b}_B \]

model.classic <- CMF(X_train, k=25, lambda=0.1, scale_lam=TRUE, verbose=FALSE)

print_rmse <- function(X_test, X_hat, model_name) {
  rmse <- sqrt(mean( (X_test@x - X_hat@x)^2 ))
  cat(sprintf("RMSE for %s is: %.4f\n", model_name, rmse))
}

pred_classic <- predict(model.classic, X_test)
print_rmse(X_test, pred_classic, "classic model")
#> RMSE for classic model is: 0.9236

model.baseline <- MostPopular(X_train, lambda=10, scale_lam=FALSE)
pred_baseline <- predict(model.baseline, X_test)
print_rmse(X_test, pred_baseline, "non-personalized model")
#> RMSE for non-personalized model is: 0.9460

model.improved <- CMF(X_train, k=25, lambda=0.1, scale_lam=TRUE,
                      add_implicit_features=TRUE, w_main=0.75, w_implicit=0.25,
                      use_cg=FALSE, niter=30, verbose=FALSE)
pred_improved <- predict(model.improved, X_test)
print_rmse(X_test, pred_improved, "improved classic model")

#> RMSE for improved classic model is: 0.9126

Adding side information

Collective matrix factorization extends the classical model by incorporating side information about users/items into the formula, which is done by also factorizing the side information matrices, sharing the same latent components that are used for factorizing the \(\mathbf{X}\) matrix: \[ \mathbf{X} \approx \mathbf{A} \mathbf{B}^T + \mu + \mathbf{b}_A + \mathbf{b}_B \] \[ \mathbf{U} \approx \mathbf{A} \mathbf{C}^T + \mu_U \] \[ \mathbf{I} \approx \mathbf{B} \mathbf{D}^T + \mu_I \] \[ \mathbf{I}_x \approx \mathbf{A} \mathbf{B}^T_i \:\:\:\: \mathbf{I}^T_x \approx \mathbf{B} \mathbf{A}^T_i \] Where:

• \(\mathbf{U}\) is a matrix representing side information about users, with each user being a row, and columns corresponding to their attributes.
• \(\mathbf{I}\) is similarly a matrix representing side information about items.
• \(\mathbf{C}\) and \(\mathbf{D}\) are new latent factor matrices used for factorizing the side information matrices, but are not used directly for \(\mathbf{X}\).
• \(\mu_U\) and \(\mu_I\) are column means for the attributes, which are used in order to center them.

Informally, the latent factors now need to explain both the interactions data as well as the side information, thereby making them generalize better to unseen data. This library in addition allows controlling aspects such as the weight that each factorization has in the optimization objective, different regularization for each matrix, having factors that are not shared, among others.

Fetching the side information from recommenderlab:

U <- MovieLenseUser
U$id      <- NULL
U$zipcode <- NULL
U$age2    <- U$age^2
### Note that `cmfrec` does not standardize features beyond mean centering
U$age     <- (U$age - mean(U$age)) / sd(U$age)
U$age2    <- (U$age2 - mean(U$age2)) / sd(U$age2)
U <- model.matrix(~.-1, data=U)

I <- MovieLenseMeta
I$title <- NULL
I$url   <- NULL
I$year  <- ifelse(is.na(I$year), median(I$year, na.rm=TRUE), I$year)
I$year2 <- I$year^2
I$year  <- (I$year - mean(I$year)) / sd(I$year)
I$year2 <- (I$year2 - mean(I$year2)) / sd(I$year2)
I <- as.coo.matrix(I)

cat(dim(U), "\n")
#> 943 24
cat(dim(I), "\n")
#> 1664 21

Now fitting the model:

model.w.sideinfo <- CMF(X_train, U=U, I=I, NA_as_zero_item=TRUE,
                        k=25, lambda=0.1, scale_lam=TRUE,
                        niter=30, use_cg=FALSE, include_all_X=FALSE,
                        w_main=0.75, w_user=0.5, w_item=0.5, w_implicit=0.5,
                        center_U=FALSE, center_I=FALSE,
                        verbose=FALSE)
pred_side_info <- predict(model.w.sideinfo, X_test)
print_rmse(X_test, pred_side_info, "model with side info")

#> RMSE for model with side info is: 0.9117

Summary:

library(kableExtra)

calc_rmse <- function(X_test, X_hat) {
    return(sqrt(mean( (X_test@x - X_hat@x)^2 )))
}
results <- data.frame(
    NonPersonalized = calc_rmse(X_test, pred_baseline),
    ClassicalModel = calc_rmse(X_test, pred_classic),
    ClassicPlusImplicit = calc_rmse(X_test, pred_improved),
    CollectiveModel = calc_rmse(X_test, pred_side_info)
)
results <- as.data.frame(t(results))
names(results) <- "RMSE"
results %>%
    kable() %>%
    kable_styling()

Important to keep in mind:

• These RMSEs have high standard errors due to the small amount of data used here.
• The model hyperparameters are not particularly tuned, and a proper tuning should use a validation split too.
• The test split is using users and items which might not have been in the training set.
• While it looks like the difference from adding side information is very small, it also comes with the side effect of being able to recommend items based on attributes.
• RMSE as a metric can hide overfitting in models that tend to recommend items with too few ratings/interactions - these models in particular will tend to recommend many movies with only a handful ratings, which is typically undesirable. A model with higher regularization that shows higher test RMSE might in practice produce better quality recommendations (see the introductory Python notebook for more examples on this topic).
• The models evaluated so far have all used dynamic/scaled regularization (as proposed in Large-scale Parallel Collaborative Filtering for the Netflix Prize), save for the baseline most-popular model - this means that the regularization for each user and item is scaled by the number of present entries for it. This setting tends to produce dubious recommendations in small datasets like the MovieLens100k, even if it makes it look like it improves RMSE.

Generating Top-N recommendations

The goal behind building a collaborative filtering model is typically to be able to make top-N recommended lists for users or to obtain latent factors for an unseen user given its current data. cmfrec has many prediction functions for these purposes depending on what specifically one wants to do, supporting both warm-start and cold-start recommendations.

### Re-fitting the earlier model to all the data,
### this time *without* scaled regularization
model.classic <- CMF(X, k=20, lambda=10, scale_lam=FALSE, verbose=FALSE)
model.w.sideinfo <- CMF(X, U=U, I=I, k=20, lambda=10, scale_lam=FALSE,
                        w_main=0.75, w_user=0.125, w_item=0.125,
                        verbose=FALSE)

user_to_recommend <- 10
### Note: slicing of 'X' is provided by 'MatrixExtra',
### returning a 'sparseVector' object as required by cmfrec
topN(model.classic, user=user_to_recommend, n=10,
     exclude=X[user_to_recommend, , drop=TRUE])
#>  [1] 316 424 511 311 271 314 405  79 524 190

### A handy function for visualizing recommendations
movie_names <- colnames(X)
n_ratings <- colSums(as.csc.matrix(X, binary=TRUE))
avg_ratings <- colSums(as.csc.matrix(X)) / n_ratings
print_recommended <- function(rec, txt) {
    cat(txt, ":\n",
        paste(paste(1:length(rec), ". ", sep=""),
              movie_names[rec],
              " - Avg rating:", round(avg_ratings[rec], 2),
              ", #ratings: ", n_ratings[rec],
              collapse="\n", sep=""),
        "\n", sep="")
}
recommended <- topN(model.w.sideinfo, user=user_to_recommend, n=5,
                    exclude=X[user_to_recommend, , drop=TRUE])
print_recommended(recommended, "Recommended for user_id=10")
#> Recommended for user_id=10:
#> 1. Schindler's List (1993) - Avg rating:4.47, #ratings: 298
#> 2. To Kill a Mockingbird (1962) - Avg rating:4.29, #ratings: 219
#> 3. Close Shave, A (1995) - Avg rating:4.49, #ratings: 112
#> 4. Boot, Das (1981) - Avg rating:4.2, #ratings: 201
#> 5. Titanic (1997) - Avg rating:4.25, #ratings: 350

recommended_new <- topN_new(model.w.sideinfo, n=5,
                            exclude=X[user_to_recommend, , drop=TRUE],
                            X=X[user_to_recommend, , drop=TRUE],
                            U=U[user_to_recommend, , drop=TRUE])
print_recommended(recommended_new, "Recommended for user_id=10 as new user")
#> Recommended for user_id=10 as new user:
#> 1. Schindler's List (1993) - Avg rating:4.47, #ratings: 298
#> 2. To Kill a Mockingbird (1962) - Avg rating:4.29, #ratings: 219
#> 3. Close Shave, A (1995) - Avg rating:4.49, #ratings: 112
#> 4. Boot, Das (1981) - Avg rating:4.2, #ratings: 201
#> 5. Titanic (1997) - Avg rating:4.25, #ratings: 350

recommended_new <- topN_new(model.w.sideinfo, n=5,
                            exclude=X[user_to_recommend, , drop=TRUE],
                            X=X[user_to_recommend, , drop=TRUE])
print_recommended(recommended_new, "Recommended for user_id=10 as new user (NO sideinfo)")
#> Recommended for user_id=10 as new user (NO sideinfo):
#> 1. Schindler's List (1993) - Avg rating:4.47, #ratings: 298
#> 2. To Kill a Mockingbird (1962) - Avg rating:4.29, #ratings: 219
#> 3. Close Shave, A (1995) - Avg rating:4.49, #ratings: 112
#> 4. Boot, Das (1981) - Avg rating:4.2, #ratings: 201
#> 5. Titanic (1997) - Avg rating:4.25, #ratings: 350

recommended_cold <- topN_new(model.w.sideinfo, n=5,
                             exclude=X[user_to_recommend, , drop=TRUE],
                             U=U[user_to_recommend, , drop=TRUE])
print_recommended(recommended_cold, "Recommended for user_id=10 as new user (NO ratings)")
#> Recommended for user_id=10 as new user (NO ratings):
#> 1. Schindler's List (1993) - Avg rating:4.47, #ratings: 298
#> 2. Close Shave, A (1995) - Avg rating:4.49, #ratings: 112
#> 3. Wrong Trousers, The (1993) - Avg rating:4.47, #ratings: 118
#> 4. Good Will Hunting (1997) - Avg rating:4.26, #ratings: 198
#> 5. Wallace & Gromit: The Best of Aardman Animation (1996) - Avg rating:4.45, #ratings: 67