Case Study: Monte-Carlo Simulations
Here, we will provide an example of porting a Monte-Carlo random walk simulation to the GPU, such as the one shown in Figure 1. As with the earlier advice in this course, we will write a function that performs a single step in the random walk.
function mc_random_walk_step(x, sigma)
return x + randn(typeof(x)) * sigma
endNow imagine that we want to study some population level statistics, by running many of these walkers in parallel. Despite us choosing a very simple example, Monte-Carlo simulations are an extremely useful tool for computational science. Our aim for this exercise is to perform some number of steps of this Monte-Carlo update for many independent walkers and obtain an array with their final positions in.
We can implement a non-allocating array version of our desired algorithm:
function mc_random_walk!(y, x, sigma, steps)
# Copy the initial values from x into y
y .= x
for t in 1:steps
y .= mc_random_walk_step.(y, sigma)
end
return nothing
endWe can run our algorithm on the CPU easily:
n=2048; x = zeros(Float32, n); y = similar(x);
sigma = 1.0f0; steps=100;
mc_random_walk!(y, x, sigma, steps);We can extend this to run on the GPU just by changing the types:
x_gpu = to_gpu(x); y_gpu = similar(x_gpu);
mc_random_walk!(y_gpu, x_gpu, sigma, steps);Let’s benchmark the CPU version:
display(@benchmark mc_random_walk!($y, $x, $sigma, $steps))BenchmarkTools.Trial: 4828 samples with 1 evaluation per sample.
Range (min … max): 898.110 μs … 2.498 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 963.186 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.035 ms ± 145.548 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▇██▆▅▆▆▅▃▂▁▁▁ ▁ ▁ ▁▄▅▅▅▃▂▁ ▂
▆██████████████████████▇▇▇██▇▇█▇█▇▇▇▇▇▇▇▇▇███████████▆▅█▇▄▆▃▄ █
898 μs Histogram: log(frequency) by time 1.39 ms <
Memory estimate: 0 bytes, allocs estimate: 0.And the GPU version:
display(@benchmark begin
mc_random_walk!($y_gpu, $x_gpu, $sigma, $steps)
gpu_synchronize($y_gpu)
end)BenchmarkTools.Trial: 3143 samples with 1 evaluation per sample.
Range (min … max): 1.336 ms … 7.695 ms ┊ GC (min … max): 0.00% … 65.61%
Time (median): 1.496 ms ┊ GC (median): 0.00%
Time (mean ± σ): 1.590 ms ± 417.815 μs ┊ GC (mean ± σ): 1.42% ± 4.72%
▁ ▁█▆ ▁▃
▃▇█████▆▅▄▄▄▄▇███▅▄▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▁▂▁▂▂▂▂▁▂▂▂ ▃
1.34 ms Histogram: frequency by time 2.62 ms <
Memory estimate: 295.11 KiB, allocs estimate: 8081.We see that our GPU version is actually faster (for this number of walkers), however, there is actually a bit of a performance bug under hood that should be addressed. Upon each of our for loops, we are calling a new kernel. Instead, we should try to fuse these kernels together:
function mc_random_walk(x, sigma, steps)
for t in 1:steps
x = mc_random_walk_step(x, sigma)
end
return x
end
function mc_random_walk_fused!(y, x, sigma, steps)
y .= mc_random_walk.(x, sigma, steps)
return nothing
endNow we can try and benchmark this again:
display(@benchmark begin
mc_random_walk_fused!($y_gpu, $x_gpu, $sigma, $steps)
gpu_synchronize($y_gpu)
end)BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 230.100 μs … 3.355 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 377.750 μs ┊ GC (median): 0.00%
Time (mean ± σ): 353.879 μs ± 66.223 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▁ ▄▇█▅▃▁
▁▁▂▂▂▃▆██▆▄▄▅▅▄▄▄▄▃▃▃▃▃▃▂▂▂▃▃▄▄▄▅████████▆▅▄▂▂▁▂▁▁▁▁▁▁▁▁▁▁▁▁ ▃
230 μs Histogram: frequency by time 490 μs <
Memory estimate: 3.02 KiB, allocs estimate: 81.We have fused all the kernel calls together. This minimises the amount of overhead for scheduling, and allows the GPU cores to stay busy for the duration of the computation. This small change allowed us to dramatically improve the performance of our GPU code.
If the kernel calls are very large (e.g. a large matrix multiply), the overhead in calling multiple kernels is only a very small portion of the total time.
Custom Kernel
It is a good exercise to write a custom kernel and compare the execution to the broadcasted notation. Our kernel will be straightforward:
@kernel function mc_random_walk_kernel!(y, x, increments, sigma, steps)
i = @index(Global, Linear)
if i <= length(y)
@inbounds pos = x[i]
for t in 1:steps
@inbounds pos += increments[i, t] * sigma
end
@inbounds y[i] = pos
end
end
function mc_random_walk_gpu!(y, x, sigma, steps)
@assert length(y) == length(x)
increments = to_gpu(randn(eltype(x), length(x), steps))
backend = get_backend(y)
mc_random_walk_kernel!(backend, 256)(y, x, increments, sigma, steps; ndrange=length(y))
return nothing
endWe can now benchmark on the same data:
display(@benchmark begin
mc_random_walk_gpu!($y_gpu, $x_gpu, $sigma, $steps)
gpu_synchronize($y_gpu)
end)BenchmarkTools.Trial: 4253 samples with 1 evaluation per sample.
Range (min … max): 948.450 μs … 3.572 ms ┊ GC (min … max): 0.00% … 31.77%
Time (median): 1.113 ms ┊ GC (median): 0.00%
Time (mean ± σ): 1.175 ms ± 396.332 μs ┊ GC (mean ± σ): 4.09% ± 7.69%
█▇▆▅█▇▆▅▄▂▁▁ ▂
████████████▇▇▆▆▅▆▃▅▃▁▃▃▁▁▄▃▄▄▁▃▁▁▁▁▁▃▁▁▁▃▁▁▄▁▁▅▁▅▆▇▅▇▇▇▇█▇██ █
948 μs Histogram: log(frequency) by time 3.09 ms <
Memory estimate: 802.62 KiB, allocs estimate: 74.We can see that our custom kernel did not perform as well as our simpler approach. It is clear that writing a custom kernel is not necessary to achieve performance gains, as we can rely on the compiler to generate fast code for the GPU, using the array notation. It is entirely possible to rewrite our kernel to be of a similar performance to our previous implementation, but this would require more effort on our part.