在 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
缺点是这会失去并行工作者。好处是我们减少了内存减速,因为我们预先分配了很多我们需要的东西。
我正在实施详细 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
缺点是这会失去并行工作者。好处是我们减少了内存减速,因为我们预先分配了很多我们需要的东西。