futuremice scales poorly on 64/128 cores machines

[vkhodygo]

Describe the bug
futuremice runs just fine when I do 10 or 20 imputations at the same time. When I increase the number to say 50 or 100 while keeping all other parameters the same, it just sits there indefinitely.

To Reproduce
Not sure that any code would be suitable here.

Expected behavior
Running 10, 20, 50, 100 imputations at once should take roughly the same amount of time.

[stefvanbuuren]

This sounds more like a resource problem than a bug.

[vkhodygo]

@stefvanbuuren How is that so?

[stefvanbuuren]

It would be useful if we can have a reprex somehow, otherwise it is very hard for us to chase. Did you try setting m = 200 with a small problem?

[vkhodygo]

I can give it a go, but that might take some time.

[vkhodygo]

@stefvanbuuren I tried to come up with an example that closely matches my data, but that's cumbersome. Instead, I used the code from the futuremice vignette:

library(mice)
version()

set.seed(123)

n_features = 20
small_covmat <- diag(n_features)
small_covmat[small_covmat == 0] <- 0.5
small_data <- MASS::mvrnorm(10000, 
                      mu = c(1:n_features) * 0,
                      Sigma = small_covmat)

small_data_with_missings <- ampute(small_data, prop = 0.8, mech = "MCAR")$amp

n_streams <- 5
start_time <- Sys.time()

imp <- futuremice(small_data_with_missings,
                  parallelseed = 123,
                  n.core = n_streams,
                  m = n_streams,
                  maxit = 1,
                  method = "rf",
                  ntrees=10)

end_time <- Sys.time()
end_time - start_time

$ Rscript main.R

Attaching package: ‘mice’

The following object is masked from ‘package:stats’:

    filter

The following objects are masked from ‘package:base’:

    cbind, rbind

[1] "mice 3.16.0 2023-05-24 /home/software/.local/easybuild/software/R/4.2.0-foss-2021b/lib/R/library"
Time difference of 24.37422 secs

This is what I get when n_streams == 50:

$ Rscript main.R

Attaching package: ‘mice’

The following object is masked from ‘package:stats’:

    filter

The following objects are masked from ‘package:base’:

    cbind, rbind

[1] "mice 3.16.0 2023-05-24 /home/software/.local/easybuild/software/R/4.2.0-foss-2021b/lib/R/library"
Time difference of 2.550763 mins

and when n_features = 100:

$ Rscript main.R

Attaching package: ‘mice’

The following object is masked from ‘package:stats’:

    filter

The following objects are masked from ‘package:base’:

    cbind, rbind

[1] "mice 3.16.0 2023-05-24 /home/software/.local/easybuild/software/R/4.2.0-foss-2021b/lib/R/library"
Time difference of 2.003762 mins

Real-life numbers are much-much worse since I work mostly with categories, and their number is significantly higher. As this number goes up, literally every process starts spawning threads like there is no tomorrow. I understand that there is some overhead, but that's a bit too much. This automatically results in 100% load even when n_streams is low.

Just to show what I have to deal with: the same code with the actual data and n_streams == 5 needs about 25 minutes to finish the very first iteration on my laptop. On the cluster with two CPUs, 32 cores each, the code with n_streams == 25 has been running for an hour, and I have no idea when it'll be done.

[vkhodygo]

@stefvanbuuren Got it done, at least something:
The same code with n_streams == 5 on the cluster:

Time difference of 2.386232 hours

and with n_streams == 25:

Time difference of 12.78157 hours

I'd blame Intel MKL or something futureverse/future#405 , but those are AMD machines.

[stefvanbuuren]

Real-life numbers are much-much worse since I work mostly with categories, and their number is significantly higher.

Perhaps the problem is not with futuremice() but caused by a large number of categories. If you have 1001 categories, then mice tries to create a 1000 dummy variables...

What happens if you specify method = "pmm"?

[vkhodygo]

caused by a large number of categories. If you have 1001 categories, then mice tries to create a 1000 dummy variables...

The number of categories is high, that's true. However, this should not affect parallel and independent imputations.
I can run a handful of them just fine, but any further than that and it takes too much time.

Anyway, real data with pmm:

Error in (function (.x, .f, ..., .progress = FALSE)  : 
  ℹ In index: 1.
Caused by error in `chol.default()`:
! the leading minor of order 1 is not positive
Calls: futuremice ... resolve.list -> signalConditionsASAP -> signalConditions
Execution halted

The reprex works just fine.

[stefvanbuuren]

Thanks. On my desktop with 9 free cores, I found your reprex executes uses 15 seconds (n = 5), 14 seconds (n = 9), 21 seconds (n = 18), 47 seconds (n = 50) and 1.37 minutes (n = 100). I think this is as it should be.

I am not sure what causes the leading minor of order 1 is not positive error, but I have seen this error appearing when there are a lot of collinear variables. mice() tries very hard to remove these during the iteration by means of the internal remove.lindep() function. This checking process is - however - inefficient and can sometimes take >99% of the processor time.

Random forests (method = "rf") is quite robust against collinear variables, so remove.lindep() may be an overkill for your setup. It is possible to bypass remove.lindep() by adapting your call to mice(..., eps = 0). Could you gives this a try on the real data?

#225 #306

[vkhodygo]

I think I'm onto something here.

I modified the code a bit:

library(mice)
version()

library(dplyr)

set.seed(123)
n_features = c(5, 10, 20, 30, 40, 50)
available_strategies <- c("multisession")
testing_methods <- c("pmm", "cart", "rf")

df <- data.frame(row.names=c("iterations", "streams", "features",
                             "strategy", "imputation_method", "time"))

for (i in c(1, 5, 10)) {
  for (n_streams in c(2, 4, 5, 8)) {
    for (f in n_features) {
      
      small_covmat <- diag(f)
      small_covmat[small_covmat == 0] <- 0.5
      small_data <- MASS::mvrnorm(10000, 
                                  mu = c(1:f) * 0,
                                  Sigma = small_covmat)
      
      small_data_with_missings <- ampute(small_data,
                                         prop = 0.8,
                                         mech = "MCAR")$amp
      
      for (s in available_strategies) {
        for (m1 in testing_methods) {
          
          start_time <- Sys.time()
          imp <- futuremice(small_data_with_missings,
                            parallelseed = 123,
                            n.core = n_streams,
                            m = n_streams,
                            maxit = i,
                            method = m1,
                            future.plan = s,
                            future.debug = TRUE
          )
          end_time <- Sys.time()
          t <- as.double(difftime(end_time,
                                  start_time,
                                  units="secs"))
                         
          row <- data.frame(iterations = i,
                            streams = n_streams,
                            features = f,
                            strategy = s,
                            imputation_method = m1,
                            time = t
                            )

          df <- rbind(df, row)
          
          cat("Iterations:", i,
              "; streams:", n_streams,
              "; features:", f,
              "; strategy:", s,
              "; imputation method:", m1,
              '; time:', t,
              "\n")
          
        }
      }
    } 
  }
}
write.csv(df, "benchmark.csv")

benchmark_cluster.csv
benchmark_laptop.csv

My laptop is about 2 times faster than a cluster node, and I didn't go beyond the limits (8 threads at most). One would expect the laptop to finish much earlier, but that's not the case here.

Also, even with a very limited number of streams the cluster gets fully congested, 100% load even when it's supposed to be a bit above 6%. That means something goes really wrong. Every run created a set of processes/forks, and every one of them in its turn created as many threads as it could.

After reading futureverse/future#405 #225 and #306 I tried running the script with export MKL_VERBOSE=1 for rf on the laptop since I'm only interested in this method. This produces a lot of output, the most relevant part here is NThr:4 by default. export MKL_NUM_THREADS=1 does actually affect this number, however, the behaviour stayed the same.

This is the point where it becomes interesting. For the sake of completeness I decided to try and run this on the cluster for all three methods and export OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 MKL_NUM_THREADS=1. rf behaved almost like before, however, pmm & cart started acting completely differently. Now, every session needs only 2 threads and uses one at all times which is what one would expect.

I feel like this should be mentioned somewhere since this is not an intended behaviour.

The question of rf still remains. It looks like ranger spawns its own workers, however, I'm not sure how to approach this without breaking the existing functionality. We could limit the code to 1 thread per future as in #572 and revert back to randomForest as the default option and pass extra parameters to ranger thanks to #563.

[vkhodygo]

New benchmarks which seems to be noticeably better for certain cases.

benchmark_cluster_new.csv

[stefvanbuuren]

It could also help to simplify your model, e.g., by using quickpred() to select the most important predictors for each variable that you want to impute.

[vkhodygo]

I planned to do a more detailed analysis later, however, I already do not impute complete features; that's about a third of of them.

[stefvanbuuren]

Converting to discussion because I do not believe this is caused by a bug.

[stefvanbuuren]

Thanks. I ran your simulation script on my desktop, and found all behaved well. The results were similar to your cluster simulations, apart from being about 2 times faster. pmm is about 2-3 as fast as rf, and about 10-15 times as fast as cart across all conditions. I also tried setting eps = 0 and pred = quickpred(data), but those settings had only marginal effects.

I feel like this should be mentioned somewhere since this is not an intended behaviour.

I have no idea why your laptop simulations show problematic behaviour. As the timing depends on the configuration, I do not think there is something in mice we could do to improve. Please let me know if you find code tweaks that are worthwhile to incorporate.

[vkhodygo]

apart from being about 2 times faster

That's because of higher CPU clocks, that's what I also observed.

laptop simulations show problematic behaviour

It is not my laptop only, this problem persists across different machines.

I managed to catch one of the moments ("rf") I describe in my previous answer.

What you're supposed to see is constant CPU load, instead you observe spikes. When you run multiple R sessions, each of them considers all available resources as its own only.
Let's say n_streams == 5, that's 5 R sessions. Now, every one of them uses BLAS one way or another, and spawns multiple threads to do linear algebra calculations. Suddenly, instead of 5 loaded cores you see 50! On a 64/128 machine that's still well within the limits, but when n_streams increases switching between threads creates too much overhead since the task is CPU-bound.

Using export resolves this issue for cases when BLAS libraries spawn too many threads. Frankly speaking, you could probably allow a few threads per one session, that depends on the problem. However, random forests are different because ranger creates extra workers on its own.

What bothers me here is not the performance itself, but the fact that multiple parts of the code compete for the same resources resulting in poor an unpredictable timings.

Please, let me know if that makes any sense to you.