4
\$\begingroup\$

I am trying in R to "center" a matrix (i.e. remove mean) group-wise (i.e. remove group means for each variable).

Input is a matrix (including the grouping variable), and output is the same matrix, where columns are now group-wise centered (possibly excluding now grouping variable).

I tried 3 solutions:

  1. Using dplyr: use group_by(cell) and mutate_all(funs(. - mean(.)))
  2. Using data.table: dt[, lapply(.SD, function(x) x - mean(x)), by= cell]
  3. Using data.table: same as above, but in 2 steps: instead of over-writing variables, 1) add mean variables as new columns, 2) then compute differences, splitting the original matrix as two different matrices.

Results in terms of speed are: 3 < 2 < 1. (3) being faster than (2) is surprising (as (2) is just overwriting variables, not adding)... I suspect it could be from the fact that using function(x) x - mean(x) prevents data.table to use an optimised version of the mean function.

My questions are:

  1. Any way I can speed up some of the code?
  2. How can I understand (3) is faster than (2)? Can I get (2) to be faster?

enter image description here

library(tidyverse)
library(data.table)

## function (1)
center_dplyr <- function(x) {
  x %>%
    group_by(cell) %>%
    mutate_all(funs(. - mean(.))) %>%
    ungroup() %>%
    dplyr::select(-cell)
}

## function (2)
center_dt_1 <- function(x) {
  x <- as.data.table(x)
  setkey(x, cell)
  res <- x[, lapply(.SD, function(x) x - mean(x)), by= cell][, -"cell"]
  as.data.frame(res)
}

## function (3)
center_dt_2 <- function(x) {
  x <- as.data.table(x)
  x_names <- colnames(x)[colnames(x) !="cell"]
  x_names_new <-  paste(x_names, "mean", sep="_")

  setkey(x, cell)
  x[,  paste(x_names, "mean", sep="_"):= lapply(.SD, mean, na.rm = TRUE), by = cell]
  res <- x[, x_names, with=FALSE] - x[, x_names_new, with=FALSE]
  as.data.frame(res)
}

## Data
T = 6; 
N = 10^4
set.seed(123)
sim_df <- data.frame(A = sample(c(0,1), N * T, replace = TRUE),
                      B1 = sample(c(0,1), N * T, replace = TRUE),
                      B2 = rnorm(N),
                      cell = rep(1:N, each = T))


ans_dplyr <- center_dplyr(x=sim_df)
ans_dt1 <- center_dt_1(sim_df)
ans_dt2 <- center_dt_2(sim_df)

all.equal(ans_dt1, ans_dplyr, check.attributes = FALSE)
all.equal(ans_dt2, ans_dplyr, check.attributes = FALSE)


### Benchmark:

library(microbenchmark)

## small
sim_df_s <- sim_df[1:1000,]
bench_small <-  microbenchmark(ans_dplyr = center_dplyr(sim_df_s),
                           ans_dt1 = center_dt_1(sim_df_s),
                           ans_dt2 = center_dt_2(sim_df_s),
                           times = 10)


bench_large <- microbenchmark(ans_dplyr = center_dplyr(sim_df),
                              ans_dt1 = center_dt_1(sim_df),
                              ans_dt2 = center_dt_2(sim_df),
                              times = 10)

bench_all <- rbind(bench_small %>% summary %>% mutate(data_size = "small"),
                   bench_large %>% summary %>% mutate(data_size = "large")) %>%
  select(data_size)

bench_all
\$\endgroup\$
  • 1
    \$\begingroup\$ Potentially some performance could be gained by avoiding as.data.frame/as.data.table and using setDF/setDT as the former copy their inputs. That said, you might need the inputs to be copied and the performance difference seems pretty small. \$\endgroup\$ – Hugh Sep 11 '18 at 3:34
  • \$\begingroup\$ oh, thanks for the hint! i tried it and indeed it makes a small performance improvement. But I am curious about the data table 1 would be slower compared to the second? \$\endgroup\$ – Matifou Sep 11 '18 at 5:55
  • \$\begingroup\$ Pretty sure your diagnosis is correct: mean is optimized but not x - mean(x). You could check it by using a different function than mean \$\endgroup\$ – Hugh Sep 11 '18 at 6:05
  • \$\begingroup\$ You talk about centering a matrix, but your sample data looks to be in long form (if I understand your sample data correctly). I suspect you may get some speed gains via an algorithm that uses matrix math (and therefore matrix libraries). The other option (or in addition to using matrix libraries) would be to write your algorithm in C++ and access it via Rcpp. Given that the operation itself is fairly simple, it shouldn't be too hard to implement via Rcpp. \$\endgroup\$ – mikeck Sep 15 '18 at 17:57
2
\$\begingroup\$

This answer only goes to how to speed up processing and doesn't get at the intricacies of data.table vs dplyr under the hood.

Below is a slightly faster version (at least for larger datasets).

Rather than create additional columns in the main DT, I created a new DT that had the mean values for each. As that had a 6-fold smaller size, the second DT has to be keyed and then joined back to the original.

This produces a bunch of extra columns (with values from x), so that's why we only selection columns of [,1:length...]. Note that subtraction of two DTs ignores column labels and assumes the same order of all the columns - i.e without the setcolorder, mean_x has "cell" being the first column, which results in the wrong columns being subtracted.

center_dt_3 <- function(x) {
  setDT(x,key="cell")
  x_fullnames <- names(x)
  mean_x <- x[,lapply(.SD,mean, na.rm=TRUE),by=cell]
  setcolorder(mean_x,x_fullnames)
  setkey(mean_x,cell)
  res<-x-mean_x[x,][,1:length(names(mean_x))]
  res[,cell:=NULL]
  setDF(res)
  return(res)
}

Perhaps your package version of microbenchmark is different than mine

bench_all <- rbind(bench_small %>% summary %>% mutate(expr = paste0(expr,"_small")),
                   bench_large %>% summary %>% mutate(expr = paste0(expr,"_large")))

bench_all
             expr        min        lq       mean     median         uq        max neval
1 ans_dplyr_small  14.016470  14.63970  16.057286  14.919529  15.477332  26.772788    10
2   ans_dt1_small   2.033062   2.13456   2.349156   2.189441   2.638336   2.754370    10
3   ans_dt2_small   2.217810   2.27483   2.438935   2.407972   2.628640   2.661144    10
4   ans_dt3_small   2.925432   3.00355   3.896746   3.093642   3.514454  10.440436    10
5 ans_dplyr_large 660.563300 677.73756 709.718870 695.432125 707.791010 848.521906    10
6   ans_dt1_large  54.344144  55.17835  59.058532  56.304646  62.696206  69.227890    10
7   ans_dt2_large  21.906666  22.25135  23.184892  22.716924  23.423408  28.072566    10
8   ans_dt3_large   8.380004   8.94536   9.065873   9.136237   9.218774   9.665814    10

Doing this a bit more natively in data.table (unfortunately I couldn't figure out how to dereference in x and i, as with=FALSE assumes the name of each list is e.g. x.x_names)

center_dt_4 <- function(t) {
  setDT(t,key="cell")
  x_fullnames <- names(t)
  x_names <- x_fullnames[which(x_fullnames != "cell")]
  mean_x <- t[,lapply(.SD,mean,na.rm=TRUE),by=cell]
  res<-t[mean_x,.(A=x.A-i.A,B1=x.B1-i.B1,B2=x.B2-i.B2),on="cell"]
  setDF(res)
  return(res)
}  
> bench_all
                        expr        min         lq       mean     median         uq        max neval
1            ans_dplyr_small  14.366576  14.494300  15.392998  14.925517  15.385244  20.094272    10
2              ans_dt1_small   2.029644   2.216670   2.384766   2.381316   2.554228   2.729282    10
3  ans_dt1_precompiled_small   3.341396   3.495352   3.804031   3.781878   4.003972   4.439324    10
4              ans_dt2_small   2.155372   2.321018   2.463510   2.405263   2.584166   2.936550    10
5              ans_dt3_small   2.638620   2.760072   2.859517   2.873259   2.965346   3.063422    10
6              ans_dt4_small   1.925294   1.930712   2.961612   2.119307   2.510894  10.003090    10
7            ans_dplyr_large 672.502800 677.593296 692.047965 685.753199 706.827080 716.847560    10
8              ans_dt1_large  55.923322  56.490104  59.640225  60.677824  62.484374  63.247306    10
9  ans_dt1_precompiled_large 135.373544 139.275732 141.147599 141.397891 142.505084 150.496206    10
10             ans_dt2_large  21.825126  22.328332  23.785973  22.710226  23.469594  29.990162    10
11             ans_dt3_large   8.644578   8.912002   9.775805   9.052987   9.578002  15.339058    10
12             ans_dt4_large   7.787560   8.263110   8.517451   8.644150   8.788840   8.970734    10
\$\endgroup\$
  • \$\begingroup\$ oh thanks! That's cool, a third way to do it. Actually, I think the merge might not be necessary. one could just use x's index? mean_x <- x[, lapply(.SD,mean, na.rm=TRUE), by=cell][,cell:=NULL][x$cell,] This seems to be faster :-) Thanks!! \$\endgroup\$ – Matifou Sep 15 '18 at 2:49
  • \$\begingroup\$ i think there's no guarantee with that method that you get the proper mapping of mean_x to x. As cell is the only thing in common between the two tables, your elimination of it with [,cell:=NULL] makes it so you're binding a 10000 item to a 60000 by repeating the 10000 element end-to-end (10000*6), while I think you want a repeat of each item 6 times before moving on to the next (6*10000). \$\endgroup\$ – mpag Sep 15 '18 at 14:34
  • \$\begingroup\$ I think it does! It's a little tricky, as it looks like a joint, but I think it does simply basic R indexing, and the by argument of data.table preserves initial ordering, so we are just asking to repeat each value 6 times, as int he raw dataset. It seems to be unofrmly faster, see the full simulations in github.com/MatthieuStigler/Misconometrics/blob/master/…, where 3a is your first aog, 3b your second, and 3c the same with indexing instead of joining \$\endgroup\$ – Matifou Sep 15 '18 at 23:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.