diff --git a/R/MatrixFactorizationRecommender.R b/R/MatrixFactorizationRecommender.R
index 9a5d844..2498c5f 100644
--- a/R/MatrixFactorizationRecommender.R
+++ b/R/MatrixFactorizationRecommender.R
@@ -8,6 +8,9 @@ MatrixFactorizationRecommender = R6::R6Class(
     components = NULL,
     #' @field global_bias global mean (for centering values in explicit feedback)
     global_bias = 0.,
+    #' @field global_bias_base Pre-calculated `-(factors*global_bias)` (for centering values in
+    #' implicit feedback when not using user/item biases)
+    global_bias_base = numeric(),
     #' @description recommends items for users
     #' @param x user-item interactions matrix (usually sparse - `Matrix::sparseMatrix`).Users are
     #' rows and items are columns
diff --git a/R/RcppExports.R b/R/RcppExports.R
index c5315cc..ffa5d30 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -117,19 +117,19 @@ als_explicit_float <- function(m_csc_r, X_, Y_, cnt_X_, lambda, n_threads, solve
     .Call(`_rsparse_als_explicit_float`, m_csc_r, X_, Y_, cnt_X_, lambda, n_threads, solver, cg_steps, dynamic_lambda, with_biases, 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_double <- function(m_csc_r, X, Y, XtX, lambda, n_threads, solver, cg_steps, with_biases, is_x_bias_last_row, global_bias, global_bias_base, initialize_bias_base) {
+    .Call(`_rsparse_als_implicit_double`, m_csc_r, X, Y, XtX, lambda, n_threads, solver, cg_steps, with_biases, is_x_bias_last_row, global_bias, global_bias_base, initialize_bias_base)
 }
 
-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_implicit_float <- function(m_csc_r, X_, Y_, XtX_, lambda, n_threads, solver, cg_steps, with_biases, is_x_bias_last_row, global_bias, global_bias_base_, initialize_bias_base) {
+    .Call(`_rsparse_als_implicit_float`, m_csc_r, X_, Y_, XtX_, lambda, n_threads, solver, cg_steps, with_biases, is_x_bias_last_row, global_bias, global_bias_base_, initialize_bias_base)
 }
 
-initialize_biases_double <- function(m_csc_r, m_csr_r, user_bias, item_bias, lambda, dynamic_lambda, non_negative, calculate_global_bias = FALSE) {
-    .Call(`_rsparse_initialize_biases_double`, m_csc_r, m_csr_r, user_bias, item_bias, lambda, dynamic_lambda, non_negative, calculate_global_bias)
+initialize_biases_double <- function(m_csc_r, m_csr_r, user_bias, item_bias, lambda, dynamic_lambda, non_negative, calculate_global_bias = FALSE, is_explicit_feedback = FALSE, initialize_item_biases = FALSE) {
+    .Call(`_rsparse_initialize_biases_double`, m_csc_r, m_csr_r, user_bias, item_bias, lambda, dynamic_lambda, non_negative, calculate_global_bias, is_explicit_feedback, initialize_item_biases)
 }
 
-initialize_biases_float <- function(m_csc_r, m_csr_r, user_bias, item_bias, lambda, dynamic_lambda, non_negative, calculate_global_bias = FALSE) {
-    .Call(`_rsparse_initialize_biases_float`, m_csc_r, m_csr_r, user_bias, item_bias, lambda, dynamic_lambda, non_negative, calculate_global_bias)
+initialize_biases_float <- function(m_csc_r, m_csr_r, user_bias, item_bias, lambda, dynamic_lambda, non_negative, calculate_global_bias = FALSE, is_explicit_feedback = FALSE, initialize_item_biases = FALSE) {
+    .Call(`_rsparse_initialize_biases_float`, m_csc_r, m_csr_r, user_bias, item_bias, lambda, dynamic_lambda, non_negative, calculate_global_bias, is_explicit_feedback, initialize_item_biases)
 }
 
diff --git a/R/model_WRMF.R b/R/model_WRMF.R
index 1b5869a..bdee07e 100644
--- a/R/model_WRMF.R
+++ b/R/model_WRMF.R
@@ -87,18 +87,6 @@ WRMF = R6::R6Class(
       solver = match.arg(solver)
       feedback = match.arg(feedback)
 
-      if (feedback == 'implicit') {
-        # FIXME
-
-        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) {
@@ -141,7 +129,10 @@ WRMF = R6::R6Class(
             precision = private$precision,
             with_user_item_bias = private$with_user_item_bias,
             is_bias_last_row = is_bias_last_row,
-            XtX = XtX)
+            global_bias = self$global_bias,
+            initialize_bias_base = !avoid_cg,
+            XtX = XtX,
+            global_bias_base = self$global_bias_base)
         } else {
           als_explicit(
             x, X, Y, cnt_X,
@@ -162,7 +153,8 @@ WRMF = R6::R6Class(
                      initialize_biases_double,
                      initialize_biases_float)
         FUN(c_ui, c_iu, user_bias, item_bias, private$lambda, private$dynamic_lambda,
-            private$non_negative, private$with_global_bias)
+            private$non_negative, private$with_global_bias, feedback == "explicit",
+            private$solver_code != 1)
       }
 
       self$components = init
@@ -211,28 +203,43 @@ WRMF = R6::R6Class(
         private$U = large_rand_matrix(private$rank, n_user)
         # for item biases
         if (private$with_user_item_bias) {
-          private$U[1, ] = rep(1.0, n_user)
+          private$U[1L, ] = rep(1.0, n_user)
         }
       } else {
         private$U = flrnorm(private$rank, n_user, 0, 0.01)
         if (private$with_user_item_bias) {
-          private$U[1, ] = float::fl(rep(1.0, n_user))
+          private$U[1L, ] = float::fl(rep(1.0, n_user))
         }
       }
 
       if (is.null(self$components)) {
-        if (private$precision == "double") {
-          self$components = large_rand_matrix(private$rank, n_item)
-          # for user biases
-          if (private$with_user_item_bias) {
-            self$components[private$rank, ] = rep(1.0, n_item)
+
+        if (private$solver_code == 1L) { ### <- cholesky
+          if (private$precision == "double") {
+            self$components = matrix(0, private$rank, n_item)
+            if (private$with_user_item_bias) {
+              self$components[private$rank, ] = rep(1.0, n_item)
+            }
+          } else {
+            self$components = float::float(0, private$rank, n_item)
+            if (private$with_user_item_bias) {
+              self$components[private$rank, ] = float::fl(rep(1.0, n_item))
+            }
           }
-        } else {
-          self$components = flrnorm(private$rank, n_item, 0, 0.01)
-          if (private$with_user_item_bias) {
-            self$components[private$rank, ] = float::fl(rep(1.0, n_item))
+          } else {
+            if (private$precision == "double") {
+              self$components = large_rand_matrix(private$rank, n_item)
+              # for user biases
+              if (private$with_user_item_bias) {
+                self$components[private$rank, ] = rep(1.0, n_item)
+              }
+            } else {
+              self$components = flrnorm(private$rank, n_item, 0, 0.01)
+              if (private$with_user_item_bias) {
+                self$components[private$rank, ] = float::fl(rep(1.0, n_item))
+              }
+            }
           }
-        }
       } else {
         stopifnot(is.matrix(self$components) || is.float(self$components))
         stopifnot(ncol(self$components) == n_item)
@@ -268,10 +275,23 @@ WRMF = R6::R6Class(
         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
-        c_iu@x = c_iu@x - self$global_bias
+      } else if (private$with_global_bias) {
+        if (private$feedback == "explicit") {
+          self$global_bias = mean(c_ui@x)
+          c_ui@x = c_ui@x - self$global_bias
+          c_iu@x = c_iu@x - self$global_bias
+        } else {
+          s = sum(c_ui@x)
+          self$global_bias = s / (s + as.numeric(nrow(c_ui))*as.numeric(ncol(c_ui)) - length(c_ui@x))
+        }
+      }
+
+      if (private$feedback == "implicit") {
+        size_global_bias_base = ifelse(private$with_user_item_bias, 0L, private$rank-1L)
+        if (private$precision == "double")
+          self$global_bias_base = numeric(size_global_bias_base)
+        else
+          self$global_bias_base = float(size_global_bias_base)
       }
 
       logger$info("starting factorization with %d threads", getOption("rsparse_omp_threads", 1L))
@@ -350,7 +370,7 @@ WRMF = R6::R6Class(
 
       x = as(x, "CsparseMatrix")
       x = private$preprocess(x)
-      if (self$global_bias != 0.)
+      if (self$global_bias != 0. && private$feedback == "explicit")
         x@x = x@x - self$global_bias
 
       if (private$precision == "double") {
@@ -420,7 +440,10 @@ als_implicit = function(
   precision,
   with_user_item_bias,
   is_bias_last_row,
-  XtX = NULL) {
+  initialize_bias_base,
+  global_bias = 0.,
+  XtX = NULL,
+  global_bias_base = NULL) {
 
   solver = ifelse(precision == "float",
                   als_implicit_float,
@@ -437,8 +460,15 @@ als_implicit = function(
     }
     XtX = tcrossprod(XX) + ridge
   }
+  if (is.null(global_bias_base)) {
+    global_bias_base = numeric()
+    if (precision == "float")
+      global_bias_base = float::fl(global_bias_base)
+  }
   # Y is modified in-place
-  loss = solver(x, X, Y, XtX, lambda, n_threads, solver_code, cg_steps, with_user_item_bias, 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, global_bias,
+                global_bias_base, initialize_bias_base)
 }
 
 als_explicit = function(
diff --git a/inst/include/wrmf_explicit.hpp b/inst/include/wrmf_explicit.hpp
index bd97c29..9bb5746 100644
--- a/inst/include/wrmf_explicit.hpp
+++ b/inst/include/wrmf_explicit.hpp
@@ -32,8 +32,8 @@ arma::Col<T> cg_solver_explicit(const arma::Mat<T>& X_nnz, const arma::Col<T>& c
 
 template <class T>
 T als_explicit(const dMappedCSC& Conf, arma::Mat<T>& X, arma::Mat<T>& Y,
-               const double lambda, const unsigned n_threads, const unsigned solver,
-               const unsigned cg_steps, const bool dynamic_lambda,
+               const double lambda, const int n_threads, const unsigned int solver,
+               const unsigned int cg_steps, const bool dynamic_lambda,
                const arma::Col<T>& cnt_X, const bool with_biases,
                const bool is_x_bias_last_row) {
   /* Note about biases:
@@ -63,7 +63,7 @@ T als_explicit(const dMappedCSC& Conf, arma::Mat<T>& X, arma::Mat<T>& Y,
       x_biases = X.row(0).t();
   }
 
-  T loss = 0;
+  double loss = 0;
   size_t nc = Conf.n_cols;
 #ifdef _OPENMP
 #pragma omp parallel for num_threads(n_threads) schedule(dynamic, GRAIN_SIZE) reduction(+:loss)
diff --git a/inst/include/wrmf_implicit.hpp b/inst/include/wrmf_implicit.hpp
index 92de512..693da67 100644
--- a/inst/include/wrmf_implicit.hpp
+++ b/inst/include/wrmf_implicit.hpp
@@ -31,11 +31,68 @@ arma::Col<T> cg_solver_implicit(const arma::Mat<T>& X_nnz, const arma::Col<T>& c
   return x;
 }
 
+/* FIXME: 'global_bias_base' used like this has very poor numerical precision and ends up making things worse */
+template <class T>
+arma::Col<T> cg_solver_implicit_global_bias(const arma::Mat<T>& X_nnz, const arma::Col<T>& confidence,
+                                            const arma::Col<T>& x_old, const arma::uword n_iter,
+                                            const arma::Mat<T>& XtX,   const arma::Col<T> global_bias_base,
+                                            T global_bias) {
+  arma::Col<T> x = x_old;
+  const arma::Col<T> confidence_1 = confidence - 1.0;
+
+  arma::Col<T> Ap;
+  arma::Col<T> r = X_nnz * (confidence - (confidence_1 % (X_nnz.t() * x + global_bias))) - XtX * x + global_bias_base;
+  arma::Col<T> p = r;
+  double rsold, rsnew, alpha;
+  rsold = arma::dot(r, r);
+
+  for (auto k = 0; k < n_iter; k++) {
+    Ap = XtX * p + X_nnz * (confidence_1 % (X_nnz.t() * p));
+    alpha = rsold / dot(p, Ap);
+    x += alpha * p;
+    r -= alpha * Ap;
+    rsnew = dot(r, r);
+    if (rsnew < CG_TOL) break;
+    p = r + p * (rsnew / rsold);
+    rsold = rsnew;
+  }
+  return x;
+}
+
+template <class T>
+arma::Col<T> cg_solver_implicit_user_item_bias(const arma::Mat<T>& X_nnz, const arma::Col<T>& confidence,
+                                               const arma::Col<T>& x_old, const arma::uword n_iter,
+                                               const arma::Mat<T>& XtX,   const arma::Col<T> &rhs_init,
+                                               const arma::Col<T> &x_biases, T global_bias) {
+  arma::Col<T> x = x_old;
+  const arma::Col<T> confidence_1 = confidence - 1.0;
+
+  arma::Col<T> Ap;
+  arma::Col<T> r = X_nnz * (confidence - (confidence_1 % (X_nnz.t() * x + x_biases + global_bias)))
+                    - XtX * x + rhs_init;
+  arma::Col<T> p = r;
+  double rsold, rsnew, alpha;
+  rsold = arma::dot(r, r);
+
+  for (auto k = 0; k < n_iter; k++) {
+    Ap = XtX * p + X_nnz * (confidence_1 % (X_nnz.t() * p));
+    alpha = rsold / dot(p, Ap);
+    x += alpha * p;
+    r -= alpha * Ap;
+    rsnew = dot(r, r);
+    if (rsnew < CG_TOL) break;
+    p = r + p * (rsnew / rsold);
+    rsold = rsnew;
+  }
+  return x;
+}
+
 template <class T>
 T als_implicit(const dMappedCSC& Conf, arma::Mat<T>& X, arma::Mat<T>& Y,
-               const arma::Mat<T>& XtX, double lambda, unsigned n_threads,
-               unsigned solver, unsigned cg_steps, const bool with_biases,
-               bool is_x_bias_last_row) {
+               const arma::Mat<T>& XtX, double lambda, int n_threads,
+               const unsigned int solver, unsigned int cg_steps, const bool with_biases,
+               const bool is_x_bias_last_row, double global_bias,
+               arma::Col<T> &global_bias_base, const bool initialize_bias_base) {
   // if is_x_bias_last_row == true
   // X = [1, ..., x_bias]
   // Y = [y_bias, ..., 1]
@@ -48,6 +105,12 @@ T als_implicit(const dMappedCSC& Conf, arma::Mat<T>& X, arma::Mat<T>& Y,
   arma::Col<T> x_biases;
   arma::Mat<T> rhs_init;
 
+  if (global_bias < std::sqrt(std::numeric_limits<T>::epsilon()))
+    global_bias = 0;
+
+  if (global_bias && initialize_bias_base && !with_biases)
+    global_bias_base = arma::sum(X, 1) * (-global_bias);
+
   if (with_biases) {
     if (is_x_bias_last_row)  // last row
       x_biases = X.row(X.n_rows - 1).t();
@@ -78,16 +141,27 @@ T als_implicit(const dMappedCSC& Conf, arma::Mat<T>& X, arma::Mat<T>& Y,
     //    X_nnz_user * 1 * (0 - x_biases_nnz_user) +
     //    X_nnz_user * C_nnz_user * (1 - x_biases_nnz_user)
 
-    // here we do following:
-    // drop row with "ones" placeholder for convenient ALS form
-    rhs_init = drop_row<T>(X, is_x_bias_last_row);
-    // p = 0
-    // C = 1 (so we omit multiplication on eye matrix)
-    // rhs = X * eye * (0 - x_biases) = -X * x_biases
-    rhs_init *= -x_biases;
+    if (!global_bias) {
+      // here we do following:
+      // drop row with "ones" placeholder for convenient ALS form
+      rhs_init = drop_row<T>(X, is_x_bias_last_row);
+      // p = 0
+      // C = 1 (so we omit multiplication on eye matrix)
+      // rhs = X * eye * (0 - x_biases) = -X * x_biases
+      rhs_init *= -x_biases;
+    }
+
+    else {
+      rhs_init = - (drop_row<T>(X, is_x_bias_last_row) * (x_biases + global_bias));
+    }
+  }
+
+  else if (global_bias) {
+    rhs_init = arma::Mat<T>(&global_bias_base[0], rank - (int)with_biases, 1, false, true);
   }
 
-  T loss = 0;
+
+  double loss = 0;
   size_t nc = Conf.n_cols;
 #ifdef _OPENMP
 #pragma omp parallel for num_threads(n_threads) schedule(dynamic, GRAIN_SIZE) reduction(+:loss)
@@ -114,7 +188,15 @@ T als_implicit(const dMappedCSC& Conf, arma::Mat<T>& X, arma::Mat<T>& Y,
       arma::Col<T> Y_new;
 
       if (solver == CONJUGATE_GRADIENT) {
-        Y_new = cg_solver_implicit<T>(X_nnz, confidence, init, cg_steps, XtX);
+        if (!with_biases && !global_bias)
+          Y_new = cg_solver_implicit<T>(X_nnz, confidence, init, cg_steps, XtX);
+        else if (with_biases)
+          Y_new = cg_solver_implicit_user_item_bias<T>(X_nnz, confidence, init, cg_steps, XtX,
+                                                       rhs_init, x_biases(idx), global_bias);
+        else
+          Y_new = cg_solver_implicit_global_bias<T>(X_nnz, confidence, init, cg_steps, XtX,
+                                                    rhs_init, global_bias);
+
       } else {
         const arma::Mat<T> lhs =
             XtX + X_nnz.each_row() % (confidence.t() - 1) * X_nnz.t();
@@ -137,6 +219,8 @@ T als_implicit(const dMappedCSC& Conf, arma::Mat<T>& X, arma::Mat<T>& Y,
           // expression above can be simplified further:
           rhs = rhs_init + X_nnz * (confidence - x_biases(idx) % (confidence - 1));
 
+        } else if (global_bias) {
+          rhs = X_nnz * confidence + rhs_init;
         } else {
           rhs = X_nnz * confidence;
         }
diff --git a/inst/include/wrmf_utils.hpp b/inst/include/wrmf_utils.hpp
index ca9be28..c7e5420 100644
--- a/inst/include/wrmf_utils.hpp
+++ b/inst/include/wrmf_utils.hpp
@@ -9,12 +9,33 @@ arma::Mat<T> drop_row(const arma::Mat<T>& X_nnz, const bool drop_last) {
   }
 };
 
+/* https://en.wikipedia.org/wiki/Kahan_summation_algorithm */
 template <class T>
-double initialize_biases(dMappedCSC& ConfCSC,  // modified in place
-                         dMappedCSC& ConfCSR,  // modified in place
-                         arma::Col<T>& user_bias, arma::Col<T>& item_bias, T lambda,
-                         bool dynamic_lambda, bool non_negative,
-                         bool calculate_global_bias) {
+long double compensated_sum(T *arr, int n)
+{
+  long double err = 0.;
+  long double diff = 0.;
+  long double temp;
+  long double res = 0;
+
+  for (int ix = 0; ix < n; ix++)
+  {
+    diff = arr[ix] - err;
+    temp = res + diff;
+    err = (temp - res) - diff;
+    res = temp;
+  }
+
+  return res;
+}
+
+template <class T>
+double initialize_biases_explicit(dMappedCSC& ConfCSC,  // modified in place
+                                  dMappedCSC& ConfCSR,  // modified in place
+                                  arma::Col<T>& user_bias, arma::Col<T>& item_bias, T lambda,
+                                  bool dynamic_lambda, bool non_negative,
+                                  bool calculate_global_bias)
+{
   /* Robust mean calculation */
   double global_bias = 0;
   if (calculate_global_bias) {
@@ -59,3 +80,67 @@ double initialize_biases(dMappedCSC& ConfCSC,  // modified in place
   }
   return global_bias;
 }
+
+template <class T>
+double initialize_biases_implicit(dMappedCSC& ConfCSC, dMappedCSC& ConfCSR,
+                                  arma::Col<T>& user_bias, arma::Col<T>& item_bias,
+                                  T lambda, bool calculate_global_bias, bool non_negative,
+                                  const bool initialize_item_biases)
+{
+  double global_bias = 0;
+  if (calculate_global_bias) {
+    long double s = compensated_sum(ConfCSR.values, ConfCSR.nnz);
+    global_bias = s / (s + (long double)ConfCSR.n_rows*(long double)ConfCSR.n_cols - (long double)ConfCSR.nnz);
+  }
+  if (non_negative) global_bias = std::fmax(0., global_bias); /* <- should not happen, but just in case */
+
+  user_bias.zeros();
+  item_bias.zeros();
+
+  double sweight;
+  const double n_items = ConfCSR.n_rows;
+
+  for (int row = 0; row < ConfCSR.n_cols; row++) {
+    sweight = 0;
+    for (int ix = ConfCSR.col_ptrs[row]; ix < ConfCSR.col_ptrs[row + 1]; ix++) {
+      user_bias[row] += ConfCSR.values[ix] + global_bias * (1. - ConfCSR.values[ix]);
+      sweight += ConfCSR.values[ix] - 1.;
+    }
+    user_bias[row] -= global_bias * n_items;
+    user_bias[row] /= sweight + n_items + lambda;
+    user_bias[row] /= 3; /* <- item biases are unaccounted for, don't want to assign everything to the user */
+    if (non_negative) user_bias[row] = std::fmax(0., user_bias[row]);
+  }
+
+  const double n_users = ConfCSC.n_rows;
+  for (int col = 0; col < ConfCSC.n_cols; col++) {
+    sweight = 0;
+    for (int ix = ConfCSC.col_ptrs[col]; ix < ConfCSC.col_ptrs[col + 1]; ix++) {
+      item_bias[col] += ConfCSC.values[ix] + global_bias * (1. - ConfCSC.values[ix]);
+      sweight += ConfCSC.values[ix] - 1.;
+    }
+    item_bias[col] -= global_bias * n_users;
+    item_bias[col] /= sweight + n_users + lambda;
+    item_bias[col] /= 3; /* <- user biases are unaccounted for */
+    if (non_negative) item_bias[col] = std::fmax(0., item_bias[col]);
+  }
+
+  return global_bias;
+}
+
+
+template <class T>
+double initialize_biases(dMappedCSC& ConfCSC,  // modified in place
+                         dMappedCSC& ConfCSR,  // modified in place
+                         arma::Col<T>& user_bias, arma::Col<T>& item_bias, T lambda,
+                         bool dynamic_lambda, bool non_negative,
+                         bool calculate_global_bias, bool is_explicit_feedback,
+                         const bool initialize_item_biases) {
+  if (is_explicit_feedback)
+    return initialize_biases_explicit(ConfCSC, ConfCSR, user_bias, item_bias,
+                                      lambda, dynamic_lambda, non_negative,
+                                      calculate_global_bias);
+  else
+    return initialize_biases_implicit(ConfCSC, ConfCSR, user_bias, item_bias, lambda,
+                                      calculate_global_bias,non_negative, initialize_item_biases);
+}
diff --git a/man/FTRL.Rd b/man/FTRL.Rd
index bad1c20..81f472a 100644
--- a/man/FTRL.Rd
+++ b/man/FTRL.Rd
@@ -16,8 +16,7 @@ y = sample(c(0, 1), 1000, TRUE)
 x = sample(c(-1, 1), 1000 * 100, TRUE)
 odd = seq(1, 99, 2)
 x[i \%in\% which(y == 1) & j \%in\% odd] = 1
-m = sparseMatrix(i = i, j = j, x = x, dims = c(1000, 1000), giveCsparse = FALSE)
-x = as(m, "RsparseMatrix")
+x = sparseMatrix(i = i, j = j, x = x, dims = c(1000, 1000), repr="R")
 
 ftrl = FTRL$new(learning_rate = 0.01, learning_rate_decay = 0.1,
 lambda = 10, l1_ratio = 1, dropout = 0)
@@ -26,7 +25,7 @@ ftrl$partial_fit(x, y)
 w = ftrl$coef()
 head(w)
 sum(w != 0)
-p = ftrl$predict(m)
+p = ftrl$predict(x)
 }
 \section{Methods}{
 \subsection{Public methods}{
diff --git a/man/MatrixFactorizationRecommender.Rd b/man/MatrixFactorizationRecommender.Rd
index 7b5d1af..d74e344 100644
--- a/man/MatrixFactorizationRecommender.Rd
+++ b/man/MatrixFactorizationRecommender.Rd
@@ -13,6 +13,9 @@ All matrix factorization recommenders inherit from this class
 \item{\code{components}}{item embeddings}
 
 \item{\code{global_bias}}{global mean (for centering values in explicit feedback)}
+
+\item{\code{global_bias_base}}{Pre-calculated `-(factors*global_bias)` (for centering values in
+implicit feedback when not using user/item biases)}
 }
 \if{html}{\out{</div>}}
 }
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index ffc4b0f..ae09a0a 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -279,14 +279,14 @@ BEGIN_RCPP
 END_RCPP
 }
 // c_nnls_double
-arma::Mat<double> c_nnls_double(const arma::mat& x, const arma::vec& y, uint max_iter, double rel_tol);
+arma::Mat<double> c_nnls_double(const arma::mat& x, const arma::vec& y, unsigned int max_iter, double rel_tol);
 RcppExport SEXP _rsparse_c_nnls_double(SEXP xSEXP, SEXP ySEXP, SEXP max_iterSEXP, SEXP rel_tolSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP);
     Rcpp::traits::input_parameter< const arma::vec& >::type y(ySEXP);
-    Rcpp::traits::input_parameter< uint >::type max_iter(max_iterSEXP);
+    Rcpp::traits::input_parameter< unsigned int >::type max_iter(max_iterSEXP);
     Rcpp::traits::input_parameter< double >::type rel_tol(rel_tolSEXP);
     rcpp_result_gen = Rcpp::wrap(c_nnls_double(x, y, max_iter, rel_tol));
     return rcpp_result_gen;
@@ -445,8 +445,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, 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) {
+double als_implicit_double(const Rcpp::S4& m_csc_r, arma::mat& X, arma::mat& Y, const arma::mat& XtX, double lambda, int n_threads, const unsigned int solver, const unsigned int cg_steps, const bool with_biases, const bool is_x_bias_last_row, const double global_bias, arma::vec& global_bias_base, const bool initialize_bias_base);
+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, SEXP global_biasSEXP, SEXP global_bias_baseSEXP, SEXP initialize_bias_baseSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -455,18 +455,21 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< arma::mat& >::type Y(YSEXP);
     Rcpp::traits::input_parameter< const arma::mat& >::type XtX(XtXSEXP);
     Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP);
-    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< int >::type n_threads(n_threadsSEXP);
+    Rcpp::traits::input_parameter< const unsigned int >::type solver(solverSEXP);
+    Rcpp::traits::input_parameter< const unsigned int >::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, with_biases, is_x_bias_last_row));
+    Rcpp::traits::input_parameter< const bool >::type is_x_bias_last_row(is_x_bias_last_rowSEXP);
+    Rcpp::traits::input_parameter< const double >::type global_bias(global_biasSEXP);
+    Rcpp::traits::input_parameter< arma::vec& >::type global_bias_base(global_bias_baseSEXP);
+    Rcpp::traits::input_parameter< const bool >::type initialize_bias_base(initialize_bias_baseSEXP);
+    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, global_bias, global_bias_base, initialize_bias_base));
     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, 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) {
+double als_implicit_float(const Rcpp::S4& m_csc_r, Rcpp::S4& X_, Rcpp::S4& Y_, Rcpp::S4& XtX_, double lambda, int n_threads, const unsigned int solver, const unsigned int cg_steps, const bool with_biases, const bool is_x_bias_last_row, const double global_bias, Rcpp::S4& global_bias_base_, const bool initialize_bias_base);
+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, SEXP global_biasSEXP, SEXP global_bias_base_SEXP, SEXP initialize_bias_baseSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -475,18 +478,21 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< Rcpp::S4& >::type Y_(Y_SEXP);
     Rcpp::traits::input_parameter< Rcpp::S4& >::type XtX_(XtX_SEXP);
     Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP);
-    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< int >::type n_threads(n_threadsSEXP);
+    Rcpp::traits::input_parameter< const unsigned int >::type solver(solverSEXP);
+    Rcpp::traits::input_parameter< const unsigned int >::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, with_biases, is_x_bias_last_row));
+    Rcpp::traits::input_parameter< const bool >::type is_x_bias_last_row(is_x_bias_last_rowSEXP);
+    Rcpp::traits::input_parameter< const double >::type global_bias(global_biasSEXP);
+    Rcpp::traits::input_parameter< Rcpp::S4& >::type global_bias_base_(global_bias_base_SEXP);
+    Rcpp::traits::input_parameter< const bool >::type initialize_bias_base(initialize_bias_baseSEXP);
+    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, global_bias, global_bias_base_, initialize_bias_base));
     return rcpp_result_gen;
 END_RCPP
 }
 // initialize_biases_double
-double initialize_biases_double(const Rcpp::S4& m_csc_r, const Rcpp::S4& m_csr_r, arma::Col<double>& user_bias, arma::Col<double>& item_bias, double lambda, bool dynamic_lambda, bool non_negative, bool calculate_global_bias);
-RcppExport SEXP _rsparse_initialize_biases_double(SEXP m_csc_rSEXP, SEXP m_csr_rSEXP, SEXP user_biasSEXP, SEXP item_biasSEXP, SEXP lambdaSEXP, SEXP dynamic_lambdaSEXP, SEXP non_negativeSEXP, SEXP calculate_global_biasSEXP) {
+double initialize_biases_double(const Rcpp::S4& m_csc_r, const Rcpp::S4& m_csr_r, arma::Col<double>& user_bias, arma::Col<double>& item_bias, double lambda, bool dynamic_lambda, bool non_negative, bool calculate_global_bias, bool is_explicit_feedback, const bool initialize_item_biases);
+RcppExport SEXP _rsparse_initialize_biases_double(SEXP m_csc_rSEXP, SEXP m_csr_rSEXP, SEXP user_biasSEXP, SEXP item_biasSEXP, SEXP lambdaSEXP, SEXP dynamic_lambdaSEXP, SEXP non_negativeSEXP, SEXP calculate_global_biasSEXP, SEXP is_explicit_feedbackSEXP, SEXP initialize_item_biasesSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -498,13 +504,15 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< bool >::type dynamic_lambda(dynamic_lambdaSEXP);
     Rcpp::traits::input_parameter< bool >::type non_negative(non_negativeSEXP);
     Rcpp::traits::input_parameter< bool >::type calculate_global_bias(calculate_global_biasSEXP);
-    rcpp_result_gen = Rcpp::wrap(initialize_biases_double(m_csc_r, m_csr_r, user_bias, item_bias, lambda, dynamic_lambda, non_negative, calculate_global_bias));
+    Rcpp::traits::input_parameter< bool >::type is_explicit_feedback(is_explicit_feedbackSEXP);
+    Rcpp::traits::input_parameter< const bool >::type initialize_item_biases(initialize_item_biasesSEXP);
+    rcpp_result_gen = Rcpp::wrap(initialize_biases_double(m_csc_r, m_csr_r, user_bias, item_bias, lambda, dynamic_lambda, non_negative, calculate_global_bias, is_explicit_feedback, initialize_item_biases));
     return rcpp_result_gen;
 END_RCPP
 }
 // initialize_biases_float
-double initialize_biases_float(const Rcpp::S4& m_csc_r, const Rcpp::S4& m_csr_r, Rcpp::S4& user_bias, Rcpp::S4& item_bias, double lambda, bool dynamic_lambda, bool non_negative, bool calculate_global_bias);
-RcppExport SEXP _rsparse_initialize_biases_float(SEXP m_csc_rSEXP, SEXP m_csr_rSEXP, SEXP user_biasSEXP, SEXP item_biasSEXP, SEXP lambdaSEXP, SEXP dynamic_lambdaSEXP, SEXP non_negativeSEXP, SEXP calculate_global_biasSEXP) {
+double initialize_biases_float(const Rcpp::S4& m_csc_r, const Rcpp::S4& m_csr_r, Rcpp::S4& user_bias, Rcpp::S4& item_bias, double lambda, bool dynamic_lambda, bool non_negative, bool calculate_global_bias, bool is_explicit_feedback, const bool initialize_item_biases);
+RcppExport SEXP _rsparse_initialize_biases_float(SEXP m_csc_rSEXP, SEXP m_csr_rSEXP, SEXP user_biasSEXP, SEXP item_biasSEXP, SEXP lambdaSEXP, SEXP dynamic_lambdaSEXP, SEXP non_negativeSEXP, SEXP calculate_global_biasSEXP, SEXP is_explicit_feedbackSEXP, SEXP initialize_item_biasesSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -516,7 +524,9 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< bool >::type dynamic_lambda(dynamic_lambdaSEXP);
     Rcpp::traits::input_parameter< bool >::type non_negative(non_negativeSEXP);
     Rcpp::traits::input_parameter< bool >::type calculate_global_bias(calculate_global_biasSEXP);
-    rcpp_result_gen = Rcpp::wrap(initialize_biases_float(m_csc_r, m_csr_r, user_bias, item_bias, lambda, dynamic_lambda, non_negative, calculate_global_bias));
+    Rcpp::traits::input_parameter< bool >::type is_explicit_feedback(is_explicit_feedbackSEXP);
+    Rcpp::traits::input_parameter< const bool >::type initialize_item_biases(initialize_item_biasesSEXP);
+    rcpp_result_gen = Rcpp::wrap(initialize_biases_float(m_csc_r, m_csr_r, user_bias, item_bias, lambda, dynamic_lambda, non_negative, calculate_global_bias, is_explicit_feedback, initialize_item_biases));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -551,10 +561,10 @@ static const R_CallMethodDef CallEntries[] = {
     {"_rsparse_deep_copy", (DL_FUNC) &_rsparse_deep_copy, 1},
     {"_rsparse_als_explicit_double", (DL_FUNC) &_rsparse_als_explicit_double, 11},
     {"_rsparse_als_explicit_float", (DL_FUNC) &_rsparse_als_explicit_float, 11},
-    {"_rsparse_als_implicit_double", (DL_FUNC) &_rsparse_als_implicit_double, 10},
-    {"_rsparse_als_implicit_float", (DL_FUNC) &_rsparse_als_implicit_float, 10},
-    {"_rsparse_initialize_biases_double", (DL_FUNC) &_rsparse_initialize_biases_double, 8},
-    {"_rsparse_initialize_biases_float", (DL_FUNC) &_rsparse_initialize_biases_float, 8},
+    {"_rsparse_als_implicit_double", (DL_FUNC) &_rsparse_als_implicit_double, 13},
+    {"_rsparse_als_implicit_float", (DL_FUNC) &_rsparse_als_implicit_float, 13},
+    {"_rsparse_initialize_biases_double", (DL_FUNC) &_rsparse_initialize_biases_double, 10},
+    {"_rsparse_initialize_biases_float", (DL_FUNC) &_rsparse_initialize_biases_float, 10},
     {NULL, NULL, 0}
 };
 
diff --git a/src/rsparse.h b/src/rsparse.h
index 7a35051..c562e72 100644
--- a/src/rsparse.h
+++ b/src/rsparse.h
@@ -3,6 +3,8 @@
 #include "mapped_csr.hpp"
 
 #include <stdexcept>
+#include <limits>
+#include <cmath>
 
 #ifdef _OPENMP
 #include <omp.h>
diff --git a/src/wrmf_implicit.cpp b/src/wrmf_implicit.cpp
index 84843fb..3defc3c 100644
--- a/src/wrmf_implicit.cpp
+++ b/src/wrmf_implicit.cpp
@@ -3,24 +3,29 @@
 
 // [[Rcpp::export]]
 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) {
+                           const arma::mat& XtX, double lambda, int n_threads,
+                           const unsigned int solver, const unsigned int cg_steps, const bool with_biases,
+                           const bool is_x_bias_last_row, const double global_bias,
+                           arma::vec& global_bias_base, const bool initialize_bias_base) {
   const dMappedCSC Conf = extract_mapped_csc(m_csc_r);
   return (double)als_implicit<double>(Conf, X, Y, XtX, lambda, n_threads, solver,
-                                      cg_steps, with_biases, is_x_bias_last_row);
+                                      cg_steps, with_biases, is_x_bias_last_row,
+                                      global_bias, global_bias_base, initialize_bias_base);
 }
 
 // [[Rcpp::export]]
 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) {
+                          Rcpp::S4& XtX_, double lambda, int n_threads,
+                          const unsigned int solver, const unsigned int cg_steps, const bool with_biases,
+                          const bool is_x_bias_last_row, const double global_bias,
+                          Rcpp::S4& global_bias_base_, const bool initialize_bias_base) {
   const dMappedCSC Conf = extract_mapped_csc(m_csc_r);
   // get arma matrices which share memory with R "float" matrices
   arma::fmat X = extract_float_matrix(X_);
   arma::fmat Y = extract_float_matrix(Y_);
   arma::fmat XtX = extract_float_matrix(XtX_);
+  arma::fvec global_bias_base = extract_float_vector(global_bias_base_);
   return (double)als_implicit<float>(Conf, X, Y, XtX, lambda, n_threads, solver, cg_steps,
-                                     with_biases, is_x_bias_last_row);
+                                     with_biases, is_x_bias_last_row,
+                                     global_bias, global_bias_base, initialize_bias_base);
 }
diff --git a/src/wrmf_init.cpp b/src/wrmf_init.cpp
index 686e7c3..7310746 100644
--- a/src/wrmf_init.cpp
+++ b/src/wrmf_init.cpp
@@ -7,18 +7,23 @@ double initialize_biases_double(const Rcpp::S4& m_csc_r, const Rcpp::S4& m_csr_r
                                 arma::Col<double>& user_bias,
                                 arma::Col<double>& item_bias, double lambda,
                                 bool dynamic_lambda, bool non_negative,
-                                bool calculate_global_bias = false) {
+                                bool calculate_global_bias = false,
+                                bool is_explicit_feedback = false,
+                                const bool initialize_item_biases = false) {
   dMappedCSC ConfCSC = extract_mapped_csc(m_csc_r);
   dMappedCSC ConfCSR = extract_mapped_csc(m_csr_r);
   return initialize_biases<double>(ConfCSC, ConfCSR, user_bias, item_bias, lambda,
-                                   dynamic_lambda, non_negative, calculate_global_bias);
+                                   dynamic_lambda, non_negative, calculate_global_bias,
+                                   is_explicit_feedback, initialize_item_biases);
 }
 
 // [[Rcpp::export]]
 double initialize_biases_float(const Rcpp::S4& m_csc_r, const Rcpp::S4& m_csr_r,
                                Rcpp::S4& user_bias, Rcpp::S4& item_bias, double lambda,
                                bool dynamic_lambda, bool non_negative,
-                               bool calculate_global_bias = false) {
+                               bool calculate_global_bias = false,
+                               bool is_explicit_feedback = false,
+                                const bool initialize_item_biases = false) {
   dMappedCSC ConfCSC = extract_mapped_csc(m_csc_r);
   dMappedCSC ConfCSR = extract_mapped_csc(m_csr_r);
 
@@ -27,5 +32,6 @@ double initialize_biases_float(const Rcpp::S4& m_csc_r, const Rcpp::S4& m_csr_r,
 
   return initialize_biases<float>(ConfCSC, ConfCSR, user_bias_arma, item_bias_arma,
                                   lambda, dynamic_lambda, non_negative,
-                                  calculate_global_bias);
+                                  calculate_global_bias, is_explicit_feedback,
+                                  initialize_item_biases);
 }