在 R 中矢量化这些嵌套的 for 循环答案

【问题标题】：Vectorize these nested for loops in R在 R 中矢量化这些嵌套的 for 循环
【发布时间】：2023-03-21 08:24:01
【问题描述】：

我通常只要稍加思考就能弄清楚如何进行矢量化，但是尽管阅读了一堆 StackOverflow 问答，我仍然感到困惑！我想用合适的 apply 函数替换这些嵌套的 for 循环，但如果有一些明显不同的方法可以解决我错过的整个问题，请随时告诉我！

在测试的上下文中考虑这个示例，其中第一行是关键，随后的每一行都是学生的答案。作为输出，我想要一个数组，每个正确答案为 1，每个错误答案为 0。 for 循环可以工作，但是当您扩展到数千行和列时会非常慢。

这是我的可重现示例，提前感谢您的帮助！

   #build sample data
    dat <- array(dim=c(9,6))
    for (n in 1:9){
       dat[n,1:6] <- c(paste("ID00",n,sep=""),
           sample(c("A","B","C","D"), size=5, replace=TRUE))}
    dat[3,4]<-NA
    key<-c("key","A","B","B","C","D")
    dat <- rbind(key,dat)

>dat
[,1]    [,2] [,3] [,4] [,5] [,6]
"key"   "A"  "B"  "B"  "C"  "D" 
"ID001" "B"  "A"  "D"  "B"  "C" 
"ID002" "C"  "C"  "C"  "B"  "B" 
"ID003" "A"  "C"  NA   "D"  "D" 
"ID004" "D"  "B"  "D"  "A"  "A" 
"ID005" "A"  "C"  "A"  "C"  "A" 
"ID006" "D"  "D"  "B"  "B"  "A" 
"ID007" "B"  "D"  "A"  "D"  "A" 
"ID008" "D"  "D"  "B"  "D"  "A" 
"ID009" "D"  "C"  "B"  "D"  "D" 

    #score file
    dat2 <- array(dim=c(9,5))
    for (row in 2:10){
      for (column in 2:6){
        if (is.na(dat[row,column])){
          p <- NA
        }else if (dat[row,column]==dat[1,column]){
          p <- 1
        }else p <- 0
        dat2[row-1,column-1]<-p
      }
    }
> dat2
      [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    0    0    0    0    0
[3,]    1    0   NA    0    1
[4,]    0    1    0    0    0
[5,]    1    0    0    1    0
[6,]    0    0    1    0    0
[7,]    0    0    0    0    0
[8,]    0    0    1    0    0
[9,]    0    0    1    0    1

【问题讨论】：

我坚信解决这个问题的方法是彻底重新思考你的数据结构......我会试着编造一个例子。
我跑题了，你得到了其他答案……没关系。

标签： r for-loop apply

【解决方案1】：

为可重复性设置种子：

set.seed(1)
dat <- array(dim=c(9,6))
for (n in 1:9){
   dat[n,1:6] <- c(paste("ID00",n,sep=""),
       sample(c("A","B","C","D"), size=5, replace=TRUE))}
dat[3,4]<-NA
key<-c("key","A","B","B","C","D")
dat <- rbind(key,dat)

这样就可以了：

key <- rep(dat[1, -1], each = nrow(dat) - 1L)  ## expand "key" row
dummy <- (dat[-1, -1] == key) + 0L  ## vectorized / element-wise "=="

基本上我们想要一个矢量化的"=="。但我们首先需要将dat[1,-1] 扩展为与dat[-1,-1] 相同的维度。最后将+ 0L 强制TRUE / FALSE 矩阵转换为1 / 0 矩阵。

#  [,1] [,2] [,3] [,4] [,5]
#    0    1    0    0    0
#    0    0    0    1    0
#    1    0   NA    0    1
#    0    0    0    0    1
#    0    0    0    0    0
#    0    0    1    0    0
#    0    0    1    0    1
#    0    0    0    1    0
#    0    0    0    1    0

我还没有检查 Gregor 的基准脚本。但这是我的。

set.seed(1)
dat <- matrix(sample(LETTERS[4], 1000 * 1000, TRUE), 1000)
key <- sample(LETTERS[1:4], 1000, TRUE)
microbenchmark(rep(key, each = 1000) == dat, t(t(dat) == key))

#Unit: milliseconds
#                         expr      min       lq     mean   median       uq
# rep(key, each = 1000) == dat 32.16888 34.01138 42.61639 35.57526 40.27944
#             t(t(dat) == key) 50.93348 52.96008 63.74475 56.04706 60.38750
#       max neval cld
#  81.96044   100  a 
# 106.54916   100   b

我的方法和 Gregor 的唯一区别是 rep(, each) 扩展 vs. rep_len 扩展。两种扩展都消耗相同数量的内存，并且在扩展之后，"==" 以列方式完成。我预测额外的开销将由两个t() 引起，基准测试结果似乎证明了这一点。希望结果不依赖于平台。

【讨论】：

【解决方案2】：

这和哲元的回答基本一样（靠向量化的==再强制回数值），我只是先转置矩阵而不是展开键。

由于矩阵是按列而不是按行存储/操作的，因此如果键是一列并且每个学生也是一列，那么向量回收就可以了。

在生成数据之前使用set.seed(1)...

key = dat[1, -1]
tdat = t(dat[-1, -1])
t((tdat == key) + 0L)
 # [,1] [,2] [,3] [,4] [,5]
 #    0    1    0    0    0
 #    0    0    0    1    0
 #    1    0   NA    0    1
 #    0    0    0    0    1
 #    0    0    0    0    0
 #    0    0    1    0    0
 #    0    0    1    0    1
 #    0    0    0    1    0
 #    0    0    0    1    0

如果您改为将第一列更改为行名称，则可以轻松保留它们，而不会有将学生 ID 标记为不正确的风险，因为它们不是 'key'。这也使得最后总结事情变得更好：

row.names(dat) = dat[, 1]
dat = dat[, -1]
key = dat[1, ]   

tdat = t(dat[-1, ])
result = t((tdat == key) + 0)
result
#       [,1] [,2] [,3] [,4] [,5]
# ID001    0    1    0    0    0
# ID002    0    0    0    1    0
# ID003    1    0   NA    0    1
# ID004    0    0    0    0    1
# ID005    0    0    0    0    0
# ID006    0    0    1    0    0
# ID007    0    0    1    0    1
# ID008    0    0    0    1    0
# ID009    0    0    0    1    0

rowSums(result)
# ID001 ID002 ID003 ID004 ID005 ID006 ID007 ID008 ID009 
#     1     1    NA     1     0     1     2     1     1

简化输入并对中等规模的数据运行基准测试，两者都非常快。双转置要快一些。

gregor = function(key, dat) {
    t(t(dat) == key)
}

zheyuan = function(key, dat) {
    dat == rep(key, each = nrow(dat))
}

library(microbenchmark)
nr = 10000
nc = 1000
key = sample(1:10, nc, replace = T)
dat = matrix(sample(1:10, nr * nc, replace = T), nrow = nr)

print(microbenchmark(gregor(key, dat), zheyuan(key, dat)), signif = 4)
# Unit: milliseconds
#               expr   min    lq     mean median    uq   max neval cld
#   gregor(key, dat) 104.5 113.2 135.5970  128.2 144.5 336.2   100  a 
#  zheyuan(key, dat) 196.0 202.8 215.7822  207.0 224.9 394.4   100   b

identical(gregor(key, dat), zheyan(key, dat))
# [1] TRUE

【讨论】：

是的，我意识到我们都在不同的地方使用回收。
谢谢大家，我知道我遗漏了一些明显的东西！！
我很好奇它在更大的数据上会是什么 - 当然转置会引入开销，但创建与原始维度相同的 key 矩阵也是如此。
嗯，很高兴知道。我很惊讶，因为在简化的基准测试中，双重转置似乎更快

【解决方案3】：

如果您想在一行中不使用for 或apply，请尝试类似

dat2 <- matrix(as.numeric(dat==rep(dat[1,],each=nrow(dat))),nrow=nrow(dat))[-1,-1]

【讨论】：