Multithreading in Julia

Fortunately, since Julia v1.3, multithreading has been fairly simple to use inside of Julia. Performance of multithreading has also increased as Julia develops. For this reason, results shown in this book may vary greatly, depending on when it is read, and the current version of Julia you are using.

This section will cover a basic use of multithreading to speed up an embarrassingly parallel algorithm. An embarrassingly parallel problem is one where we have a one-to-one relationship between our inputs (parameters) and our desired outputs (results). This sort of problem can be distilled down to repeating the following operation $n$n times:

$r_i = f(x_i)$ri​=f(xi​)

where $x_i$xi​ represents one element of the $n$n elements that are to be processed, $f$f is the common function that performs the mapping and $r_i$ri​ represents the $i^{\text{th}}$ith element of the results vector. To use Julia specific language, any operation that can be written in vector format (broadcasted) is likely a good candidate to be embarrassingly parallel.

Monte Carlo Estimation of π

For a case study, let’s look at a way that we can calculate $\pi$π. We will choose a Monte-Carlo method of estimation, since this will prove to be an excellent case study for multithreading. Monte-Carlo methods are sampling techniques that can be used to estimate quantities via random simulations. They are often very easy to implement and can be used as a simple technique to provide numerical estimates for quantities that difficult or impossible to calculate analytically. In the case of calculating $\pi$π, one can estimate its value by playing darts on a square block, with a circle touching the edge of the square. If the radius of the circle is $r$r, then the area is $\pi r^2$πr2 and the area of the square it lies within is $4r^2$4r2. We know that the ratio of the area of the circle compared with the square is $\rho = \frac{\pi r^2}{4 r^2}=\frac{\pi}{4}$ρ=4r2πr2​=4π​. If we can calculate an estimate for $\rho$ρ, we can approximate $\pi$π as $\pi = 4 \rho$π=4ρ.

We can now randomly play darts on this setup, making sure that each point in the square is equally likely to be somewhere a dart hits. Then we can count how many of the darts land in the circle, compared with landing outside the circle and in the square. If we throw $n$n darts uniformly randomly in the square and only $n_c$nc​ land inside the circle, then we know that $\rho \approx \frac{n_c}{n}$ρ≈nnc​​, which will approach the true value of $\rho$ρ when $n\rightarrow \infty$n→∞, or in other terms:

$\pi = 4\rho = 4 \lim_{n\rightarrow \infty} {\frac{n_c}{n}}$π=4ρ=4limn→∞​nnc​​

First, let’s write an algorithm to do this in serial.

function est_pi_mc_serial(n)
    n_c = zero(typeof(n))
    for _ in 1:n
        # Choose random numbers between -1 and +1 for x and y
        x = rand() * 2 - 1
        y = rand() * 2 - 1
        # Work out the distance from origin using Pythagoras
        r2 = x*x+y*y
        # Count point if it is inside the circle (r^2=1)
        if r2 <= 1
            n_c += 1
        end
    end
    return 4 * n_c / n
end

One can see that the estimate for $\pi$π improves with the number of darts thrown, by using a much higher number of points:

function est_pi_mc_serial(n)
    n_c = zero(typeof(n))
    for _ in 1:n
        # Choose random numbers between -1 and +1 for x and y
        x = rand() * 2 - 1
        y = rand() * 2 - 1
        # Work out the distance from origin using Pythagoras
        r2 = x*x+y*y
        # Count point if it is inside the circle (r^2=1)
        if r2 <= 1
            n_c += 1
        end
    end
    return 4 * n_c / n
end

pi
est_pi_mc_serial(100)
est_pi_mc_serial(10_000)

Additionally, one can plot the standard deviation of the relative errors of each estimate of $\pi$π against the number of darts thrown. This can be seen in Figure 1, which shows that the relative error approaches zero as $n\rightarrow \infty$n→∞.

While this is not the most efficient way of calculating $\pi$π, heavily relying on the quality of the random number generator, it does provide a case study for a parallel speed up.