Multi-Order Monte Carlo Methods

• Multi-Order Monte Carlo is a hierarchical Monte Carlo framework that builds levels by varying spatial and temporal discretization orders on a fixed mesh instead of using mesh refinements.
• It leverages control variates by treating lower-order approximations as variance reducers for higher-order solvers while integrating IMEX Runge-Kutta schemes for stiff systems.
• The method effectively reduces both sampling error and variance while maintaining asymptotic consistency, making it valuable for uncertainty quantification in multiscale hyperbolic models.

Searching arXiv for recent and foundational papers on Multi-Order Monte Carlo and closely related multilevel/multi-index Monte Carlo methods. Multi-Order Monte Carlo (MOMC) is a hierarchical Monte Carlo framework for uncertainty quantification in time-dependent partial differential equations in which the hierarchy is built by changing the spatial and temporal discretization orders of accuracy, rather than by constructing a mesh-refinement ladder. In the formulation introduced for multiscale hyperbolic systems, MOMC uses numerical solvers of increasing order on the same mesh, interprets lower-order approximations as control variates for higher-order ones, and is designed to work naturally with Implicit-Explicit Runge-Kutta schemes that preserve asymptotic limits in stiff regimes. This places MOMC within the broader family of hierarchical Monte Carlo methods while giving it a distinct operational role in settings where hierarchical re-meshing is costly or awkward (Bertaglia et al., 27 Aug 2025).

1. Position within hierarchical Monte Carlo methods

The immediate background to MOMC is the multilevel Monte Carlo (MLMC) philosophy. In MLMC, the expectation of a finest-level approximation is written through a telescoping representation,

$$\mathbb{E}[\widehat{P}_L] = \mathbb{E}[\widehat{P}_0] + \sum_{l=1}^L \mathbb{E}[\widehat{P}_l - \widehat{P}_{l-1}],$$

so that many cheap coarse approximations and fewer expensive fine approximations can be combined to reduce computational cost for a prescribed accuracy (Higham, 2015). In PDE settings, this same principle has been adapted to finite element hierarchies, random coefficients, optimal control, ensemble time-stepping, and other structured discretizations (Sun et al., 2015, Luo et al., 2018).

Multi-Index Monte Carlo (MIMC) extends the MLMC idea from a single discretization axis to several axes simultaneously. Its defining feature is the use of high-order mixed differences rather than first-order differences, together with structured index sets such as full tensor or total degree sets; under standard mixed-regularity assumptions, this can yield improved complexity and memory scaling in high dimensions (Haji-Ali et al., 2014). Dimension-adaptive MIMC further replaces a priori index-set design with profit-based adaptive refinement across discretization directions such as space and Karhunen–Loève truncation levels (Robbe et al., 2017).

Within this lineage, MOMC differs in a specific way: its levels are orders of accuracy rather than meshes. The hierarchy is therefore constructed by varying both spatial and temporal discretization orders on a fixed mesh, and the estimator acts as a generalized control-variate construction between adjacent numerical orders (Bertaglia et al., 27 Aug 2025).

Method	Hierarchy	Characteristic mechanism
MLMC	Mesh refinements	Telescoping sum of level differences
MIMC	Multi-index discretization directions	High-order mixed differences
MOMC	Discretization order on the same mesh	Control variates between solver orders

This distinction matters most when the dominant numerical burden comes from multiscale stiffness, complex geometry, or solver infrastructure for which nested remeshing is undesirable. A plausible implication is that MOMC targets a different implementation bottleneck from classical MLMC, even though both methods belong to the same variance-reduction tradition.

2. Estimator construction and variance reduction

In the MOMC setting, the goal is to estimate statistics such as $\mathbb{E}[u](x,t)$ for a parametric PDE solution $u(x,t,z)$, where $z$ is the uncertain parameter vector. The standard Monte Carlo estimator is

$E_M[u](x,t) = \frac{1}{M} \sum_{k=1}^M u(x,t,z_k).$

MOMC replaces this single-model average with a hierarchy of order-$l$ solvers evaluated on shared random inputs (Bertaglia et al., 27 Aug 2025).

For a two-level hierarchy with a higher-order solver $L$ and a lower-order solver $L-1$, the construction described in the literature uses $M_L$ expensive samples at order $L$ and $\mathbb{E}[u](x,t)$0 cheaper samples at order $\mathbb{E}[u](x,t)$1, with $\mathbb{E}[u](x,t)$2 matched to the cost ratio. The corresponding estimator is

$\mathbb{E}[u](x,t)$3

The parameter $\mathbb{E}[u](x,t)$4 is a control-variate weight. Its optimal value is

$\mathbb{E}[u](x,t)$5

and, with this choice, the minimized variance becomes

$\mathbb{E}[u](x,t)$6

where $\mathbb{E}[u](x,t)$7 is the Pearson correlation coefficient between adjacent orders (Bertaglia et al., 27 Aug 2025).

The recursive version extends the same idea across multiple orders $\mathbb{E}[u](x,t)$8. Each level uses empirical variances and covariances from shared samples to set quasi-optimal weights, and lower-order models absorb a substantial fraction of the sampling burden because they are cheaper. The method is therefore explicitly a generalized control-variate scheme rather than a pure telescoping-difference estimator in the strict MLMC sense (Bertaglia et al., 27 Aug 2025).

This control-variate viewpoint links MOMC to a wider multifidelity literature. In particular, generalized multifidelity Monte Carlo estimators for higher-order statistics, vector-valued quantities of interest, and nonlinearly statistically dependent models also optimize weights and sample allocations from empirical variance-covariance information (Quaglino et al., 2018). This suggests a conceptual overlap between MOMC and multifidelity control-variate design, even though the MOMC hierarchy is defined by numerical order rather than by a generic collection of low-fidelity surrogates.

3. IMEX Runge-Kutta hierarchies and asymptotic-preserving structure

A defining feature of the multi-order formulation for multiscale hyperbolic systems is its use of Implicit-Explicit Runge-Kutta (IMEX) integrators. IMEX schemes split stiff and non-stiff terms, treating the hyperbolic component explicitly and the stiff relaxation or source term implicitly. In the MOMC framework, they are used at all hierarchy levels so that changes in discretization order do not destroy the asymptotic behavior of the underlying multiscale model (Bertaglia et al., 27 Aug 2025).

The relevant property is asymptotic preserving (AP). In the MOMC literature, the hierarchy is designed so that the AP property is satisfied across different discretization orders of accuracy. An AP-MOMC variant enforces that the AP property of the underlying numerical method is preserved at all hierarchy levels, making the estimator itself asymptotic preserving in the stiff limit (Bertaglia et al., 27 Aug 2025).

This is particularly important for hyperbolic systems with stiff relaxation, kinetic equations, and low Mach number flows, which are identified as principal application classes. In these settings, the same fixed mesh may be retained while the spatial reconstruction and time integrator order are varied, allowing variance reduction without the costly hierarchical re-meshing associated with traditional MLMC (Bertaglia et al., 27 Aug 2025).

The framework also admits a bi-fidelity AP-MOMC variant in which the lowest level is not merely a lower-order discretization of the full model, but a reduced or asymptotic model. In the Jin–Xin example summarized in the literature, the upper levels use full AP-IMEX schemes of different orders, while the lowest level computes the equilibrium or limit solution. The stated purpose is to maximize variance reduction for minimal additional cost (Bertaglia et al., 27 Aug 2025).

4. Error decomposition and computational logic

The MOMC error estimate separates Monte Carlo sampling error from deterministic discretization error. For a deterministic solver of order $\mathbb{E}[u](x,t)$9, the reported bound is

$u(x,t,z)$0

For the multi-order estimator with variance-minimizing weights,

$u(x,t,z)$1

where the first term is the statistical contribution and the second is the deterministic error from the highest-fidelity method (Bertaglia et al., 27 Aug 2025).

The computational logic is therefore the same broad logic that governs MLMC and related methods: use many cheap samples where variance reduction is still effective, and reserve expensive samples for the most accurate approximation. In MLMC, this is formalized through level-dependent sample counts that balance variance $u(x,t,z)$2 and cost $u(x,t,z)$3, often yielding complexity regimes of the form

$u(x,t,z)$4

depending on the relationship between variance-decay and work-growth exponents (Higham, 2015). MOMC does not reuse the mesh-based hierarchy behind those formulas, but it follows the same principle of variance-cost balancing through lower-cost auxiliary models (Bertaglia et al., 27 Aug 2025).

Related methods illustrate alternative ways of altering this balance. Projected MLMC replaces coarse solves with projections of fine-grid solutions onto coarser spaces, reducing cost while retaining the multilevel correction structure (Kim et al., 2015). Unbiased multi-index Monte Carlo removes discretization bias entirely by randomizing the truncation level of the multi-index sum (Crisan et al., 2017). Adaptive MLMC for probabilities refines only ambiguous samples near a discontinuity threshold so as to recover the complexity rates usually associated with smooth functionals (Haji-Ali et al., 2021). These variants clarify that the central design space in hierarchical Monte Carlo concerns not only which hierarchy is used, but also how corrections are coupled, weighted, and terminated.

5. Applications in multiscale hyperbolic systems

The reported applications of MOMC and AP-MOMC are concentrated in multiscale, time-dependent PDEs. The cited examples include prototype nonlinear conservation laws such as inviscid Burgers’ equation, hyperbolic systems with stiff relaxation such as the shallow water equations, and kinetic or multiscale blood-flow models, including one-dimensional blood flow with a viscoelastic vessel wall (Bertaglia et al., 27 Aug 2025).

The deterministic building blocks in these experiments include finite volume IMEX Runge-Kutta schemes of order 1, 2, and 3, spatial reconstructions of different order such as TVD and WENO, and explicit schemes where appropriate. On this basis, the numerical experiments are reported to show that MOMC and AP-MOMC consistently reduce both error and variance of estimated quantities of interest, for the same sample budget or computational cost, when compared with standard Monte Carlo or MLMC (Bertaglia et al., 27 Aug 2025).

Several application-specific findings are emphasized. For stiff relaxation problems, AP-MOMC allows large time steps while avoiding the time-step restriction of explicit schemes. For the shallow water equations, IMEX-based MOMC gives effective variance reduction for quantities such as water depth across randomness in the Froude regime. For the blood-flow model, AP-IMEX hierarchies and bi-fidelity correction deliver variance-reduced estimates for flow and pressure under uncertain vessel properties, with the reduced model becoming especially useful at low error tolerances (Bertaglia et al., 27 Aug 2025).

The 2025 abstract characterizes the overall outcome more generally: the method achieves substantial reduction of both error and variance while maintaining asymptotic consistency in the asymptotic limit (Bertaglia et al., 27 Aug 2025). This formulation is significant because it identifies asymptotic consistency as a coequal design objective with variance reduction, rather than as a secondary numerical property.

6. Distinctions, related developments, and common confusions

A recurring source of confusion is the proximity of MOMC to several neighboring ideas whose names also invoke levels, indices, order, or fidelity. The literature distinguishes these approaches by the object that defines the hierarchy.

MIMC uses a multi-dimensional hierarchy of discretization directions and high-order mixed differences. Under mixed-regularity assumptions, properly chosen total degree index sets can make computational complexity effectively independent of the number of directions up to logarithmic terms (Haji-Ali et al., 2014). Dimension-adaptive MIMC then selects indices greedily by a profit criterion,

$u(x,t,z)$5

and applies this logic, for example, to spatial refinement and multiple Karhunen–Loève truncations in a heat-exchanger model with lognormal random fields (Robbe et al., 2017). MOMC, by contrast, does not organize the approximation space by a multi-index of mesh or truncation levels, but by order of accuracy on a fixed mesh (Bertaglia et al., 27 Aug 2025).

Higher-order multifidelity Monte Carlo addresses a different notion of order. In that setting, “higher-order” refers to target statistics such as variance and sensitivity indices, together with vector-valued quantities of interest and nonlinear statistical dependence between models. The generalized estimators optimize weights and workload allocation analytically or numerically from small batches of high- and low-fidelity runs (Quaglino et al., 2018). MOMC is related through its control-variate interpretation, but its hierarchy is numerical-order based rather than fidelity-generic.

Other recent developments show how the hierarchical Monte Carlo program has diversified. Budgeted MLMC reverses the classical objective by fixing a computational budget on an HPC system and minimizing the mean squared error via dynamic programming following Bellman’s optimality principle (Baumgarten et al., 2023). Nested MLMC with biased and antithetic sampling addresses nested expectations of the form

$u(x,t,z)$6

and obtains order $u(x,t,z)$7 cost under stronger convergence assumptions, or order $u(x,t,z)$8 under more general ones (Haji-Ali et al., 2023). Advanced MLMC methods further embed multilevel structure into MCMC and SMC when exact sampling and exact coupling across levels are unavailable (Jasra et al., 2017, Beskos et al., 2015).

The available literature therefore places Multi-Order Monte Carlo as one specialized branch of a broader hierarchical sampling family: it retains the variance-reduction rationale of MLMC, but shifts the hierarchy from resolution to numerical order, and combines that shift with IMEX/AP structure for multiscale hyperbolic systems (Bertaglia et al., 27 Aug 2025).