initial work on biases for model with implicit feedback (#54)

* initial work on biases for model with implicit feedback

* more work on implicit with biases - now rhs takes into account unobserved interactions.

* use precomputed XtX in transform

* fix bug in transform when modelling with biases

* patches wrathematics/float#36

* add tests for models with biases
* fixed bug in transform when using dynamic lambda
* added extra transform() call in fit_transform() to ensure results from transform() and fit_transform() are identical

Author: dselivanov

Diff for: R/RcppExports.R

@@ -73,12 +73,12 @@ dense_csc_prod <- function(x_r, y_csc_r, num_threads = 1L) {
     .Call(`_rsparse_dense_csc_prod`, x_r, y_csc_r, num_threads)
 }
 
-als_implicit_double <- function(m_csc_r, X, Y, XtX, lambda, n_threads, solver, cg_steps, is_x_bias_last_row) {
-    .Call(`_rsparse_als_implicit_double`, m_csc_r, X, Y, XtX, lambda, n_threads, solver, cg_steps, is_x_bias_last_row)
+als_implicit_double <- function(m_csc_r, X, Y, XtX, lambda, n_threads, solver, cg_steps, with_biases, is_x_bias_last_row) {
+    .Call(`_rsparse_als_implicit_double`, m_csc_r, X, Y, XtX, lambda, n_threads, solver, cg_steps, with_biases, is_x_bias_last_row)
 }
 
-als_implicit_float <- function(m_csc_r, X_, Y_, XtX_, lambda, n_threads, solver, cg_steps, is_x_bias_last_row) {
-    .Call(`_rsparse_als_implicit_float`, m_csc_r, X_, Y_, XtX_, lambda, n_threads, solver, cg_steps, is_x_bias_last_row)
+als_implicit_float <- function(m_csc_r, X_, Y_, XtX_, lambda, n_threads, solver, cg_steps, with_biases, is_x_bias_last_row) {
+    .Call(`_rsparse_als_implicit_float`, m_csc_r, X_, Y_, XtX_, lambda, n_threads, solver, cg_steps, with_biases, is_x_bias_last_row)
 }
 
 als_explicit_double <- function(m_csc_r, X, Y, cnt_X, lambda, n_threads, solver, cg_steps, dynamic_lambda, with_biases, is_x_bias_last_row) {

Diff for: R/model_WRMF.R

@@ -85,15 +85,22 @@ WRMF = R6::R6Class(
                           ...) {
       stopifnot(is.null(init) || is.matrix(init))
       solver = match.arg(solver)
-      private$non_negative = ifelse(solver == "nnls", TRUE, FALSE)
       feedback = match.arg(feedback)
 
       if (feedback == 'implicit') {
         # FIXME
-        # now only support bias for explicit feedback
-        with_user_item_bias = FALSE
+
+        if (solver == "conjugate_gradient" && with_user_item_bias == TRUE) {
+          msg = paste("'conjugate_gradient' is not supported for a model",
+            "`with_user_item_bias == TRUE`. Setting to 'cholesky'."
+          )
+          warning(msg)
+          solver = "cholesky"
+        }
         with_global_bias = FALSE
       }
+      private$non_negative = ifelse(solver == "nnls", TRUE, FALSE)
+
       if (private$non_negative && with_global_bias == TRUE) {
         logger$warn("setting `with_global_bias=FALSE` for 'nnls' solver")
         with_global_bias = FALSE
@@ -257,9 +264,10 @@ WRMF = R6::R6Class(
           item_bias = float(n_item)
         }
 
-        self$global_bias = private$init_user_item_bias(c_ui, c_iu, user_bias, item_bias)
+        global_bias = private$init_user_item_bias(c_ui, c_iu, user_bias, item_bias)
         self$components[1L, ] = item_bias
         private$U[private$rank, ] = user_bias
+        if(private$with_global_bias) self$global_bias = global_bias
       } else if (private$feedback == "explicit" && private$with_global_bias) {
         self$global_bias = mean(c_ui@x)
         c_ui@x = c_ui@x - self$global_bias
@@ -282,15 +290,23 @@ WRMF = R6::R6Class(
         cnt_u = float::fl(cnt_u)
         cnt_i = float::fl(cnt_i)
       }
+      private$cnt_u = cnt_u
 
       # iterate
       for (i in seq_len(n_iter)) {
+
         # solve for items
-        loss = private$solver(c_ui, private$U, self$components, TRUE, cnt_X=cnt_i)
+        loss = private$solver(c_ui, private$U, self$components,
+                              is_bias_last_row = TRUE,
+                              cnt_X = cnt_i)
+        logger$info("iter %d (items) loss = %.4f", i, loss)
+
         # solve for users
-        loss = private$solver(c_iu, self$components, private$U, FALSE, cnt_X=cnt_u)
+        loss = private$solver(c_iu, self$components, private$U,
+                              is_bias_last_row = FALSE,
+                              cnt_X = cnt_u)
+        logger$info("iter %d (users) loss = %.4f", i, loss)
 
-        logger$info("iter %d loss = %.4f", i, loss)
         if (loss_prev_iter / loss - 1 < convergence_tol) {
           logger$info("Converged after %d iterations", i)
           break
@@ -299,31 +315,28 @@ WRMF = R6::R6Class(
         loss_prev_iter = loss
       }
 
-      rank_ = ifelse(private$with_user_item_bias, private$rank - 1L, private$rank)
-      ridge = fl(diag(x = private$lambda, nrow = rank_, ncol = rank_))
-
-      X = if (private$with_user_item_bias) tcrossprod(self$components[-1L, ]) else self$components
-      private$XtX = tcrossprod(X) + ridge
-
       if (private$precision == "double")
         data.table::setattr(self$components, "dimnames", list(NULL, colnames(x)))
       else
         data.table::setattr(self$components@Data, "dimnames", list(NULL, colnames(x)))
 
-      res = t(private$U)
-      private$U = NULL
 
-      if (private$precision == "double")
-        setattr(res, "dimnames", list(rownames(x), NULL))
-      else
-        setattr(res@Data, "dimnames", list(rownames(x), NULL))
-      res
+      rank_ = ifelse(private$with_user_item_bias, private$rank - 1L, private$rank)
+      ridge = fl(diag(x = private$lambda, nrow = rank_, ncol = rank_))
+      XX = if (private$with_user_item_bias) self$components[-1L, , drop = FALSE] else self$components
+      private$XtX = tcrossprod(XX) + ridge
+
+      # call extra transform to ensure results from transform() and fit_transform()
+      # are the same (due to avoid_cg, etc)
+      # this adds some extra computation, but not a big deal though
+      self$transform(x)
     },
     # project new users into latent user space - just make ALS step given fixed items matrix
     #' @description create user embeddings for new input
     #' @param x user-item iteraction matrix
     #' @param ... not used at the moment
     transform = function(x, ...) {
+
       stopifnot(ncol(x) == ncol(self$components))
       if (private$feedback == "implicit" ) {
         logger$trace("WRMF$transform(): calling `RhpcBLASctl::blas_set_num_threads(1)` (to avoid thread contention)")
@@ -346,7 +359,19 @@ WRMF = R6::R6Class(
         res = float(0, nrow = private$rank, ncol = nrow(x))
       }
 
-      loss = private$solver(t(x), self$components, res, FALSE, private$XtX, avoid_cg=TRUE)
+      if (private$with_user_item_bias) {
+        res[1, ] = if(private$precision == "double") 1.0 else float::fl(1.0)
+      }
+
+      loss = private$solver(
+        t(x),
+        self$components,
+        res,
+        is_bias_last_row = FALSE,
+        XtX = private$XtX,
+        cnt_X = private$cnt_u,
+        avoid_cg = TRUE
+      )
 
       res = t(res)
 
@@ -367,6 +392,7 @@ WRMF = R6::R6Class(
     dynamic_lambda = FALSE,
     rank = NULL,
     non_negative = NULL,
+    cnt_u = NULL,
     # user factor matrix = rank * n_users
     U = NULL,
     # item factor matrix = rank * n_items
@@ -404,15 +430,15 @@ als_implicit = function(
     rank = ifelse(with_user_item_bias, nrow(X) - 1L, nrow(X))
     ridge = fl(diag(x = lambda, nrow = rank, ncol = rank))
     if (with_user_item_bias) {
-      index_row_to_discard = ifelse(is_bias_last_row, rank, 1L)
-      XtX = tcrossprod(X[-index_row_to_discard, ])
+      index_row_to_discard = ifelse(is_bias_last_row, nrow(X), 1L)
+      XX = X[-index_row_to_discard, , drop = FALSE]
     } else {
-      XtX = tcrossprod(X)
+      XX = X
     }
-    XtX = XtX + ridge
+    XtX = tcrossprod(XX) + ridge
   }
   # Y is modified in-place
-  loss = solver(x, X, Y, XtX, lambda, n_threads, solver_code, cg_steps, is_bias_last_row)
+  loss = solver(x, X, Y, XtX, lambda, n_threads, solver_code, cg_steps, with_user_item_bias, is_bias_last_row)
 }
 
 als_explicit = function(

Diff for: src/RcppExports.cpp

@@ -246,8 +246,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // als_implicit_double
-double als_implicit_double(const Rcpp::S4& m_csc_r, arma::mat& X, arma::mat& Y, const arma::mat& XtX, double lambda, unsigned n_threads, unsigned solver, unsigned cg_steps, bool is_x_bias_last_row);
-RcppExport SEXP _rsparse_als_implicit_double(SEXP m_csc_rSEXP, SEXP XSEXP, SEXP YSEXP, SEXP XtXSEXP, SEXP lambdaSEXP, SEXP n_threadsSEXP, SEXP solverSEXP, SEXP cg_stepsSEXP, SEXP is_x_bias_last_rowSEXP) {
+double als_implicit_double(const Rcpp::S4& m_csc_r, arma::mat& X, arma::mat& Y, const arma::mat& XtX, double lambda, unsigned n_threads, unsigned solver, unsigned cg_steps, const bool with_biases, bool is_x_bias_last_row);
+RcppExport SEXP _rsparse_als_implicit_double(SEXP m_csc_rSEXP, SEXP XSEXP, SEXP YSEXP, SEXP XtXSEXP, SEXP lambdaSEXP, SEXP n_threadsSEXP, SEXP solverSEXP, SEXP cg_stepsSEXP, SEXP with_biasesSEXP, SEXP is_x_bias_last_rowSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -259,14 +259,15 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< unsigned >::type n_threads(n_threadsSEXP);
     Rcpp::traits::input_parameter< unsigned >::type solver(solverSEXP);
     Rcpp::traits::input_parameter< unsigned >::type cg_steps(cg_stepsSEXP);
+    Rcpp::traits::input_parameter< const bool >::type with_biases(with_biasesSEXP);
     Rcpp::traits::input_parameter< bool >::type is_x_bias_last_row(is_x_bias_last_rowSEXP);
-    rcpp_result_gen = Rcpp::wrap(als_implicit_double(m_csc_r, X, Y, XtX, lambda, n_threads, solver, cg_steps, is_x_bias_last_row));
+    rcpp_result_gen = Rcpp::wrap(als_implicit_double(m_csc_r, X, Y, XtX, lambda, n_threads, solver, cg_steps, with_biases, is_x_bias_last_row));
     return rcpp_result_gen;
 END_RCPP
 }
 // als_implicit_float
-double als_implicit_float(const Rcpp::S4& m_csc_r, Rcpp::S4& X_, Rcpp::S4& Y_, Rcpp::S4& XtX_, double lambda, unsigned n_threads, unsigned solver, unsigned cg_steps, bool is_x_bias_last_row);
-RcppExport SEXP _rsparse_als_implicit_float(SEXP m_csc_rSEXP, SEXP X_SEXP, SEXP Y_SEXP, SEXP XtX_SEXP, SEXP lambdaSEXP, SEXP n_threadsSEXP, SEXP solverSEXP, SEXP cg_stepsSEXP, SEXP is_x_bias_last_rowSEXP) {
+double als_implicit_float(const Rcpp::S4& m_csc_r, Rcpp::S4& X_, Rcpp::S4& Y_, Rcpp::S4& XtX_, double lambda, unsigned n_threads, unsigned solver, unsigned cg_steps, const bool with_biases, bool is_x_bias_last_row);
+RcppExport SEXP _rsparse_als_implicit_float(SEXP m_csc_rSEXP, SEXP X_SEXP, SEXP Y_SEXP, SEXP XtX_SEXP, SEXP lambdaSEXP, SEXP n_threadsSEXP, SEXP solverSEXP, SEXP cg_stepsSEXP, SEXP with_biasesSEXP, SEXP is_x_bias_last_rowSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -278,8 +279,9 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< unsigned >::type n_threads(n_threadsSEXP);
     Rcpp::traits::input_parameter< unsigned >::type solver(solverSEXP);
     Rcpp::traits::input_parameter< unsigned >::type cg_steps(cg_stepsSEXP);
+    Rcpp::traits::input_parameter< const bool >::type with_biases(with_biasesSEXP);
     Rcpp::traits::input_parameter< bool >::type is_x_bias_last_row(is_x_bias_last_rowSEXP);
-    rcpp_result_gen = Rcpp::wrap(als_implicit_float(m_csc_r, X_, Y_, XtX_, lambda, n_threads, solver, cg_steps, is_x_bias_last_row));
+    rcpp_result_gen = Rcpp::wrap(als_implicit_float(m_csc_r, X_, Y_, XtX_, lambda, n_threads, solver, cg_steps, with_biases, is_x_bias_last_row));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -538,8 +540,8 @@ static const R_CallMethodDef CallEntries[] = {
     {"_rsparse_cpp_glove_partial_fit", (DL_FUNC) &_rsparse_cpp_glove_partial_fit, 6},
     {"_rsparse_csr_dense_tcrossprod", (DL_FUNC) &_rsparse_csr_dense_tcrossprod, 3},
     {"_rsparse_dense_csc_prod", (DL_FUNC) &_rsparse_dense_csc_prod, 3},
-    {"_rsparse_als_implicit_double", (DL_FUNC) &_rsparse_als_implicit_double, 9},
-    {"_rsparse_als_implicit_float", (DL_FUNC) &_rsparse_als_implicit_float, 9},
+    {"_rsparse_als_implicit_double", (DL_FUNC) &_rsparse_als_implicit_double, 10},
+    {"_rsparse_als_implicit_float", (DL_FUNC) &_rsparse_als_implicit_float, 10},
     {"_rsparse_als_explicit_double", (DL_FUNC) &_rsparse_als_explicit_double, 11},
     {"_rsparse_als_explicit_float", (DL_FUNC) &_rsparse_als_explicit_float, 11},
     {"_rsparse_initialize_biases_double", (DL_FUNC) &_rsparse_initialize_biases_double, 8},