Julia Optimization

user_14415173 0
Viewed 129

I am trying to run a bootstrap in Julia and have coded up a working bootstrap function. However, it is slow and the same code in R runs in half the time. I am sure that there must be some inefficiencies in my code, I am very new to using Julia. I was wondering if anyone could provide me with some advice/suggestions.

Here is the fully reproducible code

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);

This takes around 68 seconds to run it takes 35 seconds in R.

Your advice is greatly appreciated. Thank you.

julia dataframe
1 Answers

Using DataFrames.jl you can do e.g.:

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)]

Which should be faster than R. The changes are:

  • major: use views instead of copying everything in every iteration
  • minor: do not create a data frame for each bootstrap replicate but rather a vector and push! it instead of append!ing it (this saves time of creation and validation of data frame objects)

(I have made only the major optimizations of the code; there are some additional minor optimizations that could be made, but they should not affect the run time in a significant way)

Also note that you are close to the maximum execution speed as:

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);

Is roughly the lower bound of execution time you can expect to have and it is not much faster than the code above (it is faster of course, by around 25%-30%, as it does less work and is more CPU cache friendly).

As a small comment showing how details matter in such cases (I think it is interesting, although it is a minor optimization, so I left it out).

Instead of rand(1:nrow(data), nrow(data)) if you use sort!(rand(1:nrow(data), nrow(data))) you save an additional 1 second. The reason is that in this way you ensure that you access data sequentially when you calculate mean (which is more CPU cache friendly and mean is unaffected by observation order).

A second comment like this is that on a multi-CPU machine (and started Julia with -t switch selecting to use more than one thread) one could use threading to speed things up like this (again - I did not optimize things out here to the very last possible tweak, but rather wanted to show the main idea):

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

This is much faster and easy to do in Julia (while doable, but not so easy in R).

Regarding views you can read about them here.

user_12695670
Related