Determinism, Randomness, and Ensemble Runs

“Reproducible” does not always mean “the output files are identical.” A serial deterministic update on fixed hardware may reasonably be tested bit for bit. A parallel reduction may be correct when it agrees only within a numerical tolerance. A Monte Carlo calculation is expected to produce different samples and must be compared through a declared estimand and sampling distribution. Trust begins by choosing the right equivalence relation before the rerun, then recording enough environment and random-stream information to test it.

Required background. Calculus and Linear Algebra for Relativity supplies norms, orthogonal modes, and conserved quadratic forms; Differential Equations for Relativity supplies time evolution and convergence tests.

Helpful background. Formulations of General Relativity connects evolved fields to their physical constraints; Probability and Inference for Gravity supplies estimands, sampling errors, and interval coverage.

Three meanings of repeated agreement

A deterministic mathematical problem can have a nondeterministic implementation. Thread scheduling changes the order of floating-point sums; accelerator kernels may use unordered atomic updates; adaptive mesh or solver decisions may branch on last-bit differences. Conversely, a program can be bitwise repeatable while computing the wrong equation. Determinism is an implementation property, not a proof of accuracy.

Choose one of three comparison levels for each product:

Level	Appropriate object	Required comparison
bitwise	canonical inputs, exact integer streams, fixed-platform regression products	identical bytes and digest
numerical	deterministic fields or observables with legal operation-order variation	normed discrepancy below a justified tolerance
statistical	random ensembles or stochastic algorithms	estimand, uncertainty, coverage, and distributional agreement

The level can change along a workflow. A random generator’s raw integer output can be bitwise defined, each evolved field can be compared numerically, and the ensemble conclusion can be statistical. Calling the entire pipeline merely “deterministic” or “random” loses this structure.

For a numerical observable $Q$, a tolerance test should scale with its error budget, for example

$$\frac{|Q_A-Q_B|}{Q_{\rm scale}} \le \tau_{\rm impl} \ll \tau_{\rm discretization}+\tau_{\rm model}.$$

The scale and tolerances must be stated. “Close” is not an equivalence relation, and matching printed digits discards information about both magnitude and uncertainty.

Random streams are indexed scientific inputs

A single seed is not a complete stream specification. It omits the generator and version, seed-to-state transformation, mapping of draws to ensemble members, distribution algorithm, and draw order. Adding one diagnostic draw can shift every subsequent value in a stateful stream. Dividing one shared generator among threads makes the mapping depend on scheduling.

A counter-based generator instead evaluates a keyed function

$$R_{m,q,s}=b\!\left(K(m,q,s),C(m,q,s)\right),$$

where immutable identifiers—such as ensemble member $m$, Fourier mode $q$, and sample component $s$—determine the key $K$ and counter $C$. No mutable global state is required. Salmon and collaborators develop this construction, its period and stream separation, and its statistical tests in Salmon et al. 2011, § 2, §§ 2.1–2.2, and § 4. The identifiers and generator specification must be preserved beside the physical input.

This construction makes scheduling independence testable: generate the same set of $(m,q,s)$ tuples with one thread, two threads, and a shuffled task order. The raw integers must be identical for every tuple. It does not make the generator physically appropriate, ensure independence for a flawed key design, or validate the probability law. Those remain separate checks.

Even a widely used library may intentionally change algorithms across releases. NumPy’s revised policy explains why a fixed seed does not guarantee an unchanged stream across versions, builds, or machines in NumPy Developers 2019, NEP 19, §§ “Status quo” and “Proposal”. Preserve the bit generator, library version, distribution routine, parameters, and call mapping; archive raw variates when they are small enough and scientifically consequential.

A coupled scalar-wave ensemble

Use the periodic wave equation on $0\le x<2\pi$ and define member $m$ by

$$u_m(0,x)=\sin(2x) +\varepsilon\sum_{q=1}^{16} \left[a_{mq}\sin(qx)+b_{mq}\cos(qx)\right], \qquad \Pi_m(0,x)=-\partial_xu_m(0,x), \quad \Phi_m(0,x)=\partial_xu_m(0,x),$$

with $\partial_tu_m=\Pi_m$, $\varepsilon=10^{-4}$, and independent standard-normal $a_{mq}$ and $b_{mq}$. Generate each coefficient from a counter keyed by the immutable experiment ID, $m$, $q$, and a sine/cosine label. Because the random field is defined spectrally rather than at grid points, the same member exists at $N=64,128,256$; adding resolution does not silently change the continuum input.

For this right-moving exact solution, the continuum energy is

$$\mathscr H_m =4\pi+8\pi\varepsilon a_{m2} +\pi\varepsilon^2\sum_{q=1}^{16}q^2(a_{mq}^2+b_{mq}^2),$$

and orthogonality makes the ensemble mean

$$\mathbb E[\mathscr H] =4\pi+2\pi\varepsilon^2\sum_{q=1}^{16}q^2 =4\pi+2992\pi\varepsilon^2.$$

The cross term survives only for $q=2$, because that is the deterministic wave’s wavenumber. It averages to zero, but it dominates the fluctuations at small $\varepsilon$. In fact,

$$\operatorname{Var}(\mathscr H) =64\pi^2\varepsilon^2 +4\pi^2\varepsilon^4\sum_{q=1}^{16}q^4.$$

For $\varepsilon=10^{-4}$, the standard deviation is approximately $8\pi\varepsilon=0.00251$, while the increase in the mean above $4\pi$ is only $9.40\times10^{-5}$. Thus even $256$ members need not resolve that small increase: the standard error is approximately $1.57\times10^{-4}$. A noisy mean is not an energy-conservation failure.

This supplies analytic mean and variance benchmarks independent of the evolution code. Evolve $n=32,64,128,256$ members on $N=64,128,256$ point grids with the fourth-order centered stencil, classical fourth-order Runge–Kutta, maximum Courant factor $1/4$, and final time $t=1$, storing the input digest and coefficient-stream identifier with every field. Sampling the continuum $\Phi_m$ makes $\Phi_m-D_4u_m$ nonzero initially by $O(\Delta x^4)$; the semidiscrete equations preserve that known residual, which must converge rather than being compared with zero. At each $n$ and grid resolution, report the sample mean $\bar H_n$, sample standard deviation $s_H$, Monte Carlo standard error $s_H/\sqrt n$, this field-constraint residual, and discretization error relative to the member’s exact continuum Fourier evolution.

The predeclared checks are:

• raw coefficient integers agree bitwise across one, two, and eight threads and across a shuffled member schedule;
• each member approaches fourth-order field convergence until roundoff or time error dominates;
• the energy sample mean is consistent with the analytic expectation within three estimated standard errors, with the actual standardized discrepancy reported;
• doubling $n$ reduces the typical Monte Carlo standard error approximately as $n^{-1/2}$;
• a clean rerun from the archived environment, source, checksum-bearing inputs, and stream map reproduces raw variates bitwise and numerical fields at the declared level.

The environment record includes operating system or container, compiler or interpreter, numerical libraries, precision, thread and process counts, hardware, generator name and version, key/counter encoding, distribution transform, source revision, and every input checksum. Repeat the energy reduction with a reproducible accumulation or high-precision summation as an independent numerical check.

Ensembles test an estimand, not a collection of plots

Before sampling, specify the population, estimand, stopping rule, and interval method. The ensemble mean of energy answers a different question from a tail probability, maximum constraint violation, or waveform credible band. Increasing the number of members until an interval happens to exclude an inconvenient value invalidates its nominal interpretation.

Seeds should be assigned through stable member identities, never recycled because a run failed. If member 17 is retried, it retains the same random input and receives a new execution record. If a genuinely independent replacement is required, it receives a new identity and the exclusion rule is recorded. Record seeds and stream identities as recommended by Sandve et al. 2013, pp. 1–4, Rule 6, but also retain the generator and mapping details needed to interpret them.

Randomness introduced only by a numerical method must not be confused with physical uncertainty. An ensemble over stochastic rounding estimates sensitivity to that arithmetic scheme. An ensemble over uncertain masses or equations of state represents a declared input distribution. Neither automatically estimates model-form error.

Failure modes and rerun sequence

The primary failure modes are a shared stateful generator whose assignment changes with scheduling, seed reuse across supposedly independent members, an unspecified generator or distribution transform, grid-dependent noise that changes under refinement, a statistical claim based on bitwise identity, and a deterministic claim assessed only by overlapping error bars. Also reject an ensemble that reports mean behavior while hiding rare constraint or positivity failures.

A fully specified rerun first verifies the environment and canonical input digest; generates and hashes the indexed coefficient table; checks its schedule independence; evolves every member at all resolutions; computes memberwise convergence, constraints, and energy with an independent analysis path; forms the $n$-sequence of estimators; applies the bitwise, numerical, and statistical criteria; and repeats from a clean directory. It emits machine-readable discrepancies and failed criteria rather than only final figures. The scalar-wave benchmark is linear, smooth, and Gaussian-driven; it does not validate turbulent, chaotic, heavy-tailed, or adaptive production ensembles.

Check your understanding

Four threads draw from one locked, stateful generator, so the sequence of generated numbers is always the same. The operating system assigns completed draws to ensemble members in whichever order threads acquire the lock. Is the ensemble member map reproducible?

Solution

No. The global sequence may be repeatable while the association between a draw and member depends on scheduling. Give each member and sample coordinate an immutable counter or an independently derived stream key. Then the random input is a function of scientific identifiers rather than execution order, and the one-thread and shuffled-schedule tests can verify the mapping bitwise.

References

• NumPy Developers. “NEP 19 — Random Number Generator Policy.” NumPy Enhancement Proposals, 2019. Specification. Current compatibility policy.
• Salmon, John K., Mark A. Moraes, Ron O. Dror, and David E. Shaw. “Parallel Random Numbers: As Easy as 1, 2, 3.” In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. ACM, 2011. Open PDF. DOI.
• Sandve, Geir Kjetil, Anton Nekrutenko, James Taylor, and Eivind Hovig. “Ten Simple Rules for Reproducible Computational Research.” PLOS Computational Biology 9, no. 10 (2013): e1003285. DOI.