Confidence

Now we will create cofidence matrices for train and cross-validation:

make_confidence = function(x, alpha) {
  x_confidence = x
  stopifnot(inherits(x, "sparseMatrix"))
  x_confidence@x = 1 + alpha * x@x
  x_confidence
}

Creating matrices:

alpha = 0.1
X_train_conf = make_confidence(X_train, alpha)
X_cv_history_conf = make_confidence(X_cv_history, alpha)

It make sense to try different values of alpha – we will try it later.

Note that we don’t want to touch X_cv_future since we will use raw listening counts as relevance scores in ndcg calculation. If we will transform them we won’t be able to compare ndcg of the models wich use different confidence transformations because ndcg values will depend on (alpha).

Adding scoring metrices to model – ndcg@10 and map@10:

model$add_scorer(x = X_cv_history_conf, y = X_cv_future, name = "ndcg-10", metric = "ndcg@10")
model$add_scorer(x = X_cv_history_conf, y = X_cv_future, name = "map-10", metric = "map@10")

Training

When using reco with optimized multithreaded BLAS such as Apple VecLib, OpenBLAS, Intel MKL it is very important to disable BLAS parallelism since reco itself already uses multiple threads at higher level. Due to thread contention and nesting such programs can be slower then single thread program by a factor of 1000 and more! In order to set number of BLAS threads to one I will use RhpcBLASctl package.

RhpcBLASctl::blas_set_num_threads(1)
user_embeddings = model$fit_transform(X_train_conf, n_iter = 10L, n_threads = 8)

Check trace of the training:

library(ggplot2)
trace = attr(user_embeddings, "trace")
g = ggplot(trace) +  geom_line(aes(x = iter, y = value, col = scorer))
plotly::ggplotly(g, width = 9, height = NULL)

As we can see algorithm quickly converges (only about 5 iterations) and loss convergence is highly correlated to ndcg and map convergence. More interesting it seems algorithm is quite resistant to overfitting.

Let’s create grid search in order to find optimal hyperparameters of (alpha) and rank. Usually models with higher rank work better. But also the larger the rank the stronger should be regularization (lambda). Here we omitting tuning of (lambda) for simplicity.

futile.logger::flog.threshold(futile.logger::ERROR)
RhpcBLASctl::blas_set_num_threads(1)
trace = NULL
alpha = c(0.01, 1)
rank = c(16, 8)
n_iter_max = 10L
n_threads = 8L
grid = expand.grid(alpha = alpha, rank = rank)
convergence_tol = 0.01
for(k in seq_len(nrow(grid))) {
  alpha = grid$alpha[[k]]
  rank = grid$rank[[k]]
  # futile.logger::flog.info("alpha = %.3f, rank = %d", alpha, rank)
  model = ALS$new(rank = rank)
  X_train_conf = make_confidence(X_train, alpha)
  X_cv_history_conf = make_confidence(X_cv_history, alpha)
  model$add_scorer(x = X_cv_history_conf, y = X_cv_future, name = "ndcg-10", metric = "ndcg@10")
  model$add_scorer(x = X_cv_history_conf, y = X_cv_future, name = "map-10", metric = "map@10")
  # fit model
  user_embeddings = model$fit_transform(X_train_conf, n_iter = n_iter_max, 
                                        convergence_tol = convergence_tol, n_threads = n_threads)
  # store strace
  grid_trace = attr(user_embeddings, "trace")
  grid_trace$param_set = sprintf("alpha=%.3f; rank=%d", alpha, rank)
  trace = c(trace, list(grid_trace))
}
trace = rbindlist(trace)

Plotting traces:

g = ggplot(trace) + 
  geom_line(aes(x = iter, y = value, col = param_set)) +
  facet_wrap( ~ scorer, scales = "free") +
  theme(legend.position="bottom")
plotly::ggplotly(g, width = 9, height = NULL)

Calculation of map and ndcg for each iteration is not always necessary because as we’ve seen above it is highly correlated with loss. Also it takes quite long time to compute such metrices. So usually it worth just to check them after loss converges. This is easy with $ndcg_k() and $ap_k() methods:

ndcg_k = model$ndcg_k(x = X_cv_history_conf, y = X_cv_future, k = 10)
mean_ndcg_k = mean(ndcg_k)
ap_k = model$ap_k(x = X_cv_history_conf, y = X_cv_future, k = 10)
mean_ap_k = mean(ap_k)

Recommendations for new users

This question is not particularly well explained in articles or blog posts I found. Moreover it looks like big source of confusion – a lot of questions at StackOverflow: 1, 2, 3, 4, 5.

One obvious solution will be to add new users interactions and run matrix factoriation again. However such method won’t work for real world problems – full factorization is very expensive operation and we need recommendations in near real-time.

But we can make very good approximation of the idea above. Remember that after convergence we got 2 matrices (U) and (I) – users and items embeddings. They are calculated from large amount of user-item interactions and each of them unlikely will change too much after addition of couple of users or items.

Reminder – our goal is to create recommendations for user which was not in our training data. In order to do this we need to:

1. Take embeddings (I) from full matrix factorization
2. Somehow obtain user embeddings (U_j)
3. Calculate (score = U_j I) and sort items according to scores values

So the task will be following: given a new user interactions and precalculated fixed item embeddings (from full matrix factorization) we need to obtain user embeddings. How to do that? Little linear algebra will help – we need just to run single iteration of our ALS algorithm. This is exactly how we minimized cost during full matrix factorization – alternated between fixing user factors and solving for item factors and fixing item factors and solving for user factors. So given item factor we can calculate new users factors.

And this is how reco’s transform method works. predict method takes user embeddings from transform method, multiplies with item embeddings and takes top of them. This is exactly how we made recommendations for our cross-validation.

Lets go through example:

Obtain new user embeddings:

model = ALS$new(rank = 16)
set.seed(1)
alpha = 0.01
X_train_conf = make_confidence(X_train, alpha)
X_cv_history_conf = make_confidence(X_cv_history, alpha)
user_embeddings = model$fit_transform(X_train_conf, n_iter = 10L, n_threads = 8)
new_user_embeddings = model$transform(X_cv_history_conf)

And let’s check how to make predictions for a single user and how much time it will take:

library(microbenchmark)
new_user_1 = X_cv_history_conf[1:1, , drop = FALSE]
microbenchmark({new_user_predictions = model$predict(new_user_1, k = 10)}, times = 10)

As we can see it is single digit millisecond, so we can serve them real-time. In the next post I will describe how to build scalable high-performance RESTful service for serving recommendations in real-time.

Similar items

Another very useful type of recommendations is item-to-item recommendations. A good example is Amazon which show similar items on a product page. And interestingly we can create such recommendations by using item embeddings and calculationg cosine distance between them:

# transpose for convenience of cosine distance calculation
artist_emb = t(model$components)
artist_query = artist_emb["the beatles", , drop = FALSE]
artist_query_sim = text2vec::sim2(artist_query, artist_emb, method = "cosine")
# flatten 1-row matrix
artist_query_sim = artist_query_sim[1, ]
# check similar artists
sort(artist_query_sim, decreasing = T)[1:10]

I can’t say why “Bob Dylan” or “The Rolling Stones” are ahead of “John Lennon”. Keep in mind that we didn’t solve optimization task for item “similarity” (and notion of similarity itself is ambigious). But they look reasonalbe and we got them complimentary to our user-item recommendations.

Here is another example for “The Offspring”:

# transpose for convenience of cosine distance calculation
artist_emb = t(model$components)
artist_query = artist_emb["the offspring", , drop = FALSE]
artist_query_sim = text2vec::sim2(artist_query, artist_emb, method = "cosine")
# flatten 1-row matrix
artist_query_sim = artist_query_sim[1, ]
# check similar artists
sort(artist_query_sim, decreasing = T)[1:10]