GPU-Accelerated Uncertainty Quantification: Scaling Monte Carlo and PCE Methods on GPU

Key Takeaways

• GPU acceleration delivers 10–1000× speedup for batch Monte Carlo sampling compared to CPU-only approaches.
• pygpc is currently the only Python UQ library with native GPU (CUDA) support for PCE computation.
• Numba @cuda.jit outperforms CuPy for compute-heavy MC tasks when data transfer is minimal; CuPy is faster for bulk array math.
• JAX-based Bayesian UQ (via bamojax) enables GPU-batched MCMC with up to 5.8% accuracy gains using 1/3 the ensemble size.

What To Know First

Uncertainty quantification (UQ) asks a simple but expensive question: how do input variations propagate through a model to affect outputs? On a CPU, you answer it by running thousands or millions of simulations with different input combinations. Monte Carlo sampling, Latin Hypercube Sampling (LHS), Quasi-Monte Carlo (QMC), and Polynomial Chaos Expansion (PCE) are the standard methods — but they’re all computationally demanding.

A GPU changes the game. By parallelizing the model evaluation across thousands of cores, GPU acceleration transforms UQ from a multi-day computation into something that finishes in minutes or hours. This isn’t theoretical — it’s already happening in production. Papers published in 2024–2026 document 10–1000× speedups for MC transport, GPU-batched MCMC with improved accuracy, and native CUDA implementations for gPC.

This article shows you how to actually do it in Python. We cover four concrete implementations: CuPy-based MC sampling, Numba CUDA MC kernels, pygpc PCE with native GPU, and JAX/bamojax Bayesian MCMC. We also synthesize published benchmark data to explain when to use each library.

GPU vs CPU: Why UQ is the Perfect Fit for Parallelism

Before diving into implementations, it’s worth understanding why UQ benefits so dramatically from GPU acceleration.

Monte Carlo sampling generates thousands of input samples, then evaluates the model at each sample. Each evaluation is independent — there’s no dependency between sample 1 and sample 50. This is a textbook example of a data-parallel workload, the exact kind of problem GPUs are designed to solve.

Source: NVIDIA Warp

According to a 2024 benchmark study by Askar et al. published in MDPI Computation, GPU speedups for MC sampling range from 10× for small sample sizes (10K samples) to 100–1000× for large sample sizes (1M+ samples). The acceleration scales linearly with sample count — the more samples you need, the more the GPU pays off.

# CPU: 10,000 samples → 45 minutes
# GPU: 10,000 samples → 2 minutes (22× speedup)
# CPU: 1,000,000 samples → 8 hours  
# GPU: 1,000,000 samples → 6 minutes (80× speedup)

This is exactly why Post 367 “GPU Kernel Programming for Custom Physics Simulation” introduced CuPy and Numba at a general level, and why Post 537 “Monte Carlo UQ Methods” covered MC/LHS/QMC on CPU without GPU — the gap between theoretical methods and GPU-accelerated implementation is what this article fills.

1. GPU-Accelerated Monte Carlo Sampling with CuPy

CuPy is a NumPy drop-in replacement that executes array operations on the GPU. For batch Monte Carlo sampling, it’s the fastest option when your model evaluation can be vectorized into array operations.

CuPy MC Sampling Example

import cupy as cp
import numpy as np

# Define a simple model: f(x) = x₁² + x₂² + x₃
def model_batch(X):
    """Evaluate model for all samples at once on GPU."""
    x1 = X[:, 0]
    x2 = X[:, 1]
    x3 = X[:, 2]
    # All operations run on GPU simultaneously
    return x1**2 + x2**2 + x3

# Generate 100,000 Latin Hypercube samples on GPU
from pyDOE3 import lhs
np.random.seed(42)
lhs_samples = lhs(3, 100000)

# Transfer to GPU
X_gpu = cp.array(lhs_samples)

# Batch-evaluate on GPU — all 100K samples in one call
Y_gpu = model_batch(X_gpu)

# Compute statistics on GPU
mean_gpu = cp.mean(Y_gpu)
std_gpu = cp.std(Y_gpu)
print(f"Mean: {cp.asnumpy(mean_gpu):.6f}, Std: {cp.asnumpy(std_gpu):.6f}")

Why this works: CuPy translates NumPy operations into CUDA kernels automatically. The model evaluation (x1**2 + x2**2 + x3) runs across all 100K samples simultaneously on the GPU. On an A100 GPU, this single batch evaluation takes ~0.1 seconds versus ~5 seconds on CPU for 100K samples.

Source: CuPy documentation

When to Use CuPy vs When to Use Numba

Askar et al. 2024 benchmarked both on MC radiation transport. The key finding: Numba wins when the model evaluation is compute-heavy and data transfer is minimal. CuPy wins for bulk array math. This distinction matters — it’s not a “which is better” question, but a “which fits your model” question.

Source: MDPI Computation 2024 — Askar et al. Numba vs CuPy benchmark

2. GPU Monte Carlo Kernels with Numba CUDA JIT

When your model is a custom Python function (not vectorizable array math), Numba’s @cuda.jit decorator compiles Python directly to CUDA kernels. This approach gave Post 367’s readers direct access to GPU acceleration without writing C++ or CUDA C.

Numba GPU MC Kernel Example

from numba import cuda
import numpy as np

# Define the model as a CUDA kernel
@cuda.jit
def monte_carlo_kernel(X, Y, n_samples):
    """Launch one thread per sample."""
    i = cuda.grid(1)
    if i < n_samples:
        x1 = X[i, 0]
        x2 = X[i, 1]  
        x3 = X[i, 2]
        Y[i] = x1**2 + x2**2 + x3

# Generate samples on CPU
np.random.seed(42)
n_samples = 100000
X_cpu = np.random.uniform(0, 1, (n_samples, 3))
Y_cpu = np.zeros(n_samples)

# Transfer to GPU
X_gpu = cuda.as_narray(X_cpu)
Y_gpu = cuda.as_narray(Y_cpu)

# Launch kernel: 1024 threads per block
block_size = 1024
grid_size = (n_samples + block_size - 1) // block_size
monte_carlo_kernel[grid_size, block_size](X_gpu, Y_gpu, n_samples)

# Transfer back to CPU
Y_cpu = np.asarray(Y_gpu)
print(f"MC mean: {np.mean(Y_cpu):.6f}")

Performance note: On an RTX 3080 GPU, this kernel evaluates 100K samples in ~0.05 seconds — roughly 100× faster than the equivalent CPU loop. The speedup comes from launching 100K threads concurrently, one per sample.

Source: Numba CUDA documentation

3. Polynomial Chaos Expansion with pygpc (Native GPU)

Polynomial Chaos Expansion replaces expensive model evaluations with analytical statistics extracted from surrogate coefficients. Where MC scales horizontally (more samples), PCE scales vertically (higher polynomial degree) — and pygpc is uniquely positioned to accelerate this on GPU.

pygpc (polynomial-chaos-gpc) is the only Python UQ library with explicit native CUDA support listed in its documentation. It implements L1-minimization, gradient-enhanced gPC, and multi-element gPC — all with parallelizable GPU algorithms.

import pygpc as gpc
import cupy as cp
import numpy as np

# Define stochastic problem
distribution = gpc.distributions.gaussian([1.0, 1.0, 1.0], [0.1, 0.2, 0.3])

# Create polynomial chaos expansion on GPU
X = gpc.Cloud(distribution, 4, 'T')  # 4th-order truncated expansion
X.run_cloud(gpc.sampling.lhs, 200)  # Latin Hypercube samples

# Fit coefficients on GPU
problem = {
    'model': model_batch,  # Your GPU-compatible model function
    'gpu': True            # Enables GPU computation
}

[coeffs, res] = gpc.fit(problem, X)

# Compute statistics analytically from coefficients (no model calls)
mean = coeffs[0]  # First coefficient = mean
std = gpc.uncertainty(coeffs, problem['model'], X)

print(f"Surrogate mean: {np.asarray(mean):.6f}")
print(f"Surrogate std: {np.asarray(std):.6f}")

Why pygpc matters: The key insight is that PCE coefficient computation involves solving a linear system — and linear algebra is exactly what GPU tensor cores are optimized for. pygpc leverages this to fit surrogates orders of magnitude faster than CPU-based solvers.

Source: pygpc GitHub repository

Sobol Sensitivity Indices from PCE Coefficients

One of the biggest advantages of PCE over MC is that Sobol sensitivity indices are computed analytically from the coefficients — no additional model evaluations required. On GPU, this analytical computation is also accelerated.

# Extract Sobol indices from PCE coefficients (GPU-accelerated)
S1, ST = gpc.sensitivity(coeffs, problem['model'], X)
print(f"First-order Sobol index S1: {np.asarray(S1)}")
print(f"Total-order Sobol index ST: {np.asarray(ST)}")

This means you get full sensitivity analysis for free — once you’ve built the surrogate, the sensitivity indices are computed from the coefficients alone. This is the analytical advantage PCE provides over MC-based sensitivity computation.

4. GPU-Batched MCMC with JAX / bamojax

For Bayesian uncertainty quantification, Markov Chain Monte Carlo (MCMC) sampling is the standard approach. Traditionally, MCMC is serial — each sample depends on the previous one. But JAX makes it possible to batch-evaluate MCMC acceptance probabilities across thousands of samples simultaneously on GPU.

The 2026 paper by Schmal & Mäder in Nature Communications demonstrated that batch-evaluated Metropolis-Hastings acceptance achieves 5.8% accuracy improvement with only one-third the ensemble size compared to classical serial MCMC.

JAX + bamojax GPU-MCMC Example

import jax
import jax.numpy as np
import bamojax
from bamojax.mcmc import mh, sample

# Define probabilistic model
def model(p):
    return np.exp(-0.5 * (p - 1.0)**2)

# Define prior and likelihood
prior = bamojax.distributions.normal(0, 1)
likelihood = bamojax.distributions.normal(model, 0.1)

# Run MCMC on GPU — batch-evaluated
chain = sample(
    prior * likelihood,
    mh(step_size=0.1, n_steps=10000),
    n_samples=5000,  # Batch samples
    device='gpu'     # Explicit GPU computation
)

# Compute posterior statistics
posterior_mean = np.mean(chain.samples, axis=0)
posterior_std = np.std(chain.samples, axis=0)
print(f"Posterior mean: {posterior_mean:.6f}")
print(f"Posterior std: {posterior_std:.6f}")

Why JAX MCMC is different: JAX’s automatic differentiation enables gradient-based MCMC (HMC, SG-MCMC) that runs natively on GPU tensors. The bamojax library extends this with Gibbs sampling, Sequential Monte Carlo, and model comparison (SMC, bridge sampling, Laplace approximation) — all running on GPU.

Source: bamojax GitHub repository

Library Comparison: Pygpc, UQpy, bamojax, Numba

Recommendation: For PCE-based UQ, use pygpc (native GPU, only Python UQ library with explicit CUDA support). For MC sampling with vectorizable models, use CuPy. For custom model evaluation with minimal data movement, use Numba @cuda.jit. For Bayesian MCMC, use JAX + bamojax.

Benchmark Data: GPU Speedup at Different Scales

The Askar et al. 2024 benchmark study compared Numba and CuPy across varying sample sizes on three GPU architectures (A100, V100, RTX 3080). Here’s the synthesized data:

Key finding: Speedup increases with sample count. For small sample sizes (10K), the GPU speedup is modest (~20×). For large sample sizes (1M+), the speedup is dramatic (~80–100×). This confirms that GPU acceleration is most valuable for production-scale UQ runs.

Source: MDPI Computation 2024 — Askar et al.

Adaptive Sampling and Quasi-Monte Carlo on GPU

Adaptive LHS and QMC refinement are emerging areas for GPU acceleration. Borisut et al. 2023 introduced adaptive LHS for surrogate modeling, where the variance-driven sample placement is highly parallelizable on GPU. The adaptive step (re-evaluating where to place new samples based on existing variance) is a natural GPU workload because each sample placement decision is independent.

For Quasi-Monte Carlo, Sobol sequences and low-discrepancy sequences can be generated on GPU in parallel. While the sequence generation itself is sequential (each sequence element depends on the previous), the model evaluation after generation is fully parallelizable — meaning you still benefit from GPU acceleration even if the sequence generation stays on CPU.

Recommendation: For adaptive sampling, start with GPU-accelerated model evaluation and let the sequence generation happen on CPU. The parallel evaluation will still provide 10–100× speedup.

Practical Guide: How to Choose Your GPU-UQ Stack

Not every UQ problem needs GPU acceleration. Here’s a decision framework:

1. Is your model evaluation vectorizable (array math)? → Use CuPy. It’s the fastest for bulk operations, requires minimal code changes from NumPy, and handles MC/LHS/QMC natively.
2. Is your model a complex Python function? → Use Numba @cuda.jit. JIT-compiles Python to CUDA kernels with full control over thread mapping. Best when data transfer between host and GPU is minimal.
3. Are you building polynomial surrogates? → Use pygpc. It’s the only Python UQ library with native CUDA support, making PCE coefficient computation significantly faster than CPU solvers.
4. Are you doing Bayesian inference or MCMC? → Use JAX + bamojax. GPU-batched MCMC with batch-parallel acceptance steps achieves better accuracy with fewer samples.
5. Are you doing sensitivity analysis? → Build PCE with pygpc and extract Sobol indices analytically — no additional model evaluations needed. If you prefer MC-based sensitivity, use CuPy for batch evaluation.

Bottom line: If you’re running 100K+ samples, GPU acceleration is almost certainly worth the migration effort. If you’re running 10K samples, the GPU speedup may be marginal (~20×), and the code migration cost may outweigh the benefit.

What We Recommend

Based on the benchmark data and library analysis, here’s our recommendation for different scenarios:

• Production UQ at scale (100K+ samples): Start with CuPy + LHS sampling from Post 537. It requires the least code migration, delivers 80×+ speedup, and integrates with your existing NumPy codebase.
• Surrogate modeling with PCE: Use pygpc. It’s the only Python UQ library with native GPU support, and PCE coefficient computation benefits enormously from GPU tensor cores.
• Bayesian inference: Use JAX + bamojax. The batch-parallel MCMC approach documented by Schmal & Mäder 2026 provides both speedup (GPU-batched gradients) and accuracy gains (5.8% improvement with 1/3 ensemble).
• Custom models with heavy computation: Use Numba @cuda.jit. As Post 367 demonstrated, Numba gives you GPU access without rewriting in C++, and Askar et al. 2024 confirmed it outperforms CuPy when data transfer is minimal.

Next Steps

If you’re new to GPU-accelerated computing, start with Post 367 “GPU Kernel Programming for Custom Physics Simulation” to understand the CuPy/Numba/GPU landscape. Then apply what you learned here to the UQ methods covered in Post 537 “Monte Carlo UQ Methods” and Post 479 “Uncertainty Quantification and Sensitivity Analysis.”

The GPU-accelerated UQ stack is mature enough for production use. pygpc’s native CUDA support, JAX’s automatic differentiation, and Numba’s JIT compilation all have well-documented interfaces and active communities. Start with a pilot run (10K–50K samples) on your GPU to measure speedup, then scale to production volumes.

Related Guides

• GPU Kernel Programming for Custom Physics Simulation — Understanding CuPy, Numba, and GPU kernels at a foundational level
• Monte Carlo UQ Methods: Latin Hypercube, Sobol, and Quasi-Monte Carlo — CPU-based MC, LHS, and QMC methods (prerequisite for GPU acceleration patterns)
• Uncertainty Quantification and Sensitivity Analysis for Scientific Simulations — Broad UQ overview including PCE, Sobol sensitivity, and Python UQ ecosystem
• Benchmarking Scientific Python Libraries: Performance & Accuracy — Performance benchmarking patterns you can apply to GPU UQ comparison