循环的顺序对R中的速度是否重要

Question

我有一个问题，我想运行一个模拟研究，其中模拟取决于两个变量 x 和 y。 x 和 y 是我想在我的模拟研究中评估的潜在值的向量（所以不同的组合）。此外，对于x 和y 的每种组合，我想要多次重复（因为其中有一个随机项，x 和y 的每次运行都会有所不同）。

举一个我正在处理的例子，我有以下简化的例子：

x = 1:10
y = 11:20

iterations = 2000
iter = 1

solution = array(NA,c(length(x),3,iterations))
for(i in x){
    for(j in y){
        for(k in 1:iterations){ 
            z = rnorm(3) + c(i,j,1) 
            solution[i,,k] = z
        }
    }
}

然而，在我的实际问题中，在 for 循环中评估的代码的评估要简单得多。但是，我的输入和输出的结构是一样的。

所以我想知道的是，使用上面的示例，按该顺序设置循环是最有效的还是让k in 1:iterations 成为最外层循环并尝试使用某种outer() 在该 1 循环中的命令，因为我将在网格 x 和 y 上评估函数（在此示例中为 z）？

另外，我对完全不同的设置和设计持开放态度。归根结底，我希望能够获得基于x 和y 并在所有迭代中取平均值的解决方案，即apply(solution, c(1,2),mean)

---

编辑：

正如向我建议的那样，这是我正在使用的实际代码。

library(survival)

iter = 2000
n = 120
base = 2
hr = 0.5
m.x = 3
m.y = m.x/hr

ANS = NULL
for (vecX in c(0.3, 0.5, 0.6, 0.7)){
out = NULL

    for (vecY in c(0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95)){
        m.x.p = m.x/vecX
        m.y.p = m.y/vecX
        m.x.n = m.x
        m.y.n = m.y
        
        n.t = round(n*base/(base+1))
        n.c = n - n.t
        
        for (ii in 1:iter){
            n.t.p = rbinom(1, n.t, vecY)
            n.t.n = n.t - n.t.p
            n.c.p = rbinom(1, n.c, vecY)
            n.c.n = n.c - n.c.p
            
            S = c(rexp(n.t.p, log(2)/m.y.p), rexp(n.t.n, log(2)/m.y.n), rexp(n.c.p, log(2)/m.x.p), rexp(n.c.n, log(2)/m.x.n))
            
            data1 = data.frame(Group = c(rep("T", n.t), rep("C", n.c)), dx = c(rep("P", n.t.p), rep("N", n.t.n), rep("P", n.c.p), rep("N", n.c.n)), S)
            
            fit = survfit(Surv(data1$S)~data1$Group)
            coxfit = coxph(Surv(data1$S)~data1$Group)
            
            HR = exp(coxfit$coefficients)
            p.val=summary(coxfit)$logtest["pvalue"]
            
            out = rbind(out, c(vecX, vecY, n.t.p, n.t.n, n.c.p, n.c.n, HR, p.val))
        }
        
    }
    colnames(out) = c("vecX", "vecY", "n.t.p", "n.t.n", "n.c.p", "n.c.n", "HR", "p.val")
    
    ans = as.data.frame(out)
ANS = rbind(ANS, ans)

}

Answer 1

是的，我相信理论上它应该会有所作为（参见下面的示例）。

R 使用类似于 Fortran 的列优先排序（与 C 不同），因此为了最大限度地减少缓存未命中，您需要遍历列。因此，对于填充矩阵，最佳方法是外循环具有我们的列索引的方法。

对于 n 维数组，您也需要牢记这一点。在n = 3 的情况下，我想这意味着层是最外层的循环，然后是列，然后是行。不过，我可能在这里弄错了。

我使用5000 by 5000 矩阵运行了这个快速示例。我们看到大约 50 秒的差异，fill_matrix2() 更快。

    n <- 5000
    A <- matrix(NA, n, n)
    B <- matrix(NA, n, n)

    fill_matrix1 <- function(X, val) {
        for (i in 1:nrow(X)) {
            for (j in 1:ncol(X)) {
                X[i, j] <- val
            }
        }
        return(X)
    }

    fill_matrix2 <- function(X, val) {
        for (j in 1:ncol(X)) {
            for (i in 1:nrow(X)) {
                X[i, j] <- val
            }
        }
        return(X)
    }

    system.time(fill_matrix1(A, 0))
    system.time(fill_matrix2(B, 0))

Answer 2

循环的顺序在这里实际上是无关紧要的。如果您分析您的代码（请参阅help("Rprof")），您会发现 CPU 时间花费在像 survfit 和 coxph 这样的函数上。当然在增长out，你应该避免。将out 预先分配到其最终大小并填充它而不是增大它。