Speeding Up Monte Carlo Pi Estimation with MATLAB parfor

Why a Serial Loop is a Bottleneck

Estimating π with a Monte Carlo method requires millions of random point trials. In a standard for loop each iteration runs one after another, even on a multi‑core machine. The CPU core that is idle while the loop is waiting for the next iteration is wasted, and the total run time grows linearly with the number of samples.

When parfor Actually Helps

A parfor loop only gives a speedup when the body of the loop is *independent* and *computationally heavy* enough to outweigh the overhead of distributing work to separate MATLAB workers. Typical indicators:

• The loop body performs floating‑point math or calls external libraries.
• Each iteration takes at least a few microseconds.
• The problem size (number of iterations) is in the order of 10⁶ or larger.

If the loop body is trivial or the dataset is tiny, the cost of launching a parallel pool and moving data between workers can dominate the total time.

Preparing the Parallel Environment

1. Check for the Toolbox

   Before you call parfor you should verify that the Parallel Computing Toolbox license is available. This prevents a runtime error on machines that lack the license.

   if ~license('test','Parallel_Computing_Toolbox')
       error('Parallel Computing Toolbox license required.');
   end

2. Initialize the Random Number Stream

   To guarantee statistically independent streams across workers you can set the global seed once, before the parfor. MATLAB will split the stream automatically.

   rng('shuffle'); % uses current time as seed
   

3. Start a Parallel Pool

   Start a pool only if one is not already active. You can let MATLAB choose the default number of workers, or explicitly request a specific count.

   if isempty(gcp('nocreate'))
       parpool(); % default pool
   end

Worked Example: π Estimation with parfor

The following script demonstrates a full, runnable example that compares serial and parallel execution times. It is written for MATLAB R2023b or newer.

% ----- Monte Carlo Pi Estimation -----
% Parameters
N = 1e7; % number of random points

% Ensure Parallel Computing Toolbox is present
if ~license('test','Parallel_Computing_Toolbox')
    error('Parallel Computing Toolbox license required.');
end

% Start a parallel pool if not already running
if isempty(gcp('nocreate'))
    parpool();
end

% Seed the RNG once – MATLAB will split streams across workers
rng('shuffle');

% ---------- Serial version ----------
serialStart = tic;
serialHits = 0;
for i = 1:N
    x = rand();
    y = rand();
    if x^2 + y^2 <= 1
        serialHits = serialHits + 1;
    end
end
serialTime = toc(serialStart);
piSerial = 4 * serialHits / N;

% ---------- Parallel version ----------
% Use an array to collect per‑iteration hits – this is a valid reduction pattern
parHits = zeros(1,N);
parfor i = 1:N
    x = rand();
    y = rand();
    if x^2 + y^2 <= 1
        parHits(i) = 1; % record a hit
    end
end
parallelTime = toc(serialStart) - serialTime; % elapsed after serial
piParallel = 4 * sum(parHits) / N;

% ----- Results -----
fprintf('Serial:   π ≈ %.8f (%.2fs)\n', piSerial, serialTime);
fprintf('Parallel: π ≈ %.8f (%.2fs)\n', piParallel, parallelTime);

Key points about the code:

• We allocate a full‑size array parHits and set each element to 1 or 0. MATLAB recognises this as a valid reduction (sum) and aggregates the results across workers.
• The rng('shuffle') call ensures each worker starts with a different random seed, eliminating correlation between streams.
• We measure the serial time first; the parallel time is taken from the same wall clock reference to avoid re‑timing the pool start.

Trade‑offs and Practical Limits

• Pool Overhead – Launching a pool can take 1–3 seconds. For N < 1e6 the overhead may dominate, making parfor slower than a serial loop.
• Memory Footprint – Each worker runs as a separate MATLAB process, but read‑only data is shared via the operating system’s copy‑on‑write mechanism. Still, large per‑worker arrays can consume significant RAM.
• Variable Dependencies – A parfor cannot read or write to a variable that depends on the previous iteration (e.g., A(i)=A(i-1)+1). If such a dependency is needed, restructure the algorithm or fall back to a serial loop.
• Random Number Streams – If you need reproducible results, set a fixed seed before the loop and optionally use parfor‑specific stream handling (e.g., parallel.defaultClusterProfile or parallel.pool.Constant).

Verifying Correctness and Performance

1. Statistical Accuracy

   Run the script with increasing N (e.g., 1e5, 1e6, 1e7). The π estimate should converge toward 3.14159. If the result is far off, check that the RNG is seeded correctly and that the conditional logic is identical in both loops.

2. Timing Comparison

   Use tic and toc as shown. A typical speedup on a quad‑core machine might be 2–3× for N = 1e7, but the exact factor depends on CPU, memory bandwidth, and MATLAB’s internal scheduling.

3. Reproducibility Check

   Set a fixed seed with rng(12345) before the parfor and run the script twice. The two π estimates should be identical (within floating‑point noise). If they differ, the workers are not getting independent streams.

Actionable Takeaway

Use parfor for Monte Carlo simulations when the number of iterations is large enough that the per‑iteration cost outweighs pool startup and data movement. Always verify that your loop is truly independent, that the RNG streams are correctly seeded, and that you have a Parallel Computing Toolbox license. With those safeguards, you can reduce run times from hours to minutes on a multi‑core desktop.