在列表中随机子集向量的任何更快的方法？

Question

我正在为嵌套在列表中的随机子集向量寻找快速解决方案。

如果我们模拟以下数据，我们会得到一个列表l，其中包含 300 万个向量，每个向量的长度为 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. 是否有更快的替代方法：

   rrapply(object = l, f = randomly_subset_vec)
   

2. 或者更一般地说，有没有更快的方法从my_named_vec 开始并到达l_subsetted？

Answer 1

更新 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 在 cmets 中是正确的，上面的代码将导致二项式分布的向量长度，而 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)

Answer 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"

Answer 3

简化采样函数：

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))
}

Answer 4

似乎最大的瓶颈是运行所有sample 调用，所以我们可以尝试以下方法。一种方法是solution by 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)

然后，使用funFast 而不是sample 定义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();

Answer 5

也许我们可以用sample 和sample.int 替换randomly_subset_vec：

lapply(l, function(x) x[sample.int(5, sample(5, 1))])

Answer 6

更有效的可能是将许多单独的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 的列表元素也不能保证至少包含一个值。

Answer 7

我将其放入新答案中，以免进一步混淆我之前的答案。

我从一些 cmets 中注意到，l 中的向量旨在具有所有相同的长度 (5)，您可能根本不需要 l。您是否希望 l_subsetted 的长度介于 0 和 4 之间或介于 0 和 5 之间也有点不清楚。您似乎也对 l_subsetted 的长度分布感兴趣（均匀与二项式）。

下面是一个通用函数 if 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

Answer 8

你可以试试下面的代码

lapply(
  l,
  function(x) {
    head(sample(x), sample(length(x), 1))
  }
)