有更快的方法随机选取列表中向量的子集吗?

7

我正在寻找一种快速的解决方案,用于随机地从嵌套在列表中的向量中进行子集抽样。

如果我们模拟以下数据,将得到一个包含 300 万个向量的列表l,每个向量的长度为 5。但是,我希望每个向量的长度各不相同。所以,我想应用一个函数来随机地对每个向量进行子集抽样。问题是,这种方法不如我所希望的那样快速。

模拟数据:列表l

library(stringi)

set.seed(123)
vec_n <- 15e6
vec_vals  <- 1:vec_n
vec_names <- stringi::stri_rand_strings(vec_n, 5)

my_named_vec <- setNames(vec_vals, vec_names)

split_func <- function(x, n) {
  unname(split(x, rep_len(1:n, length(x))))
}

l <- split_func(my_named_vec, n = vec_n / 5)

head(l)
#> [[1]]
#>    HmPsw    Qk8NP    Quo3T    8f0GH    nZmjN 
#>        1  3000001  6000001  9000001 12000001 
#> 
#> [[2]]
#>    2WtYS    ZaHFl    6YjId    jbGuA    tAG65 
#>        2  3000002  6000002  9000002 12000002 
#> 
#> [[3]]
#>    xSgZ6    jM5Uw    ujPOc    CTV5F    5JRT5 
#>        3  3000003  6000003  9000003 12000003 
#> 
#> [[4]]
#>    tF2Kx    r4ZCI    Ooklo    VOLHU    M6z6H 
#>        4  3000004  6000004  9000004 12000004 
#> 
#> [[5]]
#>    tgdze    w8d1B    FYERK    jlClo    NQfsF 
#>        5  3000005  6000005  9000005 12000005 
#> 
#> [[6]]
#>    hXaH9    gsY1u    CjBwC    Oqqty    dxJ4c 
#>        6  3000006  6000006  9000006 12000006

现在我们有了l,我希望可以随机地对每个向量进行子集操作:也就是说,被子集化的元素数量(每个向量)将是随机的。因此,一种选项是设置以下实用函数:

randomly_subset_vec <- function(x) {
  my_range <- 1:length(x)
  x[-sample(my_range, sample(my_range))]
}

lapply(head(l), randomly_subset_vec)
#> [[1]]
#>   Quo3T 
#> 6000001 
#> 
#> [[2]]
#>   6YjId   jbGuA 
#> 6000002 9000002 
#> 
#> [[3]]
#>   xSgZ6   jM5Uw   ujPOc   CTV5F 
#>       3 3000003 6000003 9000003 
#> 
#> [[4]]
#>   Ooklo 
#> 6000004 
#> 
#> [[5]]
#> named integer(0)
#> 
#> [[6]]
#>    CjBwC    Oqqty    dxJ4c 
#>  6000006  9000006 12000006

但是在整个l上运行这个过程需要很长时间。我尝试使用rrapply,它是一个用于处理列表的快速包,在我的机器上只需 "110" 秒。

library(rrapply)
library(tictoc)

tic()
l_subsetted <- rrapply(object = l, f = randomly_subset_vec)
toc()
#> 110.23 sec elapsed

我很乐意做以下任何一件事:
  1. Is there a speedier alternative to:
    rrapply(object = l, f = randomly_subset_vec)
    
  2. Or more generally, is there a speedier way to start with my_named_vec and arrive at l_subsetted?
8个回答

5

更新1:修复了对于大型对象在stack中名称的行为问题。

由于您的子集不包含完整的集合,因此首先从每个向量中随机删除一个元素,然后随机保留所有其他元素:

library(stringi)

set.seed(123)
vec_n <- 15e6
vec_vals  <- 1:vec_n
vec_names <- stringi::stri_rand_strings(vec_n, 5)

my_named_vec <- setNames(vec_vals, vec_names)

split_func <- function(x, n) {
  unname(split(x, rep_len(1:n, length(x))))
}

l <- split_func(my_named_vec, n = vec_n / 5)
system.time({
  lenl <- lengths(l)
  # use stack to unlist the list while keeping the originating list index for each value
  vec_names <- names(unlist(l))
  blnKeep <- replace(sample(c(FALSE, TRUE), length(vec_names), replace = TRUE), ceiling(runif(length(l))*lenl) + c(0, head(cumsum(lenl), -1)), FALSE)
  temp <- stack(setNames(l, seq_along(l)))[blnKeep,]
  # re-list
  l_subsetted <- unname(split(setNames(temp$values, vec_names[blnKeep]), temp$ind))
})
#>    user  system elapsed 
#>  22.999   0.936  23.934
head(l_subsetted)
#> [[1]]
#>    HmPsw    nZmjN 
#>        1 12000001 
#> 
#> [[2]]
#>   2WtYS   6YjId 
#>       2 6000002 
#> 
#> [[3]]
#>   xSgZ6   jM5Uw   ujPOc 
#>       3 3000003 6000003 
#> 
#> [[4]]
#>   tF2Kx   r4ZCI 
#>       4 3000004 
#> 
#> [[5]]
#>    FYERK    NQfsF 
#>  6000005 12000005 
#> 
#> [[6]]
#>   gsY1u 
#> 3000006
Created on 2021-11-01 by the reprex package (v2.0.0)

更新2:对于均匀分布长度的向量:

@runr在评论中指出,上面的代码将导致二项式分布的向量长度,而OP的原始代码会产生均匀分布的向量长度。以下是如何使用相同的思想来获得均匀分布的向量长度的示例代码。该代码更为复杂,但运行时间似乎要快一些(可能是由于规避了stack):

library(stringi)
set.seed(123)
vec_n <- 15e6
vec_vals  <- 1:vec_n
vec_names <- stringi::stri_rand_strings(vec_n, 5)
my_named_vec <- setNames(vec_vals, vec_names)
split_func <- function(x, n) {
  unname(split(x, rep_len(1:n, length(x))))
}
l <- split_func(my_named_vec, n = vec_n / 5)

system.time({
  idx <- seq_along(l)
  lenl <- lengths(l)
  ul <- unlist(l)
  # get a random number of elements to remove from each vector
  nRemove <- ceiling(runif(length(l))*lenl)
  nRemove2 <- nRemove
  blnNotEmpty <- nRemove != lenl # will the subset vector have any elements?
  blnKeep <- rep(TRUE, length(l))
  
  # loop until the predetermined number of elements have been removed from each vector
  while (length(nRemove)) {
    # remove a random element from vectors that have too many
    ul <- ul[-(ceiling(runif(length(idx))*lenl[idx]) + c(0, head(cumsum(lenl), -1))[idx])]
    lenl[idx] <- lenl[idx] - 1L # decrement the vector lengths
    blnKeep <- nRemove != 1
    idx <- idx[blnKeep]
    nRemove <- nRemove[blnKeep] - 1L # decrement the number of elements left to remove
  }
  
  l_subsetted <- rep(list(integer(0)), length(l))
  l_subsetted[blnNotEmpty] <- unname(split(ul, rep.int(seq_along(l), lenl)))
})
#>    user  system elapsed 
#>  18.396   0.935  19.332
head(l_subsetted)
#> [[1]]
#>   Qk8NP   Quo3T   8f0GH 
#> 3000001 6000001 9000001 
#> 
#> [[2]]
#> integer(0)
#> 
#> [[3]]
#>    xSgZ6    ujPOc    CTV5F    5JRT5 
#>        3  6000003  9000003 12000003 
#> 
#> [[4]]
#>   tF2Kx   Ooklo   VOLHU 
#>       4 6000004 9000004 
#> 
#> [[5]]
#>    tgdze    w8d1B    jlClo    NQfsF 
#>        5  3000005  9000005 12000005 
#> 
#> [[6]]
#>    gsY1u    CjBwC    Oqqty    dxJ4c 
#>  3000006  6000006  9000006 12000006
# check that vector lengths are uniformly-distributed (lengths of 0-4 are equally likely)
table(lengths(l_subsetted))
#> 
#>      0      1      2      3      4 
#> 599633 599041 601209 600648 599469
Created on 2021-11-02 by the reprex package (v2.0.1)

谢谢。这是一个非常好的解决方案。但请注意,您的“l_subsetted”不包括原始的字母数字向量名称。似乎它们在过程中消失了。 - Emman
很不幸,我的输出结果与你的不同。我甚至通过 reprex() 运行了代码,但仍然没有得到相同的结果。你能否也尝试使用 reprex() 运行一下? - Emman
这很奇怪。当 vec_n <- 15e3 时,我得到了正确的字母数字名称,但是当 vec_n <- 15e6 时,这些名称被数字替换了。除了 stringi 之外,我没有安装任何其他包。我的 R 版本是 4.1.1。Windows 10 PC。 - Emman
2
@Emman 观察这种方法中长度的不同分布(相对于原本预期的结果),以确定在你的情况下是否合理。例如,调用“l_subsetted%>% lapply(。,length)%>% do.call(c,。)%>% table”并观察呈钟形的直方图,其中中位数为“2”。另一方面,OP代码中的原始实验将生成均匀分布。这可能是预期实验设计的关键差异。 - runr
@runr 是正确的。请查看我的答案中的 UPDATE 2,以获取长度均匀分布的向量。 - jblood94
显示剩余6条评论

2
简化取样函数:
randomly_subset_vec_2 <- function(x) {
  my_range <- length(x)
  x[-sample(my_range, sample(my_range, 1))]
}

仅这一点就可以显著提高速度。
虽然我没有测试过,但根据问题描述,去除一些元素(在sample之前的减号)是为了保留其他元素。为什么不提取一些元素(没有减号),从而保留它们呢?

更简单更快:直接从x中进行抽样到目前为止是最快的方法。
randomly_subset_vec_3 <- function(x) {
  sample(x, sample(length(x), 1))
}

谢谢!randomly_subset_vec_2在我的机器上将处理时间从110秒缩短到了53秒。大约快了两倍。不确定我是否理解了您关于删除与提取的问题。是的,我同意它们是一样的东西。 - Emman
@Emman 我的意思是去掉减号,请看修改。第二个函数快了35%。 - Rui Barradas

2

这段代码很粗糙,我并不是特别自豪。我相信有更加优雅的方法,但在我的机器上运行只需要几秒钟。

> # Make some fake data
> out <- lapply(1:3000000, function(i){sample(LETTERS, 5, replace = FALSE)})
> out[1:5]
[[1]]
[1] "D" "H" "C" "Y" "V"

[[2]]
[1] "M" "E" "H" "G" "S"

[[3]]
[1] "R" "P" "O" "L" "M"

[[4]]
[1] "C" "U" "G" "Q" "X"

[[5]]
[1] "Q" "L" "W" "O" "V"

> # Create list with ids to sample
> id <- lapply(1:3000000, function(i){sample(1:5, sample(1:5, 1), replace = FALSE)})
> id[1:5]
[[1]]
[1] 2

[[2]]
[1] 2 3 4 1 5

[[3]]
[1] 4

[[4]]
[1] 5

[[5]]
[1] 1 2

> # Extract the ids from the original data using the id list.
> # Like I said I'm not particularly proud of this but it gets the job
> # done quick enough on my computer
> out <- lapply(1:3000000, function(i){out[[i]][id[[i]]]})
> out[1:5]
[[1]]
[1] "H"

[[2]]
[1] "E" "H" "G" "M" "S"

[[3]]
[1] "L"

[[4]]
[1] "X"

[[5]]
[1] "Q" "L"

谢谢。虽然更新 out 只需要几秒钟,但计算 id 却占用了大部分时间。因此,在我的机器上,除了最初创建 out 之外,你建议的整个代码大约需要55秒。所以比我的原始方法快了2倍。 - Emman
我现在想知道是否有一种方法可以先创建一个包含1-5范围内随机值的矩阵作为“id”,然后将其转换为列表。 - Emman

2

看起来最大的瓶颈是运行所有的sample调用,因此我们可以尝试以下方法。一种方法是使用Julius Vainora的解决方案。首先,我们通过Rcpp生成funFast

library(inline)
library(Rcpp)
src <- 
'
int num = as<int>(size), x = as<int>(n);
Rcpp::NumericVector vx = Rcpp::clone<Rcpp::NumericVector>(x);
Rcpp::NumericVector pr = Rcpp::clone<Rcpp::NumericVector>(prob);
Rcpp::NumericVector rnd = rexp(x) / pr;
for(int i= 0; i<vx.size(); ++i) vx[i] = i;
std::partial_sort(vx.begin(), vx.begin() + num, vx.end(), Comp(rnd));
vx = vx[seq(0, num - 1)] + 1;
return vx;
'
incl <- 
'
struct Comp{
  Comp(const Rcpp::NumericVector& v ) : _v(v) {}
  bool operator ()(int a, int b) { return _v[a] < _v[b]; }
  const Rcpp::NumericVector& _v;
};
'
funFast <- cxxfunction(signature(n = "Numeric", size = "integer", prob = "numeric"),
                       src, plugin = "Rcpp", include = incl)

接下来,使用sample的替代方法funFast来定义你的randomly_subset_vec

'randomly_subset_vec_2' <- function(x) {
  range <- length(x)
  probs <- rep(1/range, range)
  
  o <- funFast(range, size = funFast(range, size = 1, prob = probs), prob = probs)
  return(x[-o])
}

tic();obj <- rrapply(object = l, f = randomly_subset_vec_2);toc();

@Emman,你尝试过这种方法吗? - runr
抱歉,是的。在我的机器上运行需要36秒,所以现在它是最快的。然而,cxxfunction()有自己的开销,因此总体上你的解决方案在我的机器上需要45秒。我需要想出一个适当的基准测试方法,而不是使用tictoc - Emman
@Emman cxxfunction只编译C++代码。您可以每个会话编译一次,或者加载已经编译的代码。 - runr
@Emman 无论如何,似乎具体来说,采样函数在这里占用了大部分时间,而不是设计或我们进行子集处理的方式等。例如,看看 profvis({ tic();obj <- lapply(l,randomly_subset_vec_2);toc(); }) 这里 profvis::profvis 提供了最大瓶颈的火焰图。简化原始函数以使采样大小固定(而不是随机),代码会更快,因为它调用了较少的 funFast。当然,我已经假设您将并行化到所有 CPU 核心,而不是当前的 1 个核心? - runr
但是,样本大小的随机性是问题固有的特征。并行化是一个好主意!我之前没有考虑过这一点。 - Emman

1
也许我们可以用samplesample.int来替换randomly_subset_vec,让代码更简单易懂。
lapply(l, function(x) x[sample.int(5, sample(5, 1))])

1
我将这部分内容放到新的答案中,以免进一步混淆之前的回答。
从某些评论中我注意到,l向量的长度都应该相同(为5),而且你可能不需要l。另外有点不太清楚你想让l_subsetted的长度在0到4之间还是在0到5之间。你似乎对l_subsetted的长度分布(均匀分布还是二项式分布)也感兴趣。
如果length(unique(lengths(l))) == 1,下面是一个通用的函数。它从my_named_vec直接进行子集化,而不创建l。它运行时间普遍在5-13秒之间。
set.seed(123)
vec_n <- 15e6L
my_named_vec <- setNames(1:vec_n, stringi::stri_rand_strings(vec_n, 5))

fSub <- function(nv, vecLen = 5L, maxLen = 5L, unif = FALSE) {
  # subset each named vector from the list l (l is not generated):
  # l <- unname(split(nv, rep_len(seq(length(nv)/vecLen), length(nv))))
  # INPUTS:
  #  nv: named vector whose length is a multiple of vecLen
  #  vecLen: the length of the vectors in l
  #  maxLen: the maximum length of the subsetted vectors
  #  unif: FALSE = binomial subset vector lengths
  #        TRUE = uniform subset vector lengths
  # OUTPUT: a list of named vectors subset from l
  
  nrw <- length(nv)%/%vecLen # length of the output list
  # get all possible logical indices for sampling each vector in l
  mKeep <- as.matrix(expand.grid(rep(list(c(TRUE, FALSE)), vecLen)), ncol = vecLen)
  nKeep <- rowSums(mKeep)
  # remove logical indices that would result in vectors greater than maxLen
  blnKeep <- nKeep <= maxLen
  mKeep <- mKeep[blnKeep,]
  nKeep <- nKeep[blnKeep]
  
  if (unif) {
    # sample mKeep with non-uniform probability in order to get uniform lengths
    iKeep <- sample(length(nKeep), nrw, replace = TRUE, prob = 1/choose(vecLen, nKeep))
  } else {
    iKeep <- sample(length(nKeep), nrw, replace = TRUE)
  }
  
  blnKeep <- c(mKeep[iKeep,])
  l <- rep(list(integer(0L)), nrw)
  l[iKeep != length(nKeep)] <- unname(split(nv[blnKeep], rep(1:nrw, vecLen)[blnKeep]))
  return(l)
}

lbinom5 <- fSub(my_named_vec) # binomial vector lengths (0 to 5)
lunif5 <- fSub(my_named_vec, unif = TRUE) # uniform vector lengths (0 to 5)
lbinom4 <- fSub(my_named_vec, maxLen = 4L) # binomial vector lenghts (0 to 4)
lunif4 <- fSub(my_named_vec, maxLen = 4L, unif = TRUE) # uniform vector lengths (0 to 4)

> microbenchmark::microbenchmark(
+   lbinom5 = {lbinom5 <- fSub(my_named_vec)},
+   lunif5 = {lunif5 <- fSub(my_named_vec, unif = TRUE)},
+   lbinom4 = {lbinom4 <- fSub(my_named_vec, maxLen = 4L)},
+   lunif4 = {lunif4 <- fSub(my_named_vec, maxLen = 4L, unif = TRUE)},
+   times = 10)
Unit: seconds
    expr      min       lq     mean    median       uq      max neval
 lbinom5 5.974837 8.060281 9.192600  9.014967 10.15609 13.01182    10
  lunif5 5.240133 6.618115 9.688577 10.799230 11.44718 12.73518    10
 lbinom4 5.082508 6.497218 8.636434  8.656817 11.40678 11.81519    10
  lunif4 5.468311 6.639423 8.310269  7.919579 10.28546 11.28075    10

1
更高效的方法可能是用一个较大的sample调用替换多个单独的sample调用。下面是一种方法,对一个大的逻辑矩阵keep进行采样(因为l最初具有矩形格式),并仅保留keep评估为TRUE的条目:
system.time({
  keep <- matrix(sample(c(TRUE, FALSE), size = vec_n, replace = TRUE), nrow = 5, ncol = length(l))
  l1 <- lapply(seq_along(l), function(i) l[[i]][keep[, i]])
})

#>    user  system elapsed 
#>   8.667   0.448   9.114

head(l1)

#> [[1]]
#>   HmPsw   Quo3T   8f0GH 
#>       1 6000001 9000001 
#> 
#> [[2]]
#>   2WtYS   ZaHFl   6YjId 
#>       2 3000002 6000002 
#> 
#> [[3]]
#>    xSgZ6    jM5Uw    ujPOc    CTV5F    5JRT5 
#>        3  3000003  6000003  9000003 12000003 
#> 
#> [[4]]
#>    M6z6H 
#> 12000004 
#> 
#> [[5]]
#>    tgdze    w8d1B    FYERK    jlClo    NQfsF 
#>        5  3000005  6000005  9000005 12000005 
#> 
#> [[6]]
#>   hXaH9   CjBwC   Oqqty 
#>       6 6000006 9000006

注意:这里 l 中的条目顺序保持不变(即没有重新采样),并且 l1 中的列表元素不能保证至少包含一个值。


另外,请注意 table(lengths(l1)) 的输出。我们得到了一个非均匀分布,表明随机性具有一定的模式(因此不完全是随机的)。请参见 @runr 的评论 - Emman
2
@Emman,它们仍然是随机子集。这将给出二项式分布的向量长度(其中n = 5)。这只取决于您希望如何采样。此答案以0.5的概率随机保留/删除每个元素,而原始帖子则随机地为每个向量抽取均匀分布的元素数量。 - jblood94

0
你可以尝试下面的代码。
lapply(
  l,
  function(x) {
    head(sample(x), sample(length(x), 1))
  }
)

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接