Multilevel Monte Carlo Methods for Stochastic Convection–Diffusion Eigenvalue Problems

We develop new multilevel Monte Carlo (MLMC) methods to estimate the expectation of the smallest eigenvalue of a stochastic convection–diffusion operator with random coefficients. The MLMC method is based on a sequence of finite element (FE) discretizations of the eigenvalue problem on a hierarchy of increasingly finer meshes. For the discretized, algebraic eigenproblems we use both the Rayleigh quotient (RQ) iteration and implicitly restarted Arnoldi (IRA), providing an analysis of the cost in each case. By studying the variance on each level and adapting classical FE error bounds to the stochastic setting, we are able to bound the total error of our MLMC estimator and provide a complexity analysis. As expected, the complexity bound for our MLMC estimator is superior to plain Monte Carlo. To improve the efficiency of the MLMC further, we exploit the hierarchy of meshes and use coarser approximations as starting values for the eigensolvers on finer ones. To improve the stability of the MLMC method for convection-dominated problems, we employ two additional strategies. First, we consider the streamline upwind Petrov–Galerkin formulation of the discrete eigenvalue problem, which allows us to start the MLMC method on coarser meshes than is possible with standard FEs. Second, we apply a homotopy method to add stability to the eigensolver for each sample. Finally, we present a multilevel quasi-Monte Carlo method that replaces Monte Carlo with a quasi-Monte Carlo (QMC) rule on each level. Due to the faster convergence of QMC, this improves the overall complexity. We provide detailed numerical results comparing our different strategies to demonstrate the practical feasibility of the MLMC method in different use cases. The results support our complexity analysis and further demonstrate the superiority over plain Monte Carlo in all cases.

The conductivity κ(x, ω) : D × Ω → R is a log-uniform random field (as used in, e.g., [19]), defined using the process convolution approach in [37], such that log κ(x, ω) = Z(x, ω) with k(x − c i ) a kernel centered at a certain number of points c i ∈ D and i.i.d.uniform random variables ω i ∼ U [0, 1].Similarly, the convection velocity a(x, ω) : D × Ω → R d can also be some bounded random variable, which also depends on uniform random variables ω i ∼ U [0, 1] and is additionally assumed to be divergence-free, i.e., The purpose of this paper is to compute the expectation of the smallest eigenvalue of (1), using multilevel Monte Carlo methods.Stochastic eigenvalue problems arise in a variety of physical and scientific applications and their numerical simulations.Factors such as measurement noise, limitations of mathematical models, the existence of hidden variables, the randomness of input parameters, and other factors contribute to uncertainties in the modelling and prediction of many phenomena.Applications of uncertainty quantification (UQ) specifically related to eigenvalue problems include: nuclear reactor criticality calculations [2,3,25], the derivation of the natural frequencies of an aircraft or a naval vessel [41], band gap calculations in photonic crystals [22,27,55], the computation of ultrasonic resonance frequencies to detect the presence of gas hydrates [51], the analysis of the elastic properties of crystals with the use of rapid measurements [52,61], or the calculation of acoustic vibrations [12,66].Stochastic convection-diffusion equations are used to describe simple cases of turbulent [24,44,54,63] or subsurface flows [64,67].
Monte Carlo sampling is one of the most popular methods for quantifying uncertainties in quantities of interest coming from stochastic PDEs.Although simple and robust, Monte Carlo methods can be severely inefficient when applied to UQ problems, because their slow convergence rate often requires a large number of samples to meet the desired accuracy.To improve the efficiency, the multilevel Monte Carlo (MLMC) method was developed, where the key idea is to reduce the computational cost by spreading the samples over a hierarchy of discretizations.The main idea was introduced by Heinrich in 2001 [36] for path integration, then generalized by Giles in 2008 [30] for SDEs.More recently, MLMC methods have been applied with great success to stochastic PDEs, see, e.g., [6,7,14,53,60,65] and [28,29] specifically for eigenproblems.A general overview of MLMC is presented by Giles in [31].
In this paper, we present a MLMC method to approximate (4), which, motivated by the use of MLMC for source problems described above, is based on a hierarchy of discretizations of the eigenvalue problem (1) and which is much more efficient in practice than a Monte carlo approximation.We consider two discretization methods, a standard Galerkin finite element method (FEM) and a streamline upwind Petrov-Galerkin (SUPG) method.The SUPG method improves the stability of the approximation for cases with high convection and also allows us to start the MLMC method from a coarser discretization.To further reduce the cost of our MLMC method, we again exploit the hierarchy of discretizations by using approximations on coarse levels as the starting values for the eigensolver on the fine level.We also present the two extensions of MLMC that aim to improve different aspects of the method.First, to improve the stability of the eigensolver for each sample we include a homotopy method for solving convection-diffusion eigenvalue problems in the MLMC algorithm.The homotopy method computes the eigenvalue of the convection-diffusion operator by following a continuous path starting from the pure diffusion operator.Second, to improve the overall complexity we present a multilevel quasi-Monte Carlo method that aims to speed up the convergence of the variance on each level by replacing the Monte Carlo samples with a quasi-Monte Carlo (QMC) quadrature rule.
The structure of the paper is as follows.Section 2 introduces the variational formulation of (1), along with necessary background material on stochastic convection-diffusion eigenvalue problems.Two discrete formulations of the eigenvalue problem are introduced: the Galerkin FEM and the SUPG method.Section 3 introduces the MLMC method and presents the corresponding complexity analysis.In particular, this section details how to efficiently use each eigensolver, the Rayleigh quotient and implicitly restarted Arnoldi iterations, within the MLMC algorithm.In Section 4, we present the two extensions of our MLMC algorithm: a homotopy MLMC and a multilevel quasi-Monte Carlo method.Section 5 presents numerical results for finding the smallest eigenvalue of the convectiondiffusion operator in a variety of settings.In particular, we present examples for difficult cases with high convection.
To ease notation, for the remainder of the paper we combine the random variables in the convection and diffusion coefficients into a single uniform random vector of dimension s < ∞, denoted by ω = (ω i ) s i=1 with ω i ∼ U [0, 1].In this case, π is the product uniform measure on Ω := [0, 1] s .

Variational formulation
The eigenvalue problem (1) needs to be discretized, because its solution is not analytically tractable for arbitrary geometries and parameters.As such, we apply the standard finite element method to (1) to obtain an approximation of the desired eigenpair (λ, u).
Before deriving the variational form of (1), we first establish certain assumptions about the problem domain, the random field κ(ω) and the velocity field a(ω) for ω ∈ Ω, which, in particular, ensure that the solution is in H 2 (D) [33] as well as incompressibility.
Assumption 1. Assume that D ⊂ R d , for d = 1, 2, or 3, is a bounded, convex domain with Lipschitz continuous boundary Γ.
A simple example of a random convection term is a homogeneous convection, a(x, ω) = [a 1 ω 1 , . . ., a d ω d ] ⊤ for a 1 , . . ., a d ∈ R, which are independent of x.Another example is the curl of random vector field, e.g., a(x, ω) = ∇ × Z(x, ω) where Z is a vector-valued random field similar to that defined in (2).Both of these examples satisfy Assumption 3.
Next we introduce the variational form of (1).Whenever it does not lead to confusion, we drop the spatial coordinate of (stochastic) functions for brevity-for example, u(x, ω) is also written as u(ω).Let V = H 1 0 (Ω) be the first-order Sobolev space of complex-valued functions with vanishing trace on the boundary with norm v V = ∇v L 2 .Then let V * denote the dual space of V .Multiplying (1) by a test function v ∈ V and then performing integration by parts, noting that we have no Neumann boundary condition term since u(x, ω) = 0 on Γ, we obtain The variational eigenvalue problem corresponding to (1) is then: Find a non-trivial eigenpair (λ(ω), where and Since the velocity a is divergence free, ∇•a = 0, the sesquilinear form in (5) is uniformly coercive, i.e., with a min > 0 independent of ω.It is also uniformly bounded, i.e., with a max < ∞ independent of ω.
In addition to the primal form (5), to facilitate our analysis later on we also consider the dual eigenproblem: Find a non-trivial dual eigenpair (λ The primal and dual eigenvalues are related to each other via λ(ω) = λ * (ω).
Finally, using classical results in Grisvard [33] it follows that where C D depends only on the domain D. Finally, substituting in the bound on f ω L 2 (11) gives the desired upper bound (10).
The result for the dual eigenfunction follows analogously.

Finite element formulation
Let {T h } h>0 be a family of (quasi-)uniform, shape-regular, conforming meshes on the spatial domain D, where each T h is parameterised by its mesh width h > 0. For h > 0, we approximate the infinite-dimensional space V by a finite-dimensional subspace V h .In this paper, we consider piecewise linear finite element (FE) spaces, but the method will work also for more general spaces.The resulting discrete variational problem is to find non-trivial primal and dual eigenpairs (λ(ω), and For each ω, it is well-known that for h sufficiently small the FE eigenvalue problem (12) admits M h := dim(V h ) eigenpairs, denoted by which approximate the first M h eigenpairs of (5).This approach is also called the Galerkin method.
In convection-dominated regions, the Galerkin method has well-known stability issues for standard (Lagrange-type) FEs, if the element size h does not capture all necessary information about the flow.The Peclet number (sometimes called the mesh Peclet number) [68] Pe governs how small the mesh size h should be in order to have a stable solution using basic (Lagrange-type) FE methods.
The error in the FE approximations ( 14) can be analysed using the Babuška-Osborn theory [4].We state the error bounds for a simple eigenpair.
Theorem 2. Let (λ(ω), u(ω)) be an eigenpair of (5) that is simple for all ω ∈ Ω, where Ω is a compact domain.Then there exist constants C λ , C u , independent of h and ω, such that and u h (ω) can be normalized such that Proof.See Babuška and Osborn [4] and the appendix, where we show explicitly that the constants are bounded uniformly in ω.

Streamline-upwind Petrov-Galerkin formulation
A sufficiently small Peclet number (15) guarantees numerical stability of the standard Galerkin method.One can either choose a small overall mesh size h or locally adapt the mesh size to satisfy the stability condition.However, globally reducing the mesh size may lead to a high computational cost, while local adaptations may need to be performed path-wise for each realisation of ω, which in turn leads to complications in the algorithmic design.In this section, we consider using the streamline-upwind Petrov-Galerkin (SUPG) method to improve numerical stability.
The SUPG method was introduced by Brooks and Hughes [10] to stabilize the finite element solution.Since then, the method has been extensively investigated and used in various applications [8,15,35,40,39,43].The SUPG method can be derived in several ways.Here, we introduce its formulation by adding a stabilization term to the bilinear form.An equivalent weak formulation can be obtained by defining a test space with additional test functions in the form v(x) = v(x) + p(x), where v(x) is a standard test function in the finite element method and p(x) is an additional discontinuous function.
We define the residual operator R as which gives the residual of the convection-diffusion equation (1) for a pair (σ, v) ∈ C × V .Then, stabilization techniques can be derived from the general formulation where |T h | is the number of elements of the mesh T h , P(ω) is some stabilization operator and τ m (ω) is the stabilization parameter acting in the mth finite element.The stabilization strategy will be determined by P(ω) and τ m (ω).
Various definitions exist for the operator P(v, ω), such as the Galerkin Least Square method [38], the SUPG method [9,10,23], the Unusual Stabilized Finite Element method [5], etc.For the SUPG method, the stablization operator P(ω) is defined as Substituting Equations ( 18) and ( 20) into (19) gives the SUPG weighted residual formulation which is equivalent to After approximating the weak form (21) by the usual finite-dimensional subspaces, we obtain the discrete variational problem: Find non-trivial (primal) eigenpairs (λ h (ω), It follows that the right-hand side matrix is no longer symmetric and is stochastic compared to the mass matrix in the standard Galerkin method.
In general, finding the optimal stabilization parameter τ m (x, ω) is an open problem, and thus it is defined heuristically [43].We employ the following stabilization parameter [8,35] However, in practical implementations the following asymptotic expressions of τ m (x, ω) are used Figure 1 shows the 20 smallest eigenvalues for a single realization of random field κ(x, ω) with velocity a(x, ω) = [50, 0] T on meshes with size h = 2 −3 , 2 −4 , 2 −5 .The standard Galerkin method has non-physical oscillations in the discretized eigenfunction for such a coarse mesh and its two smallest eigenvalues form a complex conjugate pair; this contradicts the fact that the smallest eigenvalue should be real and simple.The SUPG method, on the other hand, has a real smallest eigenvalue, indicating a stable solution.

Multilevel Monte Carlo methods
To compute E[λ], we first approximate the eigenproblem (5) for each ω ∈ Ω and then use a sampling method to estimate the expected value of the approximate eigenvalue.There are two layers of approximation: First the eigenvalue problem is discretized by a numerical method, e.g., FEM or SUPG as in Section 2.1, then the resulting discrete eigenproblem is solved by an iterative eigenvalue solver, e.g., the Rayleigh quotient method, such that 0 500 1,000 1,500 2,000 2,500 3,000 (a) SUPG eigenvalues.0 500 1,000 1,500 2,000 2,500 3,000 , where h denotes the meshwidth of the spatial discretization and K denotes the number of iterations used by the eigenvalue solver.
Applying the Monte Carlo method to λ h,K , the expected eigenvalue can be approximated by the estimator where the samples {ω n } N n=1 ⊂ Ω are i.i.d.uniformly on Ω.This introduces a third factor that influences the accuracy of the estimator in (26) in addition to h and K, namely the number of samples N .Note that we assume that the number of iterations K is uniformly bounded in ω.
The standard Monte Carlo estimator in ( 26) is computationally expensive.To measure its accuracy we use the mean squared error (MSE) where the outer expectation is with respect to the samples in the estimator Y h,K,N .Under mild conditions, the MSE can be decomposed as In this decomposition, the bias | is controlled by h and K, whereas the variance term decreases linearly with 1/N .To guarantee that the MSE remains below a threshold ε 2 , h and K need to be chosen such that the bias is O(ε 2 ), while the sample size needs to satisfy N = O(ε −2 ).Suppose K = K(h) is sufficiently large so that the bias is solely controlled by h and satisfies for some α > 0. Suppose further that the computational cost to compute λ h,K (ω) for each ω is O(h −γ ) for some γ > 0. Then the total computational complexity to achieve an MSE of ε 2 is O(ε −2−γ/α ).Note that in the best-case scenario, we have γ = d, i.e., when the computational cost of an eigensolver iteration is linear in the degrees of freedom of the discretization and the number of iterations can be bounded independently of h.Due to the quadratic convergence of algebraic eigensolvers, K is usually controlled very easily.The multilevel Monte Carlo (MLMC) method offers a natural way to reduce the complexity of the standard Monte Carlo method by spreading the samples over a hierarchy of discretizations.In our setting, we define a sequence of meshes corresponding to mesh sizes This in turn defines a sequence of discretized eigenvalues λ h 0 ,K 0 (ω), λ h 1 ,K 1 (ω), . . ., λ h L ,K L (ω) that approximate λ(ω) with increasing accuracy and increasing computational cost.The MLMC method approximates E[λ(ω)] using the telescoping sum where λ ℓ (ω) := λ h ℓ ,K ℓ (ω) is the shorthand notation for the discretized eigenvalues.Each expected value of differences in ( 27) can be estimated by an independent Monte Carlo approximation, leading to the multilevel estimator Suppose independent samples are used to compute each Y ℓ , then and the MSE of ( 28) can also be split into a bias and a variance term, i.e., and the variance, var[Y ], are both less than 1 2 ε 2 .The following theorem from [14] (see also [31]) provides bounds on the computational cost of a general MLMC estimator and applies in particular to (28).
Theorem 3. Let Q denote a random variable and Q ℓ its numerical approximation on level ℓ, and suppose C ℓ is the computational cost of evaluating one realization of the difference where then for any 0 < ε < e −1 there exist a constant c, a stopping level L, and sample sizes where the constant c is independent of α, β and γ.
For a given ε, from [14] the maximum level L in Theorem 3 is given by where c I is the implicit constant from Assumption I (convergence of bias) above.The optimal sample sizes, {N ℓ }, that minimize the computational cost of the multilevel estimator in Theorem 3 are obtained using a standard Lagrange multipliers argument as in [14] and are given by Since β > 0, Theorem 3 shows that for all cases in (31), the MLMC complexity is superior to that of Monte Carlo.When β > γ, the variance reduction rate is larger than the rate of increase of the computational cost, and thus most of the work is spent on the coarsest level.In this case, the multilevel estimator has the best computational complexity.When β < γ the total computational work of the multilevel estimator may only have a marginal improvement compared to that of the classic Monte Carlo method.
Corollary 1 (Order of convergence).For ω ∈ Ω, let h > 0 be sufficiently small and consider two finite element approximations, cf. ( 12), of the smallest eigenvalue λ(ω) of the eigenvalue problem (5) with h ℓ−1 = h and h ℓ = h/2.The expectation of their difference is bounded by while the variance of the difference is bounded by for two constants c 1 , c 2 that are independent of ω, h and ℓ.
Proof.Applying Theorem 2, since C λ is independent of ω we have and Therefore, by the triangle inequality, we have The variance reduction rate comes from the following relation and, similarly, by the Cauchy-Schwarz inequality Remark 1.In our numerical experiments, we observed that the SUPG approximation of the eigenvalue problem, cf. ( 22), has similar rates of convergence α and β in MLMC compared to the standard finite element approximation.
An important physical property of the smallest eigenvalue of ( 5) is that it is real and strictly positive.Clearly, E[λ] > 0 as well, and so we would like our multilevel approximation (28) to preserve this property.Below we show that a multilevel approximation based on Galerkin FEM with a geometrically-decreasing sequence of meshwidths is strictly positive provided that h 0 is sufficiently small.Proposition 2. Suppose that h ℓ = h 0 2 −ℓ for ℓ ∈ N with h 0 > 0 sufficiently small and let λ h ℓ (•) be the approximation of the smallest eigenvalue using the Galerkin FEM as in (12).Then, for any L ∈ N, the multilevel approximation of the smallest eigenvalue is strictly positive, i.e., Proof.First, since λ is continuous and strictly positive on Ω it can be bounded uniformly from below, i.e., there exists q λ > 0 such that For ℓ = 0, using ( 16) and ( 40) we can bound λ h 0 (ω) uniformly from below by Since this bound is independent of ω, it follows that Similarly, for ℓ ≥ 1 using ( 16) we obtain Again, this bound is independent of ω and so Finally, we bound the multilevel approximation Y from below using ( 41) and ( 42) as follows, where we have used the property that h 0 is sufficiently small, i.e., h 0 ≤ q λ/(12C λ ), to ensure Y > 0, as required.
The result above can be extended beyond the geometric sequence of FE meshwidths to a general sequence of FE meshwidths, provided that L ℓ=0 h 2 ℓ is sufficiently small.Similarly, as in Remark 1, we observe that the MLMC approximations based on SUPG are also strictly positive.
Choosing the number of iterations K ℓ such that the error of the eigensolver is of the same order as the FE error on each level, i.e., |λ h ℓ (ω) − λ h ℓ ,K ℓ (ω)| h 2 ℓ for all ℓ = 0, 1, . . ., L and ω ∈ Ω, it can similarly be shown that the multilevel approximation (28) also satisfies Y > 0.
To obtain the eigenvalue approximation on level ℓ, choosing a basis for the FE space V ℓ := V h ℓ in (12) leads to a generalized (algebraic) eigenproblem in matrix form for each sample ω, i.e., where u ℓ (ω) is the coefficient vector (with respect to the basis) and A ℓ (ω), M ℓ (ω) are the associated FE matrices corresponding to the mesh T ℓ := T h ℓ .The number of iterations K in the computational cost per sample, as well as the rate of the cost per iteration depend on the choice of the algebraic eigensolver.A variety of solvers can be applied here to solve the generalized eigenvalue problem (43), including power iteration, the QR algorithm, subspace iterations, etc.For our purposes, we only need an eigensolver that is able to compute the smallest eigenvalue, which is real and simple.As such, we consider here two eigenvalue solvers, the Rayleigh quotient iteration and the implicitly restarted Arnoldi method.

9:
i ← i + 1 10: end while 11: Output: (η, ξ, λ) We first consider the Rayleigh quotient iteration (Alg.1), introduced first by Lord Rayleigh in 1894 for a quadratic eigenproblem of oscillations of a mechanical system [57] and then extended in the 1950s and 60s to non-symmetric generalized eigenproblems [17,56].The following lemma, whose proof can be found in Crandall [17] and Ostrowski [56], establishes the error reduction rate of the Rayleigh quotient iteration, which will in turn help to bound the computational cost on each level.Lemma 1. Suppose we have an initial guess λ ℓ,0 (ω) to the eigenvalue λ ℓ (ω) at the level ℓ and |λ ℓ,0 (ω) − λ ℓ (ω)| is sufficiently small.Then the sequence λ ℓ,i (ω) converges to λ ℓ (ω) quadratically, i.e., there exists a constant Ĉ(ω) such that The computational cost of Rayleigh quotient iteration (RQI) is dominated by the cost of solving two linear systems in each iteration (cf.Lines 6 and 7 of Alg. 1).For direct solvers, such as LU decomposition, the computational cost depends on the sparsity and bandwidth of the matrices, e.g., for piecewise linear FE applied to (5) and d = 2, the cost for solving these linear systems on level ℓ is O(h −3 ℓ ) [26].However, optimal iterative solvers, such as geometric multigrid methods, are able to achieve the optimal computational complexity of (or close to) O(h −d ℓ ).All other steps in Alg. 1 are linear in the degree of freedoms, and thus O(h −d ℓ ).Hence, typically the cost per iteration grows with rate γ ≥ d, but it can be as big as γ = 3 for d = 2.The remaining factor in the computational cost is the number of iterations K for the Rayleigh quotient iteration within the MLMC estimator, but this is independent of h ℓ .
Output: λ ℓ − λ ℓ−1 12: end if Recall the MLMC estimator (28), where at each level ℓ we compute the differences λ ℓ (ω n ) − λ ℓ−1 (ω n ) for the same sample ω n .The number of RQI iterations needed for a sufficiently accurate approximation of λ ℓ (ω n ) -the more costly level ℓ computation -can be significantly reduced by using the computed approximation of the eigenvalue λ ℓ−1 (ω n ) on the coarser level as the initial guess, thus also reducing the total computational cost.In fact, we design a three-grid method, similar to the one used in [29] to implement this strategy, which uses the approximate eigenvalue λ 0 (ω n ) on level zero with mesh size h 0 as the initial guess for computing eigenvalue λ ℓ−1 (ω n ) on level ℓ − 1.Then, λ ℓ−1 (ω n ) is used as the initial guess for computing λ ℓ (ω n ); see Alg. 2 for details.
To estimate the computational cost of this three-grid method, we choose again h ℓ−1 = h = 2h ℓ and denote the exact discrete eigenvalues on level ℓ − 1 and level ℓ by λ h (ω n ) and λ h/2 (ω n ), respectively.The goal is to control the errors of the eigenvalues λ ℓ−1 (ω n ) and λ ℓ (ω n ) actually computed using Alg. 2 to be within the respective discretization errors.Due to the quadratic convergence rate of the RQI (cf.Lemma 1), often only two or three iterations are sufficient to compute a sufficiently accurate approximation λ 0 (ω n ) on Level 0 in Line 3 of Alg. 2. Similarly, in Line 5 of Alg. 2, two to three iterations of RQI are again sufficient to ensure that the error of the estimated eigenvalue λ ℓ−1 (ω n ) satisfies which is the bound on the discretization error on level ℓ − 1 in Theorem 2. When λ ℓ−1 (ω n ) is then used as the initial guess for estimating λ h/2 (ω n ), the initial error satisfies using triangle inequality and Theorem 2 again.Therefore, using Lemma 1 for sufficiently small mesh size h such that h ≤ 2 9 Ĉ(ω n )C λ −1/2 , one single iteration of RQI on level ℓ In practice, two iterations of RQI are typically used to achieve the target accuracy for λ ℓ (ω n ) in Line 10 of Alg. 2. These two calls to RQI dominate the computational cost of Alg. 2 with their four linear solves.Hence, for sparse direct solvers and d = 2, the overall computational cost of Alg. 2 is O(h −3 ℓ ) and γ = 3 in Theorem 3. The computational complexity of Alg. 2 can be further reduced using multigrid-based methods to efficiently solve the Rayleigh quotient iterations [11] that potentially offer a rate of γ = d (or close to) even in three dimensions.However, it is unclear if the same rate of convergence as for self-adjoint operators can be retained for the convection-dominated problems we are considering here.end for 10: end for We also consider the implicitly restarted Arnoldi method [1,48,58,59,62] and its implementation in the library ARPACK [49] to solve the eigenvalue problem.Compared to the Rayleigh quotient iteration, the Arnoldi method calculates a specified number of eigenpairs that depend on the dimension of the Krylov subspace.The performance of the implicitly restarted Arnoldi method is determined by several factors such as the dimension of the Krylov subspace and the initial vector.To the best of the authors' knowledge, for the eigenvalue problem (12) we are considering here, the convergence rate, and therefore the computational cost, of the implicitly restarted Arnoldi method is not yet known.As such, we numerically estimate the rate variable γ and the computational cost C ℓ for determining the optimal sample sizes in MLMC.It appears that the number of iterations grows slightly faster than O(h −1 ℓ ) leading to a similar total complexity as RQI for d = 2 of γ ≈ 3.5.

Extensions of MLMC method
In this section, we introduce two extensions of the MLMC method for convection-diffusion eigenvalue problems.First, we employ a homotopy method to add stability to the eigensolve for each sample.Second, we replace the Monte Carlo approximation of the expected value on each level in (27) with a quasi-Monte Carlo (QMC) method, which, due to the faster convergence of QMC, allows us to use less samples on each level and improves the overall complexity.

Homotopy multilevel Monte Carlo method
In Carstensen et al. [13] a homotopy method is employed to solve convection-diffusion eigenvalue problems with deterministic coefficients, using the homotopy method to derive adaptation strategies for FE methods.The authors also provided estimates on the convergence rate of the smallest eigenvalue with respect to the homotopy parameter.We aim to investigate the application of this homotopy method in the MLMC method, particularly in designing multilevel models for alleviating numerical instability (due to the high advection velocity) on coarser meshes.For eigenvalue problems, the homotopy method [50] uses an initial operator L 0 -for which the target eigenvalue is easier to compute than that of the original operator L-to form a continuation with a function f : [0; 1] → [0; 1] and f (0) = 0, f (1) = 1.For the convection-diffusion operator in (1), it is natural to set the diffusion operator as the initial operator.Here we consider a simple linear function f (t) = t to design the sequence of operators used for the homotopy.Given a sequence of homotopy parameters, 0 = t 0 < t 1 < • • • < t L = 1, the homotopy operators with stochastic coefficients define a sequence of eigenvalue problems of the form for ℓ = 0, . . ., L. The following lemma [13, Lemma 4.1] establishes the homotopy error on the smallest eigenvalue in (46) for fixed ω.
Lemma 2. Suppose the velocity field a is divergence-free and ω is fixed.The homotopy error-which is defined as the difference between the smallest eigenvalue λ(ω, t = 1) of the original operator and that of the homotopy operator in (46) satisfies for any t ∈ [0, 1] where and u * (ω, t) is the dual homotopy solution.For t sufficiently close to 1 and almost all ω ∈ Ω, C t,ω < C t for some C t < ∞ independent of ω.
Next, we bound C t,ω independently of ω.Clearly, the numerator is bounded for all t and almost all ω.Next, we show that the denominator is strictly positive.Suppose for a contradiction that u(ω, 1), u * (ω, t) = 0, then this implies that u(ω, 1), u * (ω, 1) − u * (ω, t) = u(ω, 1), u * (ω, 1) > 0, since the eigenfunction and dual eigenfunction are not orthogonal if the corresponding eigenvalues satisfy λ(ω, 1) = λ * (ω, 1).However, since u * (ω, t) → u * (ω, 1) as t → 1, the left hand side tends to zero whereas the right hand side is strictly positive and independent of t, leading to a contradiction.Hence, for t sufficiently small u(ω, 1), u * (ω, t) > 0 and similarly u(ω, t), u * (ω, 1) > 0. Thus, for t sufficiently small C t,ω < ∞.Since a(ω) along with the primal and dual eigenfunctions are continuous in ω, it follows that C t,ω is also continuous in ω and thus, can be bounded by the maximum over the compact domain Ω, With the homotopy method, the approximation error now comes from three sources: the FE discretization, the iterative eigensolver, and the value of the homotopy parameter.We suppose again that the error due to the eigensolver is bounded from above by the other two sources of error and design multilevel sequences such that the homotopy error and the discretization error are non-increasing with increasing level.Denoting the homotopy parameter and the mesh size at level ℓ by t ℓ and h ℓ , respectively, the multilevel sequence , and t L = 1.The multilevel parameters are required to be non-repetitive, i.e., (t ℓ−1 , h ℓ−1 ) = (t ℓ , h ℓ ) for all ℓ = 1, . . ., L, to ensure an asymptotically decreasing total approximation error in the sequence.However, one of these two parameters is allowed to be the same on two adjacent levels, i.e., either h ℓ−1 = h ℓ or t ℓ−1 = t ℓ is possible.This setting allows for adapting the homotopy parameter to discretisations on different meshes to satisfy the stability condition of the FE approximation.
The resulting MLMC estimator can be derived from the telescoping sum Following a similar derivation as that of Corollary 1 and based on the error bound in Lemma 2, we conjecture that the expectation and the variance of the multilevel difference with the homotopy method are bounded by respectively.This will be used as the guideline for choosing the multilevel sequences in our numerical experiments.We will also demonstrate that the above conjecture is valid in our numerical experiments.

Multilevel QMC Methods
QMC methods are a class of equal-weight quadrature rules originally designed to approximate high-dimensional integrals on the unit hypercube.A QMC approximation of the expected value of f is given by where, in contrast to Monte Carlo methods, the quadrature points {τ k } N −1 k=1 ⊂ [0, 1] s are chosen deterministically to be well-distributed and have good approximation properties in high dimensions.There are several types of QMC methods, including lattice rules, digital nets and randomised rules.The main benefit of QMC methods is that for sufficiently smooth integrands the quadrature error converges at a rate of O(N −1+δ ), δ > 0, or faster, which is better than the Monte Carlo convergence rate of O(N −1/2 ).For further details see, e.g., [20,21].
In this paper, we consider randomly shifted lattice rules, which are generated by a single integer vector z ∈ N s and a single random shift ∆ ∼ Uni[0, 1] s .The points are given by where {•} denotes taking the fractional part of each component.The benefits of random shifting are that the resulting approximation ( 52) is unbiased and that performing multiple QMC with i.i.d.random shifts provides a practical estimate for the mean-square error using the sample variance of the multiple approximations.If f is sufficiently smooth (i.e., has square-integrable mixed first derivatives) then a generating vector can be constructed such that the mean-square error (MSE) of a randomly shifted lattice rule approximation satisfies see, e.g., Theorem 5.10 in [20].I.e., for η ≈ 1/2 the convergence of the MSE is close to 1/N 2 .
Starting again with the telescoping sum (27), a multilevel QMC (MLQMC) method approximates the expectation of the smallest eigenvalue by using a QMC rule to compute the expectation on each level.MLQMC methods were first introduced in [32] for SDEs, then applied to parametric PDEs in [46,47] and elliptic eigenvalue problems in [28,29].For L ∈ N and {N ℓ } L ℓ=0 , the MLQMC approximation is given by where we apply a different QMC rule with points {τ ℓ,k } N ℓ −1 k=0 on each level, e.g., an N ℓ -point randomly shifted lattice rule (53) generated by z ℓ and an i.i.d.∆ ℓ .
The faster convergence of QMC rules leads to an improved complexity of MLQMC methods compared to MLMC, where in the best case the cost is reduced to close to ε −1 for a MSE of ε 2 .Following [46], under the same assumptions as in Theorem 3, but with Assumption II replaced by , the MSE of the MLQMC estimator ( 55) is bounded above by ε 2 and the cost satisfies The maximum level L is again given by ( 32) and {N ℓ } are given by where C ℓ is the cost per sample as in assumption III in Theorem 3 and N 0 is chosen as Verifying Assumption II(b) for the convection-diffusion EVP (1) requires performing a technical analysis similar to [28] and in particular, requires bounding the derivatives of the eigenvalue λ(ω) and its eigenfunction u(ω) with respect to ω.Such analysis is left for future work.In the numerical results, section we study the convergence of QMC and observe that II(b) holds with η ≈ 0.61.
In practice, one should perform multiple, say R ∈ N 0 , QMC approximations corresponding to i.i.d.random shifts, then take the average as the final estimate.In this way, we can also estimate the MSE by the sample variance over the different realisations.

Numerical results
In this section, we present numerical results for three test cases.The quantity of interest in all cases is the smallest eigenvalue of the stochastic convection-diffusion problem (1) in the unit domain D = [0, 1] 2 .The first two test cases use constant convection velocities at different magnitudes to benchmark the performance of eigenvalue solvers and finite element discretisation methods in the multilevel setting.In these two test cases, the random conductivity κ(x; ω) is modelled as a log-uniform random field constructed through the convolution of , where c i are the kernel centers placed uniformly on a 5 × 5 grid in the domain D. In the third test case, we also make the convection velocity a random field.Specifically, we first construct a log-uniform random field similar to that of the conductivity field using additional s a i.i.d.uniform random variables.Then, a divergence-free velocity field can be obtained by We employ the Eigen [34] library for Rayleigh quotient iteration and solve the linear systems using sparse LU decomposition with permutation from the SuiteSparse [18] library.For the implicitly restarted Arnoldi method, we use the ARPACK [49] library with the SM mode for finding the smallest eigenvalue.Random variables are generated using the standard C++ library and the pseudo-random seeds are the same across all experiments.
Numerical experiments are organized as follows.For a relatively low convection velocity a = [20; 0] T , we demonstrate the multilevel Monte Carlo (MLMC) method using the Galerkin FEM discretization.In this case, we also consider applying the homotopy method together with a geometrically refined mesh hierarchy.Then, on a test case with relatively high convection velocity a = [50; 0] T , we demonstrate the extra efficiency gain offered by the numerically more stable SUPG method, compared with the Galerkin discretization.For the third test case with a random velocity field, we apply SUPG to demonstrate the efficacy and efficiency of our multilevel method.Here we also demonstrate that quasi-Monte Carlo (QMC) samples can be used to replace Monte Carlo samples to further enhance the efficiency of multilevel methods.For all multilevel methods, we consider a sequence of geometrically refined meshes with h ℓ = h 0 × 2 −ℓ , ℓ = 0, 1, . . ., 4, and h 0 = 2 −3 .At the finest level, this gives 16129 degrees of freedom in the discretised linear system.We use 10 4 samples on each level ℓ to compute the estimates of rate variables α, β, γ in the MLMC complexity theorem (cf.Theorem 3).

Test case I
In the first experiment, we set a = [20; 0] T and use the Galerkin FEM to discretize the convection-diffusion equation.The stopping criteria for the Rayleigh quotient iteration and for the implicitly restarted Arnoldi method are set to be 10 −12 .In addition, for the implicitly restarted Arnoldi method, the Krylov subspace dimensions (the ncv values of ARPACK) are chosen empirically for each mesh size to optimize the number of Arnoldi iterations.They are m = 20, 40, 70, 70, 100 for h = 2 −3 , 2 −4 , 2 −5 , 2 −6 , 2 −7 , respectively.
We demonstrate the efficiency of four variants of the MLMC method: (i) the three-grid Rayleigh quotient iteration (tgRQI) with a model sequence defined by grid refinement; (ii) tgRQI with a model sequence defined by grid refinement and homotopy; (iii) the implicitly restarted Arnoldi method (IRAr) with a model sequence defined by grid refinement; and (iv) IRAr with a model sequence defined by grid refinement and homotopy.
(i) MLMC with tgRQI: Figure 2 illustrates the mean, the variance and the computational cost of multilevel differences λ ℓ (ω) − λ ℓ−1 (ω) of the smallest eigenvalue using tgRQI as the eigenvalue solver (without homotopy).Figure 2a also shows Monte Carlo estimates of the expected mean and variance of the smallest eigenvalue λ ℓ (ω) for each of the discretization levels.In addition to the computational cost, Figure 2b also shows the number of Rayleigh quotient iterations used at each level.We observe that the average number of iterations follows our analysis of the computational cost of tgRQI (cf.Alg. 2).From these plots, we estimate that the rate variables in the MLMC complexity theorem are α ≈ 2.0, β ≈ 4.0 and γ ≈ 2.41.Since the variance reduction rate β is larger than the cost increase rate γ, the MLMC estimator is in the best case scenario, with O(ε −2 ) complexity, as stated in Theorem 3.
(ii) MLMC with homotopy and tgRQI: Next, we consider the homotopy method in the MLMC setting together with tgRQI.We use the conjecture in (51) to set the (a) Means and variances of λ ℓ and λ ℓ − λ ℓ−1 .homotopy parameters such that 1 − t ℓ = O(h 2 ℓ ), t 0 = 0 and t L = 1.For L = 5, this results in t ℓ = {0, 3/4, 15/16, 63/64, 1}.With this choice the eigenproblem on the zeroth level contains no convection term and is thus self-adjoint.Figure 3a shows again the means and the variances of the multilevel differences λ ℓ − λ ℓ−1 in this setting, together with MC estimates of the expected means and variances of the eigenvalues for each level.The hierarchy of homotopy parameters is chosen to guarantee good variance reduction for MLMC.Indeed, the variance of the multilevel difference decays smoothly with a rate β ≈ 3.65.The expected mean of the difference, on the other hand, stagnates between ℓ = 1 and ℓ = 2.However, this initial stagnation is irrelevant for the MLMC complexity theorem; eventually for ℓ ≥ 2, the estimated means of the multilevel differences decrease again with a rate of α ≈ 2. Figure 3b shows the number of Rayleigh quotient iterations used at each level and the computational cost, which grows with a rate of γ ≈ 2.56 here.This leads to the same asymptotic complexity for MLMC, since the regime is the same, i.e., β > γ, which is the optimal regime in Theorem 3 with a complexity of O(ε −2 ).
(iii) MLMC with IRAr: Similar results are obtained by using the implicitly restarted Arnoldi eigenvalue solver (without homotopy).Since the mean and the variance of the    multilevel differences in this setting are almost identical to those of the Rayleigh quotient solver, we omit the plots here and only report the computational cost.Figure 4a shows the average number of matrix-vector products and the estimated CPU time for computing each of the multilevel differences, which grows with a rate of γ ≈ 3.5.Here, the increasing dimension of Krylov subspaces with grid refinement likely causes the higher growth rate of computational time compared to the experiment using tgRQI.Nonetheless, the MLMC estimator has again the optimal O(ε −2 ) complexity.
(iv) MLMC with homotopy and IRAr: Finally, we consider the behaviour of IRAr with homotopy, using the same sequence for the homotopy parameter t ℓ as in (ii).Again, we only focus on computational cost, showing the average number of matrix-vector products and the CPU time for computing each of the multilevel differences in Figure 4b.As in (ii), the cost grows at a rate of γ ≈ 3 leading again to the optimal O(ε −2 ) complexity for MLMC.
Overall comparison: In Figure 5, we show the CPU time versus the root mean square error for all four presented MLMC estimators together, as well as for standard Monte Carlo estimators using tgRQI (red) and IRAr (blue).The estimated complexity of standard Monte Carlo methods are O(ε −2.92 ) and O(ε −3.35 ) for tgRQI and IRAr, respectively.Overall, MLMC using tgRQI (without homotopy) outperforms all other methods, despite that all four MLMC methods achieve the optimal O(ε −2 ) complexity.Mesh level, ℓ (a) Means and variances of λ ℓ and λ ℓ − λ ℓ−1 .

Test case II
For the second experiment, we increase the velocity to a = [50; 0] T and focus on the comparison between Galerkin and SUPG discretizations.Thus, we only consider the threegrid Rayleigh quotient iteration (tgRQI) with a multilevel sequence based on geometrically refined grids without homotopy.Note that for such a strong convection, five steps in the homotopy approach are insufficient: the eigenvalues for consecutive homotopy parameters are too different to achieve variance reduction in the homotopy-based MLMC method.Its computational complexity is almost the same as the complexity of standard Monte Carlo, namely almost O(ε −3.5 ).The performance of MLMC with implicitly restarted Arnoldi on the other hand is similar to MLMC with tgRQI.Galerkin: Due to the higher convection velocity the first two levels are unstable for most of the realizations of ω as the FEM solution may exhibit non-physical oscillations.Thus, we set the coarsest level for the MLMC method to h 0 = 2 −5 here.Keeping the same finest grid level h L = 2 −7 , this means that we only use a total of three levels (L = 2) compared to the sequence in Test Case I, which had a total of five levels (L = 4).Figure 6a shows the expectation and variance of the multilevel differences.Here, we only have a couple of data points for estimating the rate variables of the MLMC complexity theorem, but the estimates are α ≈ 2 and β ≈ 4 as expected theoretically.The average number of Rayleigh quotient iterations in Figure 6b also behaves as in Test Case I with 5 iterations on the coarsest level and 2 iterations on the subsequent levels as expected for the three-grid Rayleigh quotient iteration (Alg.2) -recall that Levels 1 and 2 here correspond to Levels 3 and 4 in Figures 2b and 3b.The estimated value for γ ≈ 1.88, and thus the MLMC complexity is still O(ε −2 ).However, we cannot use as many levels due the numerical stability issues caused by the higher convection velocity, which substantially increases the prefactor in the O(ε −2 ) cost of the algorithm.SUPG: By using the SUPG discretization, we overcome the numerical stability issue and can use all five levels in MLMC, starting with h 0 = 2 −3 .As can be seen in Figure 7a, the expectation and the variance of the multilevel differences converge with the same rates as for the Galerkin FEM, namely α ≈ 2 and β ≈ 4 respectively.Also, clearly the use of SUPG leads to stable estimates even on the coarser levels.Figure 7b reports the average number of Rayleigh quotient iterations used at each level and the computational cost.We estimate that the computational cost increases at a rate of γ ≈ 2.33 here.In any case, the    use of SUPG in the MLMC also results in the optimal O(ε −2 ) complexity.
Overall comparison: Figure 8 shows CPU times versus root mean square errors for the MLMC methods (with tgRQI and without homotopy) using Galerkin FEM and SUPG discretizations.They are compared to a standard Monte Carlo method with Galerkin FEM.Although both MLMC estimates have the optimal O(ε −2 ) complexity, the stability offered by SUPG enables us to use more, coarser levels, thus leading to a smaller prefactor and a significant computational gain of a factor 10-20 over the Galerkin FEM based method.

Test Case III
In this experiment, the convection velocity becomes a divergence-free random field generated using ( 57) and (58).We discretise the eigenvalue problem using SUPG and apply the three-grid Rayleigh quotient iteration (tgRQI) without homotopy to solve multilevel eigenvalue problems.The stopping criteria for tgRQI is set to be 10 −12 .The same sequence of grid refinements, h = 2 −3 , 2 −4 , 2 −5 , 2 −6 , 2 −7 , as in previous test cases is used to construct multilevel estimators.
MLMC: Figure 9 illustrates the mean, the variance and the computational cost of multilevel differences λ ℓ (ω)− λ ℓ−1 (ω) of the smallest eigenvalue using tgRQI as the eigenvalue solver.Figure 9a also shows Monte Carlo estimates of the expected mean and variance of the smallest eigenvalue λ ℓ (ω) for each of the discretization levels.In addition to the  computational cost, Figure 9b also shows the number of Rayleigh quotient iterations used at each level.We observe that the average number of iterations follows our analysis of the computational cost of tgRQI (cf.Alg. 2).From these plots, we estimate that the rate variables in the MLMC complexity theorem are α ≈ 2.0, β ≈ 4 and γ ≈ 2.23.Since the variance reduction rate β is larger than the cost increase rate γ, the MLMC estimator is in the best case scenario, with O(ε −2 ) complexity, as stated in Theorem 3. In Figure 11, we compare the computational complexity of MLMC to that of the standard Monte Carlo.Numerically, we observe that the CPU time of MLMC is approximately O(ε −2.06 ), which is close to the theoretically predicted rate.In comparison, the CPU time of the standard MC is approximately O(ε −3.2 ) in this test case.
MLQMC: All QMC computations were implemented using Dirk Nuyens' code accompanying [45] and use a randomly shifted embedded lattice rule in base 2, as outlined in [16], with 32 i.i.d.random shifts.In Figure 10, we plot convergence of the MSE for both MC and QMC for three different cases: for λ 0 in plot (a), for the difference λ 1 − λ 0 in plot (b), and for the difference λ 2 − λ 1 in plot (c).Here the meshwidths are given by h 0 = 2 −3 , h 1 = 2 −4 and h 2 = 2 −5 .In all cases, QMC outperforms MC, where for λ 0 the MSE for QMC converges at an observed rate of −1.78, whereas MC converges with the rate −1.For the other two cases, which are MSEs of multilevel differences, the QMC converges with an approximate rate of −1.63, which is again clearly faster than the MC convergence rate of −1.This observed MSE convergence for the QMC approximations of the differences implies that II(b) holds with η ≈ 0.61.For MLQMC, to choose N ℓ we use (56) with η ≈ 0.61 and with N 0 scaled such that the overall MSE is less than ε 2 / √ 2 for each tolerance ε.Since we use a base-2 lattice rule, we round up N ℓ to the next power of 2.
The MLQMC complexity, in terms of CPU time, is plotted in Figure 11, along with the results for MC and MLMC.Comparing the three methods in Figure 11, clearly MLQMC provides the best complexity, followed by MLMC then standard MC.In this case, we have the approximate rates βη ≈ 4×0.61 = 2.44 > γ ≈ 2.23, which implies that for MLQMC we are in the optimal regime for the cost with C MLQMC (ε) ε −2η .Numerically, we observe that the rate is given by 1.28, which is very close to the theoretically predicted rate of 2η ≈ 1.22.

Conclusion
In this paper we have considered and developed various MLMC methods for stochastic convection-diffusion eigenvalue problems in 2D.First, we established certain error bounds on the variational formulation of the eigenvalue problem under assumptions such as eigenvalue gap, boundedness, and other approximation properties.Then we presented the MLMC method based on a hierarchy of geometrically refined meshes with and without homotopy.We also discussed how to improve the computational complexity of MLMC by replacing Monte Carlo samples with QMC samples.At last, we provided numerical results for three test cases with different convection velocities.
Test Case I shows that, for low convection velocity, all variants of the MLMC method (based on a Galerkin FEM discretization of the PDE) achieve optimal O(ε −2 ) complexity, including the one with homotopy.In Test Case II with a high convection velocity, the homotopy-based MLMC does not work anymore -at least without increasing the number of levels -and MLMC based on Galerkin FEM has severe stability restrictions, preventing the use of a large number of levels.This restriction can be circumvented easily by using stable SUPG discretizations.Numerical experiments suggest that MLMC with SUPG achieves the optimal O(ε −2 ) complexity and is 10-20 times faster than the Galerkin FEMbased versions for the same level of accuracy.In Test Case III, we considered both the conductivity and the convection velocity as random fields and compared the performance of MLMC and MLQMC.In this example, both MLMC and MLQMC deliver computational complexities that are close to the optimal complexities predicted by the theory, while the rate of the computational complexity of MLQMC outperforms that of MLMC.
7 Appendix: Bounding the constants in the FE error The results in Theorem 2 follow from the Babuška-Osborn theory [4].In this appendix we show that the constants can be bounded independently of the stochastic parameter.
The Babuška-Osborn theory studies how the continuous solution operators T ω , T * ω : V → V , which for f, g ∈ V are defined by A(ω; T ω f, v) = f, v for all v ∈ V, A(ω; w, T * ω g) = w, g for all w ∈ V, are approximated by the discrete operators T ω,h , T * ω,h : V h → V h , A(ω; T ω,h f, v h ) = f, v h for all v h ∈ V h , A(ω; w h , T * ω,h g) = w h , g for all w h ∈ V h .We summarize the pertinent details here.First, we introduce: The result for the eigenfunction ( 17) is given by [4, Thm.8.1], which gives for a constant C(ω) defined below.Since λ(ω) is simple, the best approximation property of V h in H 2 (D) followed by Theorem 1 gives where the best approximation constant C BAP is independent of ω.In the last inequality we have also used that λ(ω) is continuous on the compact domain Ω, thus can be bounded uniformly by λ := max is independent of ω.
For the eigenvalue error (16) we follow the proof of [4,Theorem 8.2].Since λ(ω) is simple, from Theorem 7.2 in [4], the eigenvalue error is bounded by where in the second inequality we have used (60) and the equivalent bound for the dual eigenvalue, combining the two constants into C η .By following [4], the constant C λ (ω) can be bounded independently of ω in a similar way to C u (ω).

Figure 2 -
Figure 2 -MLMC method using tgRQI for Test Case I with a = [20; 0] T and Galerkin FEM: (a) Mean (blue) and variance (red) of the eigenvalue λ ℓ (dashed) and of λ ℓ − λ ℓ−1 (solid); (b) computational times for one multilevel difference (blue) and average number of Rayleigh quotient iterations (red) on each level.Where shown, the error bars represent ± one standard deviation.

Figure 3 -
Figure 3 -MLMC method using homotopy and tgRQI for Test Case I with a = [20; 0] T and Galerkin FEM: (a) Mean (blue) and variance (red) of the eigenvalue λ ℓ (dashed) and of λ ℓ − λ ℓ−1 (solid); (b) computational times for one multilevel difference (blue) and average number of RQIs (red) on each level.Where shown, the error bars represent ± one standard deviation.

Figure 4 -
Figure 4 -MLMC method using IRAr for Test Case I with a = [20; 0] T and Galerkin FEM, both without (a) and with (b) homotopy: average computational cost (blue) and average number of matrix-vector products (red) per sample of λ ℓ − λ ℓ−1 .The error bars represent ± one standard deviation.

Figure 5 -
Figure 5 -CPU time vs. root mean square error of all estimators in Test Case I.
Computational times and average RQI.

Figure 6 -
Figure 6 -MLMC method using tgRQI for Test Case II with a = [50; 0] T and Galerkin FEM: (a) Mean (blue) and variance (red) of the eigenvalue λ ℓ (dashed) and of λ ℓ − λ ℓ−1 (solid); (b) computational time for one multilevel difference (blue) and average number of Rayleigh quotient iterations (red) on each level.Where shown, the error bars represent ± one standard deviation.

Figure 7 -
Figure 7 -MLMC method using tgRQI for Test Case II with a = [50; 0] T and SUPG discretization: (a) Mean (blue) and variance (red) of the eigenvalue λ ℓ (dashed) and of λ ℓ − λ ℓ−1 (solid); (b) computational time for one multilevel difference (blue) and average number of Rayleigh quotient iterations (red) on each level.Where shown, the error bars represent ± one standard deviation.

Figure 8 -
Figure 8 -CPU time vs. root mean square error of the estimators in Test Case II.

Figure 9 -
Figure 9 -MLMC method using tgRQI and SUPG for Test Case III with random velocity and random conductivity: (a) Mean (blue) and variance (red) of the eigenvalue λ ℓ (dashed) and of λ ℓ − λ ℓ−1 (solid); (b) computational times for one multilevel difference (blue) and average number of RQIs (red) on each level.Where shown, the error bars represent ± one standard deviation.

Figure 10 -
Figure 10 -Convergence of QMC and MC methods using tgRQI and SUPG for Test Case III with random velocity and conductivity.Plots (a), (b), (c) give the MSE of estimators versus sample sizes for grid sizes h = 2 −3 , 2 −4 , 2 −5 , respectively.Blue lines with circles and black lines with squares indicate the MSE for MC and QMC, respectively.Dashed lines and solid lines correspond to the MSE of the estimated multilevel differences and the MSE of the estimated eigenvalues, respectively.

2 − 6 2Figure 11 -
Figure 11 -CPU time vs. root mean square error of the estimators in Test Case III.