Julia优化

Question

我试图在Julia中运行一个引导程序，并编写了一个可工作的引导函数。然而，这是缓慢和相同的代码在R运行了一半的时间。我确信我的代码中一定有一些效率低下的地方，我对使用Julia非常陌生。我想知道是否有人能给我一些建议。

以下是完全可复制的代码

using DataFrames
using Statistics
using StatsBase

df = DataFrame(rand(1:9, 1000,1000), :auto); # Create data

# Bootstrap function
function bootstrap(;iters=1, data=nothing, statistic=nothing)
    statArr = DataFrame() # Init empty dataframe
    for i in 1:iters
        data_sample = data[sample(1:nrow(data), nrow(data), replace=true), :] # sample the data with replacement
        stat = statistic(data_sample)
        append!(statArr, stat) # push dataframe to empty dataframe
    end
    return statArr
end;

# Statistic function for column means
function meanmap(data)
    return mapcols(col -> mean(col), data)
end;

# Run the bootstrap on the data
@time bootDist = bootstrap(iters = 9999, data = df, statistic = meanmap);

这大约需要68秒来运行，在R中需要35秒。

我们非常感谢你的建议。谢谢。

Answer 1

使用DataFrames.jl，您可以这样做，例如：

function bootstrap(;iters=1, data=nothing, statistic=nothing)
    statArr = Float64.(empty(data)) # Init empty dataframe
    for i in 1:iters
        stat = statistic(data, rand(1:nrow(data), nrow(data)))
        push!(statArr, stat) # push row to empty dataframe
    end
    return statArr
end;

# Statistic function for column means
meanmap(data, sel) = [mean(@view x[sel]) for x in eachcol(data)]

应该比R更快的变化是：

• 专业:使用视图而不是在每次迭代中复制所有内容。
• 次要:不要为每个引导复制创建一个数据框架，而是创建一个向量和push!，而不是append!(这样可以节省创建和验证数据帧对象的时间)

(我只对代码进行了主要的优化；还可以进行一些额外的小优化，但它们不应该对运行时产生重大影响)

还请注意，您接近最大执行速度，如下所示：

julia> x = rand(1:nrow(df), nrow(df));

julia> y = df[!, 1];

julia> f(y, x) = mean(@view y[x]);

julia> g(y, x) = [f(y, x) for _ in 1:9999*1000];

julia> @time g(y, x);

大约是您可以预期的执行时间的下限，而且它并不比上面的代码快多少(当然，它比上面的代码快25%-30%，因为它的工作量更小，并且对CPU缓存更友好)。

作为一个小评论，展示了在这种情况下细节是如何重要的(我认为这很有趣，尽管这是一个小的优化，所以我忽略了它)。

与rand(1:nrow(data), nrow(data))不同，如果使用sort!(rand(1:nrow(data), nrow(data)))，则可以节省1秒。原因是这样，您确保在计算mean时按顺序访问数据(这对CPU缓存更加友好，而且mean不受观察顺序的影响)。

第二个评论是，在一台多CPU计算机上(并在-t开关中选择使用多个线程启动了Julia )，您可以使用线程来加快这样的速度(同样，我并没有在这里优化到最后一个可能的调整，而是想展示主要的想法)：

function bootstrap(;iters=1, data=nothing, statistic=nothing)
    statArr = Float64.(empty(data)) # Init empty dataframe
    tmp = Vector{Any}(undef, iters)
    Threads.@threads for i in 1:iters
        stat = statistic(data, rand(1:nrow(data), nrow(data)))
        tmp[i] = stat
    end
    for v in tmp
        push!(statArr, v) # push dataframe to empty dataframe
    end
    return statArr
end

在Julia (虽然可行，但在R中不那么容易)，这是更快、更容易做到的。

关于视图，您可以阅读有关它们的这里。