我有一個 2x4 矩陣 A 和一個 3x4 矩陣 B。我想得到一個維度為 (2, 3, 4) 的陣列 C，其中C 的條目是 A 的條目和 B的條目ijk的乘積。ikjk
有沒有一種快速的方法可以通過避免回圈來在 R 中做到這一點？下面的示例有兩種方法來計算我正在尋找的東西——這兩種方法都涉及 for 回圈
A <- matrix(1:8, 2, 4)
B <- matrix(11:22, 3, 4)
C1 <- array(NA, dim=c(2, 3, 4))
for (ii in 1:2) {
for (jj in 1:3) {
C1[ii, jj, ] <- A[ii, ] * B[jj, ]
}
}
C2 <- array(NA, dim=c(2, 3, 4))
for (ss in 1:4) {
C2[, , ss] <- outer(A[, ss], B[, ss])
}
uj5u.com熱心網友回復：
沒有太大的改進，但使用abind：
C3 <- do.call(abind::abind, c(lapply(seq(ncol(A)), function(ss) outer(A[,ss], B[,ss])), along=3))
基準
注意：我上一次運行此基準測驗是在另一臺運行 R-4.0.5 的筆記本電腦上運行的，在這種情況下，多次運行基準測驗，其C1性能與C2's 相當。我換了筆記本電腦（出于其他原因），看到了@jay.sf 的答案并想將其添加到戰斗中，現在C1明顯更快了。我無法解釋不同之處，但相關規格：
- 之前的基準測驗：Windows 10，R-4.0.5
- 這個基準：windows 11, R-4.1.2
在我看來，基準測驗中的“動蕩”是由于資料量小，受到管理開銷的嚴重影響。如果矩陣更大，我會期待更好的突破（無需驗證）。
bench::mark(
C1 = {
for (ii in 1 : 2) {
for (jj in 1 : 3) {
C1[ii, jj, ] <- A[ii, ] * B[jj, ]
}
}
as.numeric(C1)
},
C2 = {
for (ss in 1 : 4) {
C2[, , ss] <- outer(A[, ss], B[, ss])
}
as.numeric(C2)
},
C3 = as.numeric(do.call(abind::abind, c(lapply(seq(ncol(A)), function(ss) outer(A[,ss], B[,ss])), along=3))),
C4 = as.numeric(vapply(1:4, \(ss) outer(A[, ss], B[, ss]), matrix(0, 2, 3)))
)
# # A tibble: 4 x 13
# expression min median `itr/sec` mem_alloc `gc/sec` n_itr n_gc total_time result memory time gc
# <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl> <int> <dbl> <bch:tm> <list> <list> <list> <list>
# 1 C1 16.6us 52.3us 18669. 240B 0 9317 0 499ms <dbl [24]> <Rprofmem [1 x 3]> <bench_tm [9,317]> <tibble [9,317 x 3]>
# 2 C2 25.8us 85.6us 11225. 240B 2.14 5235 1 466ms <dbl [24]> <Rprofmem [1 x 3]> <bench_tm [5,236]> <tibble [5,236 x 3]>
# 3 C3 763us 797.4us 1144. 720B 0 572 0 500ms <dbl [24]> <Rprofmem [3 x 3]> <bench_tm [572]> <tibble [572 x 3]>
# 4 C4 24.8us 36.7us 21752. 240B 2.18 9999 1 460ms <dbl [24]> <Rprofmem [1 x 3]> <bench_tm [10,000]> <tibble [10,000 x 3]>
（我添加as.numeric到它們中的每一個是因為一些回傳整數，一些雙精度數，我不想as.numeric一個并且偏向基準。這對它們中的任何一個都不是嚴格要求的，但現在我們可以確信它們都是平等的，因為否則bench::mark會失敗，抱怨輸出不同。）
我喜歡@jay.sf 的答案，因為它既快速又不需要額外的軟體包（abind雖然是非標準的，但很方便）。
基準，取 2
讓我們稍微增加一下資料。
A2 <- do.call(cbind, replicate(50, do.call(rbind, replicate(50, A, simplify = FALSE)), simplify = FALSE))
Abig <- do.call(cbind, replicate(50, do.call(rbind, replicate(50, A, simplify = FALSE)), simplify = FALSE))
Bbig <- do.call(cbind, replicate(50, do.call(rbind, replicate(50, B, simplify = FALSE)), simplify = FALSE))
C1big <- C2big <- array(dim=c(dim(Abig)[1], dim(Bbig)))
dim(Abig)
# [1] 100 200
現在是一個不同的基準，現在在 Windows 11，R-4.1.2 上：
bench::mark(
C1big = {
for (ii in seq(nrow(Abig))) {
for (jj in seq(nrow(Bbig))) {
C1big[ii, jj, ] <- Abig[ii, ] * Bbig[jj, ]
}
}
as.numeric(C1big)
},
C2big = {
for (ss in seq(ncol(Abig))) {
C2big[, , ss] <- outer(Abig[, ss], Bbig[, ss])
}
as.numeric(C2big)
},
C3big = as.numeric(do.call(abind::abind, c(lapply(seq(ncol(Abig)), function(ss) outer(Abig[,ss], Bbig[,ss])), along=3))),
C4big = as.numeric(vapply(seq(ncol(Abig)), function(ss) outer(Abig[, ss], Bbig[, ss]), matrix(0, nrow(Abig), nrow(Bbig)))),
iterations = 30
)
# # A tibble: 4 x 13
# expression min median `itr/sec` mem_alloc `gc/sec` n_itr n_gc total_time result memory time gc
# <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl> <int> <dbl> <bch:tm> <list> <list> <list> <list>
# 1 C1big 100.7ms 135ms 7.38 83.5MB 2.68 22 8 2.98s <dbl [3,000,000]> <Rprofmem [75,001 x 3]> <bench_tm [30]> <tibble [30 x 3]>
# 2 C2big 22ms 24.4ms 36.0 46.8MB 7.20 25 5 694.7ms <dbl [3,000,000]> <Rprofmem [1,801 x 3]> <bench_tm [30]> <tibble [30 x 3]>
# 3 C3big 26.5ms 28.4ms 32.8 92.6MB 28.7 16 14 488.54ms <dbl [3,000,000]> <Rprofmem [1,632 x 3]> <bench_tm [30]> <tibble [30 x 3]>
# 4 C4big 10ms 17.4ms 61.2 46.7MB 18.6 23 7 375.9ms <dbl [3,000,000]> <Rprofmem [1,402 x 3]> <bench_tm [30]> <tibble [30 x 3]>
似乎我們已經達到了我的預期：jay.sf 的vapply實作應該做得很好，而且它似乎正在突破`itr/sec`并mem_alloc成為很好的指標。
uj5u.com熱心網友回復：
vapply() 應該很快。
vapply(1:4, \(ss) outer(A[, ss], B[, ss]), matrix(0, 2, 3))
# , , 1
#
# [,1] [,2] [,3]
# [1,] 11 12 13
# [2,] 22 24 26
#
# , , 2
#
# [,1] [,2] [,3]
# [1,] 42 45 48
# [2,] 56 60 64
#
# , , 3
#
# [,1] [,2] [,3]
# [1,] 85 90 95
# [2,] 102 108 114
#
# , , 4
#
# [,1] [,2] [,3]
# [1,] 140 147 154
# [2,] 160 168 176
sapply 也是可能的，但速度較慢（輸出相同）：
sapply(1:4, \(ss) outer(A[, ss], B[, ss]), simplify='array')
基準：
AA <- matrix(1:8, 1e3, 4e2)
BB <- matrix(11:22, 2e3, 4e2)
microbenchmark::microbenchmark(
`for1`={C1 <- array(NA, dim=c(1e3, 2e3, 4e2))
for (ii in 1:2) {
for (jj in 1:3) {
C1[ii, jj, ] <- AA[ii, ] * BB[jj, ]
}
}},
vapply=vapply(1:4, \(ss) outer(AA[, ss], BB[, ss]), matrix(0, 1e3, 2e3)),
abind=do.call(abind::abind, c(lapply(seq(ncol(AA)), function(ss) outer(AA[,ss], BB[,ss])), along=3)),
sapply=sapply(1:4, \(ss) outer(AA[, ss], BB[, ss]), simplify='array'),
times=3L, control=list(warmup=1e2L))
# Unit: milliseconds
# expr min lq mean median uq max neval cld
# for1 3766.87037 3792.71469 3837.77384 3818.5590 3873.2256 3927.8921 3 b
# vapply 34.04493 68.17596 98.97098 102.3070 131.4340 160.5610 3 a
# abind 11736.37882 12063.20849 12320.85601 12390.0382 12613.0946 12836.1511 3 c
# sapply 58.41669 81.65372 139.44338 104.8907 179.9567 255.0227 3 a
轉載請註明出處，本文鏈接：https://www.uj5u.com/shujuku/417899.html
標籤：