Julia 循环和 R 循环一样慢

Question

下面的 Julia 和 R 中的代码是为了表明总体方差的估计量是有偏的估计量，即它取决于样本量，并且无论我们对不同的观察平均多少次，对于少量数据点它不等于总体的方差。

Julia 大约需要 10 秒才能完成两个循环，而 R 大约需要 7 秒。 如果我将代码留在注释的循环中，那么 R 和 Julia 中的循环将花费相同的时间，并且如果我只将迭代器相加 s = s + i+ j Julia 将在 ~0.15 秒内完成，R 在 ~0.5 秒内完成。

是 Julia 循环变慢还是 R 变快了？ 如何为 Julia 提高以下代码的速度？ R代码可以变快吗？

朱莉娅：

using Plots
trials = 100000
sample_size = 10;
sd = Array{Float64}(trials,sample_size-1)

tic()
for i = 2:sample_size
    for j = 1:trials
        res = randn(i)
        sd[j,i-1] = (1/(i))*(sum(res.^2))-(1/((i)*i))*(sum(res)*sum(res))
    end
end
toc()
sd2 = mean(sd,1)
plot(sd2[1:end])

R：

trials = 100000
sample_size = 10
sd = matrix(, nrow = trials, ncol = sample_size-1)
start_time = Sys.time()
for(i in 2:sample_size){
  for(j in 1:trials){
  res <- rnorm(n = i, mean = 0, sd = 1)
  sd[j,i-1] = (1/(i))*(sum(res*res))-(1/((i)*i))*(sum(res)*sum(res))

}
}

end_time = Sys.time()
end_time - start_time
sd2 = apply(sd,2,mean)
plot(sqrt(sd2))

情节以防有人好奇！：

我可以实现更高速度的一种方法是使用在 Julia 中很容易实现的并行循环：

using Plots
trials = 100000
sample_size = 10;


sd = SharedArray{Float64}(trials,sample_size-1)

tic()
@parallel for i = 2:sample_size
    for j = 1:trials
        res = randn(i)
        sd[j,i-1] = (1/(i))*(sum(res.^2))-(1/((i)*i))*(sum(res)*sum(res))
    end
end
toc()
sd2 = mean(sd,1)
plot(sd2[1:end])

Answer 1

在 Julia 中使用全局变量通常速度很慢，并且速度应该可以与 R 相媲美。您应该将代码包装在一个函数中以使其更快。

这是我笔记本电脑上的时间（我只剪掉了相关部分）：

julia> function test()
           trials = 100000
           sample_size = 10;
           sd = Array{Float64}(trials,sample_size-1)

           tic()
           for i = 2:sample_size
               for j = 1:trials
                   res = randn(i)
                   sd[j,i-1] = (1/(i))*(sum(res.^2))-(1/((i)*i))*(sum(res)*sum(res))
               end
           end
           toc()
       end
test (generic function with 1 method)

julia> test()
elapsed time: 0.243233887 seconds
0.243233887

另外，在 Julia 中，如果您使用 randn! 而不是 randn，则可以加快速度，因为您避免重新分配 res 向量（我没有对代码进行其他优化，因为这种优化与 Julia 不同与 R 相比；此代码中所有其他可能的加速都将以类似的方式帮助 Julia 和 R）：

julia> function test2()
           trials = 100000
           sample_size = 10;
           sd = Array{Float64}(trials,sample_size-1)

           tic()
           for i = 2:sample_size
               res = zeros(i)
               for j = 1:trials
                   randn!(res)
                   sd[j,i-1] = (1/(i))*(sum(res.^2))-(1/((i)*i))*(sum(res)*sum(res))
               end
           end
           toc()
       end
test2 (generic function with 1 method)

julia> test2()
elapsed time: 0.154881137 seconds
0.154881137

最后最好使用BenchmarkTools 包来测量 Julia 中的执行时间。第一个 tic 和 toc 函数将从 Julia 0.7 中删除。其次 - 如果你使用它们，你会混合编译和执行时间（当运行 test 函数两次时，你会发现第二次运行的时间减少了，因为 Julia 不花时间编译函数）。

编辑：

您可以将trials、sample_size 和sd 保留为全局变量，但您应该在它们前面加上const。然后在这样的函数中包装一个循环就足够了：

const trials = 100000;
const sample_size = 10;
const sd = Array{Float64}(trials,sample_size-1);

function f()
    for i = 2:sample_size
        for j = 1:trials
            res = randn(i)
            sd[j,i-1] = (1/(i))*(sum(res.^2))-(1/((i)*i))*(sum(res)*sum(res))
        end
    end
end

tic()
f()
toc()

现在@parallel：

首先，您应该在@parallel 之前使用@sync 以确保所有工作正常（即，在您转到下一条指令之前所有工作人员都已完成）。要了解为什么需要这样做，请在具有多个工作人员的系统上运行以下代码：

sd = SharedArray{Float64}(10^6);
@parallel for i = 1:2
    if i < 2
        sd[i] = 1
    else
        for j in 2:10^6
            sd[j] = 1
        end
    end
end
minimum(sd) # most probably prints 0.0
sleep(1)
minimum(sd) # most probably prints 1.0

此时

sd = SharedArray{Float64}(10^6);
@sync @parallel for i = 1:2
    if i < 2
        sd[i] = 1
    else
        for j in 2:10^6
            sd[j] = 1
        end
    end
end
minimum(sd) # always prints 1.0

其次，速度的提升是由于@parallel 宏而不是SharedArray。如果您在 Julia 上与一名工作人员一起尝试您的代码，它也会更快。简而言之，原因是 @parallel 在内部将您的代码包装在一个函数中。可以使用@macroexpand查看：

julia> @macroexpand @sync @parallel for i = 2:sample_size
           for j = 1:trials
               res = randn(i)
               sd[j,i-1] = (1/(i))*(sum(res.^2))-(1/((i)*i))*(sum(res)*sum(res))
           end
       end
quote  # task.jl, line 301:
    (Base.sync_begin)() # task.jl, line 302:
    #19#v = (Base.Distributed.pfor)(begin  # distributed\macros.jl, line 172:
                function (#20#R, #21#lo::Base.Distributed.Int, #22#hi::Base.Distributed.Int) # distributed\macros.jl, line 173:
                    for i = #20#R[#21#lo:#22#hi] # distributed\macros.jl, line 174:
                        begin  # REPL[22], line 2:
                            for j = 1:trials # REPL[22], line 3:
                                res = randn(i) # REPL[22], line 4:
                                sd[j, i - 1] = (1 / i) * sum(res .^ 2) - (1 / (i * i)) * (sum(res) * sum(res))
                            end
                        end
                    end
                end
            end, 2:sample_size) # task.jl, line 303:
    (Base.sync_end)() # task.jl, line 304:
    #19#v
end