Probability

Language revisions

C++20

The <random> engines and distributions are stable and reproducible; std::uniform_random_bit_generator is a concept.

C++23

std::generator supports lazy sample streams without an allocation per sample.P2502R2

C++26

std::simd enables vectorized generation and evaluation of independent streams.P1928

Definition

Probability spaces and random variables

DefinitionProbability space

A probability space is a triple $(\Omega, \mathcal{F}, \mathbb{P})$ with $\Omega$ a set, $\mathcal{F}$ a $\sigma$-algebra on it, and $\mathbb{P} : \mathcal{F} \to [0,1]$ a countably additive measure with $\mathbb{P}(\Omega) = 1$. A random variable is a measurable map $X : \Omega \to \mathbb{R}$, and its expectation is the integral $\mathbb{E}[X] = \int_\Omega X \, d\mathbb{P}$ when it exists.

Why the measure-theoretic setting is worth keeping

For finite sample spaces the machinery is unnecessary, and for the continuous ones that occur in simulation it is not: without it there is no coherent account of conditioning, of limits of sequences of random variables, or of the stochastic integrals that stochastic differential equations are defined by. The setting also supplies the notion of a $\sigma$-algebra as the information available at a point in time, which is what makes filtrations and hence path-dependent simulation precise.

Variance and its additivity

The variance is $\operatorname{Var}(X) = \mathbb{E}[(X - \mathbb{E}X)^2]$, and for independent $X_1, \dots, X_n$ variances add. That single fact is what produces the $n^{-1/2}$ rate that governs every Monte Carlo method, and it is also what makes variance reduction the only lever available: the rate is fixed, so only the constant $\sigma^2$ can be changed.

Independence and its computational meaning

Independence is a statement about the measure. In computation it becomes a statement about streams: two sequences of pseudorandom numbers are treated as independent when no statistical test in use distinguishes them from independent samples. This is a weaker and more fragile notion than the mathematical one, and it is why parallel Monte Carlo requires deliberate stream splitting rather than sharing one generator.

Modes of convergence

Three modes are used, and they are not interchangeable.

In probability

$\mathbb{P}(|X_n - X| > \varepsilon) \to 0$ for every $\varepsilon$. This is what the weak law of large numbers gives, and it is enough to justify a single estimate at a single sample size.

Almost surely

$\mathbb{P}(X_n \to X) = 1$. This is what the strong law gives, and it is what justifies treating one long run as representative rather than needing independent repetitions.

In distribution

$\mathbb{E}[g(X_n)] \to \mathbb{E}[g(X)]$ for bounded continuous $g$. This is what the central limit theorem gives, and it is the basis of confidence intervals.

TheoremCentral limit theorem

For independent identically distributed $X_i$ with mean $\mu$ and finite variance $\sigma^2$,

$$\sqrt{n}\,\frac{\bar{X}_n - \mu}{\sigma} \;\xrightarrow{d}\; \mathcal{N}(0, 1).$$

The finite-variance hypothesis is where practical trouble enters: for heavy-tailed integrands the sample variance does not converge, the confidence interval computed from it is meaningless, and the failure is silent. The Berry-Esseen theorem quantifies the rate of approach, $O(n^{-1/2})$ with a constant involving the third absolute moment, which is why intervals from small samples of skewed distributions are unreliable even when the variance is finite.

Distributions that recur

Distribution	Where it appears	Sampling method
Uniform on $[0,1)$	The primitive from which others are built	Bit manipulation of the generator output
Normal	Errors, diffusions, random projections	Box-Muller, Ziggurat
Exponential	Waiting times, Poisson processes	Inverse transform, $-\log U$
Poisson	Event counts	Knuth's product method, or inversion
Discrete with $k$ outcomes	Resampling, mixtures	Alias method, $O(1)$ per sample
Multivariate normal	Gaussian processes, sketching	Cholesky factor times a normal vector

The alias method is worth naming because it is the case where the obvious algorithm, a linear or binary search over a cumulative distribution, is asymptotically worse than necessary: $O(k)$ or $O(\log k)$ per sample against $O(1)$ after $O(k)$ preprocessing.

Motivation

Randomness enters numerical computing in three distinct ways, and conflating them causes errors.

As a model of the problem

Stochastic differential equations, Monte Carlo integration of high-dimensional integrals, and uncertainty quantification. Here the randomness is part of the specification and the answer is a distribution or a functional of one. Reproducibility is a testing convenience, not a correctness requirement.

As an algorithmic device on a deterministic problem

Randomised low-rank approximation, sketching for least squares, randomised pivoting, hashing. Here the input is deterministic and the algorithm's guarantees are probabilistic, holding with high probability over its internal randomness.

The distinction matters for testing: such an algorithm may be rerun with a different seed, and a failure that does not recur is a probabilistic event rather than a defect. It also matters for reporting: a result without its failure probability is as incomplete as a numerical result without its error bound.

As an error model

Rounding errors treated as independent random variables give $\sqrt{n}u$ estimates instead of the worst-case $nu$ bounds of Numerical Analysis. This use is heuristic, since rounding errors are deterministic functions of the data, but it predicts observed behaviour far better than the worst-case bound, and it is the reason worst-case bounds are used to establish scaling rather than to predict magnitudes.

Formal formulation

Monte Carlo integration

Monte Carlo integration estimates $I = \mathbb{E}[f(X)]$ by

$$\hat{I}_n = \frac{1}{n} \sum_{i=1}^{n} f(X_i), \qquad \mathbb{E}[\hat{I}_n] = I, \qquad \operatorname{Var}(\hat{I}_n) = \frac{\sigma^2}{n},$$

so the root-mean-square error is $\sigma / \sqrt{n}$, independent of dimension.

The dimensional comparison

Against a tensor-product rule of order $p$ in $d$ dimensions, which needs $N = \Theta(\varepsilon^{-d/p})$ points:

$d$	Tensor product, $p = 4$	Monte Carlo
2	$\varepsilon^{-1/2}$	$\varepsilon^{-2}$
4	$\varepsilon^{-1}$	$\varepsilon^{-2}$
8	$\varepsilon^{-2}$	$\varepsilon^{-2}$
16	$\varepsilon^{-4}$	$\varepsilon^{-2}$
100	$\varepsilon^{-25}$	$\varepsilon^{-2}$

The crossover at $d = 2p$ is the whole case for Monte Carlo, and it explains why the method is standard in finance and statistical physics and rare in one-dimensional quadrature.

Quasi-Monte Carlo

Low-discrepancy sequences replace random points with deterministic ones that fill the space more evenly, giving error $O(N^{-1}(\log N)^d)$ by the Koksma-Hlawka inequality for integrands of bounded variation in the sense of Hardy and Krause. The improvement is real for moderate $d$ and smooth integrands, and it disappears as $d$ grows because of the $(\log N)^d$ factor.

Randomised quasi-Monte Carlo, which applies a random shift to the sequence, recovers an unbiased estimator with a computable error estimate while keeping most of the improvement. Without the randomisation there is no error estimate at all, which is the practical reason it is used.

Concentration inequalities

Concentration inequalities convert variance into probabilistic bounds, and they are what turn a Monte Carlo estimate into a statement with a confidence level.

The basic bounds

For independent $X_i$ with $|X_i| \le M$ and $\sum \operatorname{Var}(X_i) = v$, Bernstein's inequality gives

$$\mathbb{P}\left( \left| \sum_i (X_i - \mathbb{E}X_i) \right| \ge t \right) \le 2 \exp\left( \frac{-t^2/2}{v + Mt/3} \right),$$

and Hoeffding's inequality for $X_i \in [a_i, b_i]$ gives

$$\mathbb{P}\left( \left| \hat{I}_n - I \right| \ge t \right) \le 2\exp\left( \frac{-2n^2t^2}{\sum_i (b_i - a_i)^2} \right).$$

Bernstein's bound is the sharper of the two when the variance is small relative to the range, which is the usual case for a well-designed estimator. It is what makes variance reduction pay twice: once in the estimate and once in the width of the interval.

Why they beat Chebyshev

Chebyshev's inequality gives a bound decaying polynomially in $t$, and the exponential bounds above decay exponentially. The difference is what makes high-confidence statements affordable: to reach failure probability $\delta$, Chebyshev needs $n = \Omega(1/\delta)$ samples and Hoeffding needs $n = \Omega(\log(1/\delta))$. Confidence of $10^{-9}$ costs a billion samples under the first and a factor of about twenty under the second.

Randomised numerical linear algebra

For matrices, the corresponding tool is the matrix Bernstein inequality, which gives randomised low-rank approximation its guarantee.

TheoremRandomised range finder

Let $A \in \mathbb{R}^{m \times n}$ with singular values $\sigma_1 \ge \sigma_2 \ge \cdots$, let $\Omega \in \mathbb{R}^{n \times (k+p)}$ have independent standard Gaussian entries with oversampling $p \ge 2$, and let $Q$ be an orthonormal basis for the range of $A\Omega$. Then

$$\mathbb{E} \lVert A - QQ^{\mathsf{T}}A \rVert_F \le \left(1 + \frac{k}{p-1}\right)^{1/2} \Big( \sum_{j > k} \sigma_j^2 \Big)^{1/2}.$$

Why the structural gain exceeds the asymptotic one

The optimal rank-$k$ error is $\big(\sum_{j>k}\sigma_j^2\big)^{1/2}$ by the Eckart and Young theorem of Linear Algebra, so the randomised method is within a modest factor at cost $\Theta(mn(k+p))$ rather than $\Theta(mn^2)$.

The more important gain is structural. The randomised algorithm touches $A$ only through matrix-matrix products, so it is compute-bound and parallel, whereas a classical rank-revealing factorization is a sequence of dependent level-2 operations. It also requires only one or two passes over $A$, which makes it applicable when $A$ does not fit in memory at all.

Power iteration for slowly decaying spectra

When the singular values decay slowly the basic method is inaccurate, and the standard remedy replaces $A\Omega$ by $(AA^{\mathsf{T}})^q A\Omega$, which raises the singular values to the power $2q+1$ and separates them. Two or three power iterations suffice in practice, each costing two more passes over $A$, and the orthonormalisation between iterations is required because the powers otherwise collapse numerically onto the dominant direction.

Random projections

TheoremJohnson and Lindenstrauss

For any set of $N$ points in $\mathbb{R}^d$ and any $\varepsilon \in (0,1)$, there is a map into $\mathbb{R}^k$ with $k = O(\varepsilon^{-2}\log N)$ preserving all pairwise distances to within a factor $1 \pm \varepsilon$, and a random Gaussian or sign matrix supplies one with high probability.

The dimension of the target depends on the number of points and not on the original dimension, which is the property that makes sketching viable. The same construction underlies sketched least squares, where an $m \times n$ problem is replaced by an $O(n/\varepsilon^2) \times n$ one at the cost of a controlled loss in the residual.

C++ representation

Generators and distributions

#include <cmath>
#include <concepts>
#include <cstddef>
#include <random>

namespace hpc {

template <class G>
concept RandomBitGenerator = std::uniform_random_bit_generator<G>;

// A distribution over T: sampling requires a generator, never a global one.
template <class D, class G>
concept Distribution =
    RandomBitGenerator<G> &&
    requires(D& d, G& g) {
        typename D::result_type;
        { d(g) } -> std::same_as<typename D::result_type>;
    };

// Semantic requirements:
//   S1. Successive calls to d(g) are independent given g's state.
//   S2. The generator is not shared between threads without synchronisation.
//   S3. Reproducibility requires the seed, the engine, and the exact sequence
//       of calls, so any change to the call pattern changes the results.

}  // namespace hpc

The fragility of reproducibility

The third obligation is the one most often violated. Adding a single extra sample inside a loop, even one whose value is discarded, changes every subsequent value and therefore every result. This makes seemingly harmless refactoring a source of unreproducible output, and it is the reason counter-based generators, whose output is a pure function of a key and a counter, are preferred where reproducibility matters.

Distributions are not portable

The standard specifies the engines exactly and the distributions only by their statistical properties, so std::normal_distribution may produce different values on different implementations from the same seed. Code that must reproduce results across toolchains implements its own transformations from the uniform generator rather than relying on the library distributions.

Estimates carry their uncertainty

#include <cmath>
#include <cstddef>

namespace hpc {

struct estimate {
    double mean;
    double standard_error;
    std::size_t samples;

    // Asymptotic normal interval. Valid only when the integrand has finite
    // variance and n is large enough for the central limit theorem to apply.
    constexpr double half_width(double z = 1.96) const noexcept {
        return z * standard_error;
    }
};

template <class F, RandomBitGenerator G>
estimate monte_carlo(F&& integrand, G& generator, std::size_t n) {
    double mean = 0.0;
    double m2 = 0.0;
    for (std::size_t i = 1; i <= n; ++i) {
        const double x = integrand(generator);
        const double delta = x - mean;
        mean += delta / static_cast<double>(i);
        m2 += delta * (x - mean);
    }
    const double variance = n > 1 ? m2 / static_cast<double>(n - 1) : 0.0;
    return {mean, std::sqrt(variance / static_cast<double>(n)), n};
}

}  // namespace hpc

Why Welford rather than sums of powers

Welford's recurrence is used rather than accumulating $\sum x$ and $\sum x^2$, because the latter subtracts two large nearly equal quantities and loses all significance when the mean is large relative to the standard deviation. For $x$ near $10^6$ with standard deviation $1$, the sum of squares computation in binary64 retains about four digits of the variance; Welford's retains all of them. This is a cancellation problem of exactly the kind analysed in Numerical Analysis, appearing in a routine that looks like bookkeeping rather than numerics.

Implementation

Independent, reproducible streams

Parallel Monte Carlo requires streams that are independent and reproducible. Splitting a single generator by seeding from a counter is the reliable construction; sharing one engine across threads is a data race, and seeding from the clock destroys reproducibility.

#include <cstdint>
#include <random>
#include <vector>

namespace hpc {

// One independent, reproducible stream per worker, derived from a single root
// seed by a sequence-dependent seeding scheme.
inline std::vector<std::mt19937_64> split_streams(std::uint64_t root_seed,
                                                  std::size_t workers) {
    std::vector<std::mt19937_64> engines;
    engines.reserve(workers);
    for (std::size_t w = 0; w < workers; ++w) {
        std::seed_seq seq{root_seed, static_cast<std::uint64_t>(w)};
        engines.emplace_back(seq);
    }
    return engines;
}

}  // namespace hpc

Counter-based generation

Counter-based generators such as Philox are the better construction where they are available: the stream is a pure function of a key and a counter, so independence is by construction rather than by a seeding heuristic, there is no state to displace cache, the generator vectorises, and a simulation can be restarted at an arbitrary point without replaying the sequence. The last property is what makes checkpointing a stochastic simulation practical.

Combining partial estimates

Reduction of per-worker estimates must combine means and variances rather than concatenating samples, using the parallel form of Welford's update:

$$m = \frac{n_a m_a + n_b m_b}{n_a + n_b}, \qquad M_2 = M_{2,a} + M_{2,b} + (m_b - m_a)^2 \frac{n_a n_b}{n_a + n_b}.$$

namespace hpc {

// Associative and numerically stable: the combination is order-independent up
// to rounding, so it may be applied in any reduction tree.
constexpr estimate combine(const estimate& a, const estimate& b) noexcept;

}  // namespace hpc

The associativity is what allows the combination to be used as the operator of a parallel reduction in the sense of Algebra, so a Monte Carlo run over $p$ workers is a monoid reduction over partial estimates.

Variance reduction

Variance reduction changes the constant in a rate that cannot be improved, and it is therefore worth more than any implementation tuning.

Control variates

Subtract a correlated quantity with known mean: $\hat{I} = \overline{f(X)} - c\big(\overline{g(X)} - \mathbb{E}g\big)$, with optimal $c = \operatorname{Cov}(f,g)/\operatorname{Var}(g)$, reducing the variance by the factor $1 - \rho^2$. A correlation of $0.95$ is a tenfold reduction in variance and hence a hundredfold reduction in samples.

Antithetic variates

Pair each sample with its reflection $1 - U$, which removes the odd component of the integrand exactly and cannot increase the variance for monotone integrands. The cost is nothing beyond the pairing.

Importance sampling

Sample from a distribution concentrated where the integrand is large, weighting by the likelihood ratio. The variance can be reduced by orders of magnitude and can also be made infinite by a poor choice of proposal, which is why the weights are monitored: an effective sample size far below the nominal one is the signature of a bad proposal.

Stratification

Partition the domain and allocate samples per stratum. It cannot increase the variance, and with proportional allocation it removes the between-stratum component entirely, which for smooth integrands is most of it.

Complexity

Complexity

Operation	Time	Space	Remarks
Monte Carlo to RMSE ε	Θ(σ²/ε²) samples	Θ(1)	Dimension-independent
Quasi-Monte Carlo, smooth integrand	Θ(ε^{−1}(log n)^d)	Θ(d)	Degrades with dimension d
Tensor-product quadrature	Θ(ε^{−d/p})	Θ(1)	p = order; infeasible for d above about 6
Randomised rank-k SVD	Θ(mn(k+p))	Θ(m(k+p))	Against Θ(mn²) for a full SVD
Randomised SVD with q power iterations	Θ((2q+1)mn(k+p))	Θ(m(k+p))	Needed for slow spectral decay
Johnson–Lindenstrauss projection	Θ(Ndk)	Θ(Nk)	k = O(ε⁻² log N)
Welford update	Θ(1)	Θ(1)	Numerically stable
Sum-of-squares variance	Θ(1)	Θ(1)	Catastrophic cancellation
Mersenne Twister sample	Θ(1), large state	2.5 KB per stream	Displaces cache
Counter-based sample	Θ(1), no state	Θ(1)	Vectorises; restartable
Alias-method discrete sample	Θ(1)	Θ(k)	After Θ(k) preprocessing
Inverse-CDF discrete sample	Θ(log k)	Θ(k)	Binary search

Performance considerations

Generation often dominates

Random number generation is frequently the bottleneck in Monte Carlo codes, not the integrand. std::mt19937_64 produces one 64-bit value per call with a 2.5 KB state that displaces useful data from cache. Counter-based generators produce a value from a key and a counter with no state, vectorise across lanes, and give reproducible streams by construction.

Accuracy is expensive, so variance reduction is the lever

The $n^{-1/2}$ rate means one more digit costs a hundredfold increase in samples. Variance reduction changes $\sigma^2$ and is therefore worth more than any implementation tuning, since it improves the constant in a rate that is fixed for the method class.

Scaling results from Monte Carlo prove little

Monte Carlo parallelises with essentially no communication: each worker accumulates its own partial estimate and the reduction is $\Theta(p)$ at the end. This makes it the standard demonstration of near-linear scaling, and also the reason such demonstrations say little about a machine's behaviour on communication-bound problems.

Sampling non-uniform distributions

Sampling from non-uniform distributions is dominated by rejection rates and by transcendental function calls. The Ziggurat method for normal variates avoids both in the common case, at the cost of a table and a rare branch that breaks vectorization when it is taken. For vectorized code the Box-Muller transform is often preferable despite being slower scalar, because it is branch-free and therefore uniform across lanes.

Divergence across lanes

Any rejection method diverges: some lanes accept and some do not, so a vectorized implementation either masks and retries, wasting work, or falls back to scalar. Methods with a fixed operation count per sample are worth a substantial factor in vectorized code even when they are slower in isolation, and this is the same lane-divergence argument that appears in SIMD.

Limitations

Pseudorandom generators are deterministic. Statistical tests can detect bad ones but cannot certify a good one, and a generator adequate for simulation may be inadequate where adversarial inputs are possible.

Confidence intervals from the central limit theorem are asymptotic and are unreliable for heavy-tailed integrands, where the sample variance itself may not converge. Checking the tail behaviour of the integrand is part of using the method, not an optional refinement.

Probabilistic guarantees for randomised algorithms hold with a stated failure probability, which is a different kind of statement from a deterministic bound and must be reported as such.

Reproducibility across machines is not automatic even with a fixed seed: floating-point summation order, thread counts, vector widths and library distribution implementations all affect the result. Where bitwise reproducibility is required, the reduction tree must be fixed and the distributions implemented locally.

Importance sampling can increase the variance without any diagnostic unless the weights are monitored, and the failure is silent in exactly the regime where the method was expected to help.

Related entries

Numerical Analysis for the statistical error model; Numerical Integration for the deterministic alternatives; Linear Algebra for the approximation theorem behind randomised factorization; Algebra for the monoid structure of estimate combination; Parallelism for the reduction structure.

References

• R. Durrett, Probability: Theory and Examples, 5th ed., Cambridge University Press, 2019.
• N. Halko, P.-G. Martinsson and J. Tropp, Finding Structure with Randomness, SIAM Review 53(2), 2011.
• J. Tropp, An Introduction to Matrix Concentration Inequalities, Foundations and Trends in Machine Learning 8, 2015.
• B. Welford, Note on a Method for Calculating Corrected Sums of Squares and Products, Technometrics 4(3), 1962.
• J. Salmon, M. Moraes, R. Dror and D. Shaw, Parallel Random Numbers: As Easy as 1, 2, 3, SC11.
• A. Owen, Monte Carlo theory, methods and examples, 2013.
• G. Marsaglia and W. Tsang, The Ziggurat Method for Generating Random Variables, Journal of Statistical Software 5(8), 2000.