Up until now, every line of code you've written in R has been running on a single CPU core. For most scripts, that's plenty. But once you start dealing with massive simulations or iterating over thousands of heavy computations, you'll notice R just sitting there, chewing through one task at a time while the other seven or fifteen cores on your machine are basically idling. It feels like a waste.

The parallel package is the standard way to wake those other cores up. Today, we're going to use it to estimate the value of Pi using a Monte Carlo simulation. It's a computationally "expensive" way to find Pi—essentially throwing random darts at a square and seeing how many land in a circle—which makes it the perfect candidate for parallelization.

The bottleneck: Estimating Pi with a shotgun approach

First, let's write a simple function that performs one "batch" of these random throws. We'll generate random X and Y coordinates between 0 and 1, and check if they fall within the unit circle.

estimate_pi <- function(n) {
  x <- runif(n)
  y <- runif(n)
  # A point is inside the circle if x^2 + y^2 <= 1
  inside <- sum(x^2 + y^2 <= 1)
  return(inside / n * 4)
}

If I want to run this simulation 100 times with 1 million points each to get a stable average, doing it sequentially with lapply is slow. I can feel the lag. Instead, I want to split those 100 iterations across all my available cores.

Spinning up the worker cluster

To do this, we need to create a "cluster." Think of this as launching several invisible, separate instances of R in the background that just wait for instructions from your main session.

library(parallel)

# Detect how many cores I have
num_cores <- detectCores()
# I'll leave one core free so my computer doesn't freeze up completely
cl <- makeCluster(num_cores - 1)

Now that the cluster is live, I'll use parLapply. It's almost identical to lapply, but the "par" stands for parallel. It tells R to ship the work off to those background instances we just created.

Wait, where did my function go?

Here is where I almost always trip up the first time I use this package. I'll try to run the simulation like this:

# This is where I usually make my mistake
iterations <- 1:100
results <- parLapply(cl, iterations, function(i) estimate_pi(1e6))

If you ran that, you'd get an error saying could not find function "estimate_pi". It's frustrating because the function is right there in your environment! But remember: those worker nodes are separate R processes. They don't share your global environment. They are blank slates. They have no idea what estimate_pi is because I never told them about it.

Shipping the logic to the workers

To fix this, we have to explicitly "export" the function (and any other variables the function depends on) to every node in the cluster using clusterExport. This is the missing link.

# Export the function to the worker nodes
clusterExport(cl, "estimate_pi")

# Now we run it again
results <- parLapply(cl, iterations, function(i) estimate_pi(1e6))

# Calculate the final average of our simulations
final_pi <- mean(unlist(results))
print(final_pi)

The speed difference is immediate. Instead of one core doing 100% of the work, the load is distributed. If you have 8 cores, you're theoretically cutting your wait time down significantly.

Tearing down the cluster
One thing you cannot forget: those worker nodes don't just vanish when your script ends. They stay alive in your system memory. If you keep creating clusters without closing them, you'll eventually run out of RAM and your OS will start screaming.

stopCluster(cl)

Always pair your makeCluster with a stopCluster. I usually put the stop command in a finally block or immediately after the parallel call to make sure I don't leave "ghost" R sessions haunting my machine.

📋 Practical Task

Parallelizing Dice-Roll Simulations

You need to determine the theoretical average of a 6-sided die roll (which is 3.5) using a simulation. Instead of one giant simulation, you will run 500 separate simulations, each rolling the die 10,000 times.

Write a script that does the following:

1. Defines a function roll_dice that takes an integer n, simulates n rolls of a 6-sided die using sample(), and returns the mean() of those rolls.
2. Creates a cluster using makeCluster() based on the number of available cores.
3. Uses clusterExport() to ensure the worker nodes can see the roll_dice function.
4. Uses parLapply() to run the simulation 500 times (with 10,000 rolls per simulation).
5. Calculates the overall mean of all the simulation results.
6. Shuts down the cluster using stopCluster().