Multilevel Monte Carlo for Polytopic DG Discretisations with Spatially Represented Random Coefficients

1. Introduction

• Random-coefficient elliptic PDE statistics are expensive because realistic problems require repeated forward solves, and standard Monte Carlo becomes costly at fine spatial resolution.MLMC reduces work by using many inexpensive coarse solves and fewer expensive fine solves.Source: §1 · p. 1
• The DG-norm Strang decomposition separates forward spatial error, projected-flux consistency, and random-coefficient representation error without requiring coefficient traces.This distinction makes approximation of the data driving the forward problem explicit in the total bias.Source: §1 · p. 2
• The spectral FFT sampler drives all levels with one represented latent field on a fixed Cartesian grid, allowing unrelated forward meshes; the SPDE Matérn sampler couples agglomerated levels through restricted white-noise loads.The SPDE hierarchy is naturally nested, with the coarse load obtained by exact restriction of the fine-level load.Source: §1 · pp. 2–3
• The unified analysis relates coefficient-error bounds to latent-field representation and interpolation for the spectral sampler, and DG-resolvent approximation of white noise for the SPDE sampler.Both strategies use the same forward DG technology.Source: §1 · p. 2
• Existing work covers white-noise coupling for finite-element SPDE discretisations, SPDE-generated coefficients in multilevel quasi-Monte Carlo, and polytopal Virtual Element MLMC where coefficients are continuum realizations.This paper instead includes spatial coefficient representation in the multilevel discretisation.Source: §1 · pp. 1–2
• Experiments on Voronoi and agglomerated polytopic meshes confirm expected forward and coefficient-representation errors, strong multilevel coupling, and decay of coupled corrections.End-to-end experiments compare complete MLMC and single-level Monte Carlo work–accuracy behavior.Source: §1 · p. 3

2. Problem setup

• The setup uses the usual Lebesgue and Sobolev spaces, a complete probability space, and Bochner spaces for random variables valued in function spaces.Expectation and variance are taken in the usual sense.Source: §2 · p. 3

2.1. Elliptic model with random coefficients.

• Under uniform ellipticity and bounded reaction assumptions, the variational problem has a unique solution almost surely with a uniform H1 stability bound.The bound depends on the coefficient bounds and domain, and scales with the L2 norm of the forcing.Source: §2.1 · p. 3

2.3. Random input models.

• The random input models allow truncated lognormal diffusion and reaction coefficients generated from mean-zero Gaussian fields, which may be correlated.A single latent field may drive both coefficients; distinct or correlated fields are also admissible, provided the uniform coefficient bounds hold.Source: §2.3 · p. 4

2.5. Quantities of interest.

• Quantities of interest are functionals of the solution that are Lipschitz continuous in H1, including bounded linear functionals and regional averages.Regional averages provide a spatially localized alternative because point evaluation is not continuous on H1 when d ≥ 2.Source: §2.5 · pp. 4–5

3. Polytopic DG discretisation and Strang decomposition

• The forward DG problem uses represented coefficient fields that satisfy uniform ellipticity and positivity bounds across levels.The coefficient pair is the one actually seen by the discretisation.Source: §3 · p. 5

3.1. Polytopic meshes and broken spaces.

• The face penalty is scaled by neighboring element diameters, so stabilization does not degenerate when a polytope has small-measure faces.The elements are still subject to a uniform chunkiness condition.Source: §3.1 · p. 5

3.3. Lifting operator and the DG form. ¶

• The numerical normal flux is the L2 projection of the physical flux, avoiding coefficient traces in the method's definition.The DG form uses a lifting and symmetric interior-penalty construction.Source: §3.3 · p. 7
• For sufficiently large penalty parameter, the DG form is coercive and continuous, with constants independent of level and face count or relative face sizes.The discrete problem is uniquely solvable under the same hypotheses.Source: §3.3 · pp. 7–8

3.7. Strang decomposition.

• The abstract Strang lemma bounds the DG error by consistency terms under the coercivity and coefficient-bound hypotheses.The result applies to the discrete problem with the stated uniform coefficient bounds.Source: §3.7 · p. 8
• The decomposition distinguishes inconsistency from the lifted, projected numerical flux from error due to coefficient representation.The latter arises when (A, c) is replaced by (eAℓ, ecℓ).Source: §3.7 · p. 9
• The resulting theorem gives a Strang decomposition under the preceding assumptions.The supplied passage identifies the decomposition and its assumptions but does not provide the bound's explicit terms.Source: §3.7 · p. 9

3.10. Approximation and lifting consistency.

• Polynomial approximation on uniformly star-shaped polytopes controls solution approximation, while projected-flux consistency is controlled by best approximation of the physical flux.The regularity assumptions are stated directly in W 1+s,r.Source: §3.10 · p. 9

3.12. Coefficient representation defect.

• Lipschitz coefficient maps propagate latent-field approximation error to the represented diffusion and reaction coefficients.The bounds hold on the range relevant to the truncated latent field, whose truncations preserve uniform coefficient bounds.Source: §3.12 · p. 10, §4 · p. 11
• The error decomposition groups solution approximation with projected-flux consistency and isolates the remaining terms as coefficient-representation error.This is the decomposition used in the subsequent multilevel analysis.Source: §3.12 · p. 11
• The spectral construction shares one represented latent-field realization across forward meshes, whereas the SPDE construction couples adjacent levels through consistent restrictions of the same white noise.The two constructions are alternative couplings, not pathwise coupled to one another; their bounded-domain and periodic models differ near the boundary.Source: §4 · p. 11

4.1. Spectral sampling on a Cartesian reference grid.

• A Fourier representation of a mean-zero periodic Gaussian field is sampled on a Cartesian reference grid and interpolated for use on the domain.Its parameters determine correlation length, regularity, and variance; the computational sample is the Fourier truncation of the same field.Source: §4.1 · p. 12
• The same realized reference field drives every forward level, so both solves in a coupled MLMC correction use the same represented coefficient realization without transferring coefficient fields between meshes.The reference-grid resolution is selected for the target hierarchy's finest scale.Source: §4.1 · pp. 12–13
• Under Lipschitz coefficient maps, latent-field approximation error propagates to coefficient error, and the reference-grid approximation result can be combined with that bound.The supplied lemmas state these error bounds but do not provide their explicit rates here.Source: §4.1 · p. 13

4.4. SPDE sampling on an agglomerated hierarchy.

• The SPDE sampler solves a reaction–diffusion equation with spatial white-noise forcing using the lifting-based symmetric interior-penalty discretisation on an agglomerated hierarchy.Its homogeneous Dirichlet condition is imposed weakly and pins the latent field to zero at the boundary, modifying covariance there.Source: §4.4 · pp. 13–14

4.5. Discrete white noise and inter-level coupling.

• The discrete white-noise load vector has a centered Gaussian law with covariance given by the DG mass matrix.Because the DG mass matrix is block diagonal by elements, the load can be sampled locally with independent elemental loads.Source: §4.5 · p. 14
• Adjacent levels use restrictions of the same continuum white-noise functional, preserving the required coarse marginal while coupling the inputs.The discrete load vector represents white noise acting on the DG space, not the coefficient vector of an L2 function.Source: §4.5 · p. 14

4.6. Approximation of the SPDE latent field.

• The SPDE latent-field approximation has rates αY = 1 in d = 2 and αY = 1/2 in d = 3.The strong estimate is obtained under the stated mesh and H2 regularity assumptions.Source: §4.6 · p. 15
• The analysis propagates latent-field approximation to the represented diffusion and reaction coefficients, whose errors are uniformly bounded by deterministic coefficient bounds.For the SPDE coefficient estimate, q must be finite; one admissible two-dimensional choice gives αspde = 1/2.Source: §4.6 · p. 15

4.9. Transfer to the forward problem.

• Theorem 4.10 transfers coefficient approximation error to forward error in the DG norm, under a separate regularity assumption on the random forward problem.The stated transfer uses coefficient-error norms and forward-solution regularity.Source: §4.9 · p. 16
• For FFT sampling, the same represented reference field is used at every level, so its approximation error contributes to bias but not to a level-dependent correction defect.The coupled correction measures spatial approximation of the represented continuum output.Source: §4.9 · p. 16
• The resulting strong single-level estimates provide the bias and correction bounds used in the multilevel analysis.The coefficient-sampling estimates are combined with the lifting-based Strang bound.Source: §4.9 · pp. 15, 17

5. Multilevel estimator and total error decomposition

• The two samplers approximate coefficients differently: FFT uses one represented field throughout a hierarchy, while SPDE uses level-dependent latent fields coupled by consistent white-noise restrictions.The levels are combined using the standard MLMC bias–variance decomposition.Source: §5 · p. 17

5.1. Multilevel estimator.

• Each level correction is formed from coupled outputs: a shared reference-field realisation for FFT or a fine white-noise load and its restriction for SPDE.The estimator uses δQ0 = Q0 and δQℓ = Qℓ − Qℓ−1 for subsequent levels.Source: §5.1 · p. 17

5.2. Bias and total error.

• For FFT, the representation error controls bias relative to the continuum quantity, while the coupled correction measures spatial approximation of the represented output.Both levels approximate the same represented continuum output.Source: §5.2 · p. 18

5.6. Correction variance and computational complexity.

• The MLMC complexity result gives root-mean-square error of order ε for sampling variance of order ε², with expected work governed by the stated rate assumptions.The work bound is conditional on the theorem's assumptions on bias, correction variance, and coupled-sample cost.Source: §5.6 · pp. 18–19

6.1. Coupled coefficient fields and solutions.

• FFT coefficients are pathwise common across the hierarchy, whereas SPDE coefficient fields are distinct level approximations driven by the same underlying discrete white noise.Figure 1 shows coefficients and P1 DG solutions for 64, 256, and 1024-cell meshes.Source: §6.1 · p. 20, Figure 1 · p. 20

6.3. Coefficient representation.

• FFT latent-field and diffusion errors decay approximately first order under reference-spacing refinement, with fitted exponents close to 1.06.With a fixed under-resolved FFT field, forward-mesh refinement leaves a representation-error floor; refining the field with the mesh removes it.Source: §6.3 · p. 21, Figure 3 · p. 21

6.4. Multilevel correction decay.

• For SPDE coupling, the fitted weak, variance, and cost exponents are 1.39, 2.93, and 1.28, respectively; variance decay exceeds cost growth for both strategies.The reported exponents are bαspde = 1.39, bβspde = 2.93, and bγspde = 1.28.Source: §6.4 · p. 22
• Neighbouring-level correlations approach one on the finest pair, while the fine-output variances are 1.30 × 10^-5 for FFT and 1.42 × 10^-5 for SPDE.The correction-to-fine variance ratios and 1 − ρ decay across levels.Source: §6.4 · p. 22, Figure 4 · p. 22

6.5. End-to-end MLMC work.

• MLMC reduces measured work relative to finest-level Monte Carlo for both samplers across the tested tolerance range, with gains increasing at tighter tolerances.Figure 5 compares work–accuracy curves; Table 4 reports realised RMSE and measured work.Source: §6.5 · pp. 22–23, Figure 5 · p. 23, Table 4 · p. 23

6.2. Forward discretisation verification.

• The experiments confirm expected P1 polytopic DG rates and test spatial convergence on independently generated Voronoi meshes.The broader experiments use general Voronoi and agglomerated polytopic meshes.Source: §6 · p. 19, §6.2 · p. 21, §6.5 · p. 23

• The framework separates spatial discretisation error from error introduced by representing random coefficients in a lifting-based symmetric interior penalty DG method.Source: §7 · p. 23
• It uses a common-reference FFT construction and an SPDE-based Matérn construction coupled through projected white-noise loads on agglomerated hierarchies.Source: §7 · p. 23
• Experiments on Voronoi and agglomerated polytopic meshes confirm expected forward rates, illustrate coefficient-representation error, and show strong decay of coupled multilevel corrections.Source: §7 · pp. 23–24