在 R 中执行详尽搜索的最快方法是什么

What is the fastest way to perform an exhaustive search in R

我正在实施详细 here 的超大规模 Relieff 算法版本。

简单地说,超大规模地势将特征集 N 分成几个随机子集 Ns,其中 Ns << N。然后它计算子集中特征的 Relieff 权重 Ns。对于每个特征,最终权重将是在出现该特定特征的不同子集中分配的最高权重。

我有大约 100 个主题的大约 80000 个特征。我可以使用以下代码在合理的时间内(在 25 个内核上约 5 分钟 运行ning)计算 8000 个特征的 10000 个子集(为了更容易分析,它被缩小到 100 个特征):

library(tidyverse)
library(magrittr)
library(CORElearn)
library(doParallel)

#create fake data for example
fake_table <- matrix(rnorm(100*100), ncol = 100) %>%
  as_tibble()
outcome <- rnorm(100)
#create fake data for example

#VLSRelieff code
start_time <- Sys.time()
myCluster <- makeCluster(25, # number of cores to use
                         type = "FORK")
registerDoParallel(myCluster)
result <- foreach(x = seq(1,10000)) %dopar% {
  #set seed for results consistency among different run
  set.seed(x)
  #subsample the features table by extracting a subset of columns
  subset_index <- sample(seq(1,ncol(fake_table)),size = round(ncol(fake_table)*.01))
  subset_matrix <- fake_table[,subset_index]
  #assign the outcome as last column of the subset
  subset_matrix[,ncol(subset_matrix)+1] <- outcome
  #use the function attrEval from the CORElearn package to calculate the Relieff weights for the subset
  rf_weights <- attrEval(formula = ncol(subset_matrix), subset_matrix, estimator = "RReliefFequalK")
  #create a data frame with as many columns as features in the subset and only one row 
  #with the Relieff weigths
  rf_df <- rf_weights %>%
    unname() %>%
    matrix(., ncol = length(.), byrow = TRUE) %>%
    as_tibble() %>%
    set_colnames(., names(rf_weights))}
end_time <- Sys.time()
end_time - start_time

然而,上面的代码只完成了一半的工作:另一半是针对每个特征,进入不同重复的结果并找到获得的最大值。我已经设法编写了一个工作代码,但它非常慢(我让它 运行 在停止它之前 2 小时,尽管它在测试时使用更少的功能 - 同样,这里它被缩小到 100 个功能和应该 运行 在 ~7 秒内):

start_time <- Sys.time()
myCluster <- makeCluster(25, # number of cores to use
                         type = "FORK")
registerDoParallel(myCluster)
#get all features name
feat_names <- colnames(fake_table)
#initalize an empty vector of zeros, with the names of the features
feat_wegiths <- rep(0, length(feat_names))
names(feat_wegiths) <- feat_names
#loop in parallel on the features name, for each feature name
feat_weight_foreach <- foreach(feat = feat_names, .combine = 'c') %dopar% {
  #initalize the weight as 0
  current_weigth <- 0
  #for all element in result (i.e. repetitions of the subsampling process)
  for (el in 1:length(result)){
    #assign new weight accessing the table
    new_weigth <- result[[el]][[1,feat]]
    #skip is empty (i.e. the features is not present in the current subset)
    if(is_empty(new_weigth)){next}
    #if new weight is higher than current weight assign the value to current weight
    if (current_weigth < new_weigth){
      current_weigth <- new_weigth}}
  current_weigth
}
end_time <- Sys.time()
end_time - start_time

如果我正确理解了您要尝试做的事情,那么答案比您想象的要简单。

如果我错了请纠正我,但您正在尝试从每个特征的 attrEval 中获取最大值? 如果是这样,那么为什么不将所有结果绑定到一个数据帧(或 data.table),然后像这样获得每列的最大值:

allResults <- result %>% data.table::rbindlist(fill = TRUE)
apply(allResults, 2, max, na.rm=TRUE)

这不是最终答案,但我建议,因为这是一个数值问题,所以用 C++ 编写一个函数。这将显着提高速度,我猜是某个数量级。在我看来,将 R 用于这个非常具体的数字任务只是碰壁。

Rcpp for everyone第一章说:

第1章适合使用Rcpp的场景 R 在某些操作上很弱。如果你需要下面列出的操作,是时候考虑使用Rcpp了。

  • 后面的迭代依赖于前面的循环操作 迭代.
  • 访问vector/matrix的每个元素。
  • 循环内的循环函数调用。
  • 动态改变矢量的大小。
  • 需要高级数据结构和算法的操作。

Wickham 的 Advanced R 也有关于该主题的精彩章节。

这遵循@DS_UNI的想法,但不是绑定列表,而是从初始循环创建矩阵的方法。也就是说,一个 tibbles 列表会让我们做额外的工作。相反,我们拥有制作矩阵所需的一切:

library(tidyverse)
library(magrittr)
library(CORElearn)
library(doParallel)

nr = 50L
nc = 200L

## generate data
set.seed(123)
mat = matrix(rnorm(nr * nc), ncol = nc, dimnames = list(NULL, paste0('V', seq_len(nc))))
outcome = rnorm(nr) 

## constants for sampling
n_reps = nc
nc_sample_size = round(nc * 0.01)

## pre-allocate result
res = matrix(0, nrow = n_reps, ncol = ncol(mat), dimnames = dimnames(mat))

st = Sys.time()
for (i in seq_len(n_reps)) {
  set.seed(i) 
  
  ## similar way to do data simulations as OP
  sub_cols = sample(seq_len(nc), nc_sample_size)
  sub_mat = cbind(mat[, sub_cols], outcome)
  rf_weights = attrEval(formula = ncol(sub_mat), as.data.frame(sub_mat), estimator = 'RReliefFequalK')
  
  ## assign back to pre-allocated result
  res[i, sub_cols] = rf_weights
}

## get max of each column
apply(res, 2L, max)
et = Sys.time()
et - st

缺点是这会失去并行工作者。好处是我们减少了内存减速,因为我们预先分配了很多我们需要的东西。