Numerical solution of fractional elliptic stochastic PDEs with spatial white noise

The numerical approximation of solutions to stochastic partial differential equations with additive spatial white noise on bounded domains in $\mathbb{R}^d$ is considered. The differential operator is given by the fractional power $L^\beta$, $\beta\in(0,1)$, of an integer order elliptic differential operator $L$ and is therefore non-local. Its inverse $L^{-\beta}$ is represented by a Bochner integral from the Dunford-Taylor functional calculus. By applying a quadrature formula to this integral representation, the inverse fractional power operator $L^{-\beta}$ is approximated by a weighted sum of non-fractional resolvents $(I + t_j^2 L)^{-1}$ at certain quadrature nodes $t_j>0$. The resolvents are then discretized in space by a standard finite element method. This approach is combined with an approximation of the white noise, which is based only on the mass matrix of the finite element discretization. In this way, an efficient numerical algorithm for computing samples of the approximate solution is obtained. For the resulting approximation, the strong mean-square error is analyzed and an explicit rate of convergence is derived. Numerical experiments for $L=\kappa^2-\Delta$, $\kappa>0$, with homogeneous Dirichlet boundary conditions on the unit cube $(0,1)^d$ in $d=1,2,3$ spatial dimensions for varying $\beta\in(0,1)$ attest the theoretical results.


Introduction
A real-valued Gaussian random field u defined on a spatial domain D ⊂ R d is called a Gaussian Matérn field if its covariance function C : D × D → R is given by can be controlled: its variance, smoothness, and correlation range. Due to this flexibility, Gaussian Matérn fields are often used for modeling in spatial statistics, see, e.g., Stein (1999). However, a major drawback of this traditional covariancebased representation of Matérn fields is its high computational effort. For instance, sampling from u at n locations x 1 , . . . , x n ∈ D requires a matrix factorization of an n × n covariance matrix and, thus, in general, O(n 3 ) arithmetic operations. There are several approaches attempting to cope with this computational problem (see, e.g., Banerjee et al., 2008;Furrer et al., 2006;Nychka et al., 2015;Sun et al., 2012). One of them is based on the insight that a Gaussian Matérn field on D := R d with parameters σ, ν, κ > 0 can be seen as the statistically stationary solution u to the stochastic partial differential equation (SPDE) x ∈ D. (1.2) Here, ∆ denotes the Laplacian, 4β = 2ν + d, and W is Gaussian white noise on R d . The marginal variance of u is then given by σ 2 = Γ(2β−d/2)Γ(2β) −1 (4π) −d/2 κ d−4β . This relation between Matérn fields and SPDEs had already been noticed by Whittle (1954Whittle ( , 1963. Later on, based on a finite element discretization of (1.2), where the differential operator κ 2 − ∆ is augmented with Neumann boundary conditions, Markov random field approximations of Matérn fields on bounded domains D R d were introduced by Lindgren et al. (2011). Owing to the computational savings of this approach compared to the covariance-based approximations, these SPDE-based models have become very popular in spatial statistics (see, e.g., Bhatt et al., 2015;Lai et al., 2015), and they are still subject of current research, mainly because of the following reason: the SPDE formulation (1.2) facilitates various generalizations of approximations of Gaussian Matérn fields involving (a) other differential operators (Bolin & Lindgren, 2011;Fuglstad et al., 2015;Lindgren et al., 2011), (b) more general domains, such as the sphere (Bolin & Lindgren, 2011;Lindgren et al., 2011), and (c) non-Gaussian driving noise (Bolin, 2014;Wallin & Bolin, 2015). However, a considerable drawback of the finite element approximation proposed by Lindgren et al. (2011) is that it is only computable if 2β ∈ N. This limits flexibility, and, in particular, it implies that the method cannot be applied to the important special case of exponential covariance (ν = 1/2) in R 2 , which corresponds to β = 3/4.
With the objective of extending the approach of Lindgren et al. (2011) to all admissible values ν > 0 (i.e., β > d/4) and the generalizations mentioned above, we consider (1.2) in the more general framework of fractional order elliptic equations driven by white noise. We propose an explicit numerical scheme for generating samples of an approximation to the Gaussian solution process, which allows for any fractional power β > d/4. This method is based on: (i) a standard finite element discretization in space, (ii) a quadrature approximation of the inverse fractional elliptic differential operator proposed by Bonito & Pasciak (2015) for deterministic equations, and (iii) an approximation of the noise term on the right-hand side, whose covariance matrix after discretization is equal to the finite element mass matrix. Due to (iii) no explicit knowledge of the eigenfunctions of the differential operator is needed in contrast to approximations based on truncated Karhunen-Loève expansions of the noise term (e.g., Zhang et al., 2016).
While Zhang et al. (2016) remarked that any orthonormal basis can be used in the Karhunen-Loève expansion of the noise term before projecting it onto the finite element space, their error analysis strongly benefits from the fact that the eigenfunctions of the differential operator are used. In fact, if a general orthonormal basis was used, it would be difficult to obtain explicit rates of convergence with respect to the truncation level after truncating the expansion at a finite level. Only if the constructed basis has certain smoothness with respect to the differential operator, it is possible to obtain explicit rates of convergence (see, e.g., Kovács et al., 2011, addressing this kind of problem). Constructing such an orthonormal basis requires a lot of computational effort, in particular, for complicated domains and d = 2, 3. In contrast, the present approach is based on an expansion with respect to the standard (non-orthonormal) finite element basis for the discretized problem, which is readily available even for complex geometries.
This work follows a series of investigations of numerical methods for deterministic fractional order equations (Baeumer et al., 2015;Bonito et al., 2017;Bonito & Pasciak, 2015;Caffarelli & Silvestre, 2007;Gavrilyuk, 1996;Gavrilyuk et al., 2004Gavrilyuk et al., , 2005Jin et al., 2015;Nochetto et al., 2015;Roop, 2006), and for non-fractional elliptic SPDEs with random forcing (e.g., Babuska et al., 2004;Cao et al., 2007;Du & Zhang, 2003;Gyöngy & Martínez, 2006;Zhang et al., 2016). An essentially similar idea to our approach, which combines (iii) with Lanczos and Krylov subspace methods, was persued by Simpson (2008) for the differential operator κ 2 − ∆ in (1.2). However, neither weak nor strong convergence have been proven and the empirical results show a generally poor performance of the proposed scheme. A weak error estimate for κ 2 − ∆ and the non-fractional case β = 1 in d = 2 spatial dimensions has been derived by Simpson et al. (2012), showing quadratic weak convergence for that particular case. Since the first two moments uniquely determine the distribution of the Gaussian solution process, an alternative approach is to consider the problem of approximating the covariance function of the solution instead of solving for the solution process itself (Dölz et al., 2017). Note that this method cannot be generalized to non-Gaussian models.
The structure of this article is as follows: In §2 we present the fractional order equation with Gaussian white noise in a Hilbert space setting, the necessary assumptions, and comment on existence and uniqueness of a solution to this problem. Furthermore, we introduce the numerical approximation of the solution process and state our main result in Theorem 2.10: strong mean-square convergence of this approximation at an explicit rate. This theorem is proven in §3 by partitioning the strong mean-square error in several terms and estimating these terms one by one. In addition, a weak type convergence result is obtained in Corollary 3.4. In §4 the SPDE (1.2) is considered for numerical experiments on the unit cube D = (0, 1) d with continuous, piecewise linear finite element basis functions in d = 1, 2, 3 spatial dimensions. The outcomes of the paper are summarized in §5.

Model problem and main result
In the following, let H denote a separable Hilbert space and L : D(L) ⊂ H → H be a densely defined, self-adjoint, positive definite linear operator with a compact inverse. In this case, there exists an orthonormal basis {e j } j∈N of H consisting of eigenvectors of L. The eigenvalue-eigenvector pairs {(λ j , e j )} j∈N can be arranged such that the positive eigenvalues {λ j } j∈N are in nondecreasing order: We assume that the growth of the eigenvalues is given by λ j ∝ j α for an exponent α > 0, i.e., there exist constants c λ , C λ > 0 such that (2.1) For β > 0 and φ ∈ D(L β ) := {ψ ∈ H : j∈N λ 2β j (ψ, e j ) 2 H < ∞} the action of the fractional power operator L β : D(L β ) → H is defined by (2. 2) The subspaceḢ 2β := D(L β ) ⊂ H is itself a Hilbert space with respect to the inner product and corresponding norm given by Its dual space after identification via the inner product on H is denoted byḢ −2β . For s ≥ 0, the norm on the dual spaceḢ −s enjoys the useful representation where ·, · denotes the duality pairing betweenḢ −s andḢ s (Thomee, 2007, Proof of Lem. 5.1). We obtain the following scale of densely and continuously embedded Hilbert spaces: It is an immediate consequence of these definitions that L β is an isometric isomorphism fromḢ s toḢ s−2β for s ≥ 2β, since for φ ∈Ḣ s it holds The following lemma states that L β can be extended to a bounded linear operator betweenḢ s andḢ s−2β for all s ∈ R.
Lemma 2.1. For s ∈ R, there exists a unique continuous extension of L β defined in (2.2) to an isometric isomorphism L β :Ḣ s →Ḣ s−2β .
Proof. For s ≥ 2β the isometry property and, thus, injectivity has already been observed in (2.5). Surjectivity readily follows, since for any g ∈Ḣ s−2β the vector φ := j∈N λ −β j (g, e j ) H e j is an element ofḢ s with φ s = g s−2β and L β φ = g. Assume now that s < 2β. We obtain for φ ∈Ḣ 2β and ψ ∈Ḣ 2β−s By density ofḢ 2β →Ḣ s , it follows that there exists a unique linear continuous extension L β :Ḣ s →Ḣ s−2β . Moreover, it is an isometry, since the above estimate attains equality for φ ∈Ḣ s and ψ = L s−β φ ∈Ḣ 2β−s . Surjectivity follows similarly to s ≥ 2β from the representation of the dual norm in (2.3) applied toḢ s−2β . ] s,q denotes the real K-interpolation method (see, e.g., Thomee, 2007, Ch. 19). Furthermore, one can show (Bonito & Pasciak, 2015, Prop. 4.1) that for 1 ≤ s ≤ 2 it holds thaṫ . If D = (0, 1) d is the d-dimensional unit cube, then it is well-known (e.g., Courant & Hilbert, 1962, Ch. VI.4) that the above eigenvalue problem has the following eigenvectors The corresponding eigenvalues are given by These eigenvalues satisfy (2.1) for α = 2/d. Note that only the values of c λ and C λ in (2.1) change when considering any other bounded domain D ⊂ R d with smooth or polygonal boundary. The value α = 2/d is the same by the min-max principle.
2.1. Fractional order equation. Motivated by (1.2) we consider for g ∈ H the fractional order equation where W denotes Gaussian white noise defined on a complete probability space (Ω, A, P) with values in the separable Hilbert space H. We assume that β ∈ (0, 1) and refer to Remark 3.6 for a discussion of the generalization to β > 0. Equation (2.8) as well as all following equalities involving noise terms are understood to hold P-almost surely (P-a.s.). Note that the white noise W can formally be represented by the Karhunen-Loève expansion with respect to the orthonormal eigenbasis {e j } j∈N ⊂ H of L: where {ξ j } j∈N is a sequence of independent real-valued standard normally distributed random variables. As we will show in Proposition 2.3, this formally defined series converges in L 2 (Ω;Ḣ −s ) for any s > 1/α, which implies that realizations of W are elements ofḢ − 1 α − for any > 0 P-a.s. Related to this representation, we introduce for N ∈ N the truncated white noise (2.10) The following proposition specifies mean-square regularity of the white noise with respect to theḢ-spaces in (2.4).
Proposition 2.3. For all s ≥ 0, there exists a constant C > 0 depending only on α, c λ in (2.1), and s, such that the truncated white noise W N in (2.10) satisfies Proof. The orthonormality of the vectors {e j } along with the distribution of the random variables {ξ j } in the expansion of W N and the representation ( Choosing s = α −1 + and taking the limit N → ∞ shows the L 2 Ω;Ḣ − 1 α − regularity of W. Remark 2.4. Proposition 2.3 implies that W ∈Ḣ − 1 α − holds also P-a.s. for any > 0. Therefore, existence and uniqueness of a solution u to (2.8) with regularity u ∈Ḣ 2β− 1 α − (P-a.s.) follow from Lemma 2.1. In addition, the above results show that u ∈ L 2 Ω;Ḣ 2β− 1 α − . In particular, u ∈ L 2 (Ω; H) holds if 2αβ > 1.
Remark 2.5. The above results are in accordance with Lemma 4.2 and Theorem 4.3 by Zhang et al. (2016), where the following semilinear elliptic stochastic boundary value problem is considered on D := (0, 1) d for d ∈ {1, 2, 3}, g ∈ L 2 (D), and f : R → R: and mean-square regularity in the Sobolev space H s (D) is proven under appropriate assumptions on the nonlinearity f for (i) the spatial white noise W and s < −d/2 and (ii) the solution u and s < 2−d/2. Note that in the linear case when f ≡ κ 2 ≥ 0, this corresponds to Problem (2.8) with β = 1, L = κ 2 − ∆, and the exponent α of the eigenvalue growth in (2.1) is given by α = 2/d, see Example 2.2.
2.2. Finite element approximation. In the following, we introduce a numerical method based on a finite element discretization for solving the fractional order equation (2.8) approximately. For this purpose, we consider the family (V h ) h∈(0,1) of subspaces ofḢ 1 with finite dimensions N h := dim(V h ) < ∞. The Galerkin discretization of the operator L : These eigenvalues are again arranged in nondecreasing order: Further assumptions on the approximation properties of finite element spaces are specified below.
Assumption 2.6. The family (V h ) h∈(0,1) of finite-dimensional subspaces ofḢ 1 satisfies the following: The following example illustrates that Assumption 2.6 is, in general, satisfied for d-dimensional elliptic linear differential operators.
Example 2.7. Let D ⊂ R d be a bounded, convex, polygonal domain and assume that L : D(L) ⊂ L 2 (D) → L 2 (D) is a (strongly) elliptic linear differential operator of order 2m ∈ N, i.e., there exists a constant γ > 0 such that , and ·, · denotes the duality pairing between V and its dual V * after identification via the inner product on L 2 (D). Assume that V h ⊂ V is an admissible finite element space of polynomial degree p ∈ N with respect to a regular mesh onD. In this case, Assumption 2.6 is satisfied (Strang & Fix, 2008, Thms. 6.1 & 6.2) for H = L 2 (D),Ḣ 1 = V , and the exponents r = 2(p + 1 − m), s = min{p + 1, 2(p + 1 − m)}, and q = p+1 m . In particular, if the family of finite element spaces V h is quasi-uniform and has continuous, piecewise linear basis functions, we have for L = κ 2 − ∆ with homogeneous Dirichlet boundary conditions in Example 2.2 that r = s = q = 2.
In order to approximate the white noise W on the finite element space V h we introduce the following V h -valued random variables: (i) an expansion with respect to the discrete eigenbasis E: 14) where ξ := (ξ 1 , . . . , ξ N h ) T is the vector of the first N h independent standard normally distributed random variables in expansion (2.9) of W; (ii) an expansion with respect to any basis Φ : where the random vector ξ := ( ξ 1 , . . . , ξ N h ) T is given by It is therefore Gaussian distributed with zero mean and covariance matrix is the mass matrix with respect to the basis Φ of the finite element space V h .
The following lemma shows that the above approximations of the white noise are equal in mean-square sense.
Proof. Inserting the definitions (2.14)- By definition of the matrices R and M and due to the distribution of the random vectors ξ and ξ we have for i, j ∈ {1, . . . , N h }: where δ ij is the Kronecker delta. Thus, in terms of the trace tr of matrices we have Our numerical approach to cope with the fractional order equation (2.8) will be based on the following representation of the inverse L −β from the Dunford-Taylor calculus due to Balakrishnan (1960), see also (Yosida, 1995, §IX.11, Eq. (4)): −∞ e 2βy (I + e 2y L) −1 dy.
We choose an equidistant grid {y = k : ∈ Z, −K − ≤ ≤ K + } with step size k > 0 for y and replace the differential operator L with its discrete version L h to formulate the following quadrature method proposed by Bonito & Pasciak (2015): Exponential convergence of order O(e −π 2 /(2k) ) to L −β h with respect to the norm Note that the quadrature (2.16) involves only non-fractional resolvents of the discrete operator L h . Thus, the corresponding numerical method is readily implementable.
With the notions of the H-orthogonal projection Π h on V h , the noise approximation W Φ h in (2.15) and the quadrature Q β h,k in (2.16) at hand, we can now introduce the numerical approximation u Q h,k of the solution u to (2.8) as Remark 2.9. We emphasize the construction of the noise term W Φ h on the righthand side of (2.18): For discretizing the white noise W, it is common to project the truncated Karhunen-Loève expansion W N in (2.10) onto the finite element space V h and use the noise approximation Π h W N ∈ L 2 (Ω; V h ) (e.g., Zhang et al., 2016). Instead, we define W Φ h in (2.15) in such a way that it is mean-square equivalent to the noise approximation W E h in (2.14), which can be interpreted as V h -valued white noise expanded with respect to the H-orthonormal eigenvectors of L h . Note that the noise approximation W E h is needed only for the error analysis, but not for the actual implementation of the numerical algorithm.
This approach has the following two major advantages in practice: (i) Samples of the truncated white noise W N are not needed and, thus, neither are the exact eigenvectors of the operator L. (ii) For the computation of the approximation u Q h,k in (2.18) in practice, one has to sample from the load vector b with entries is a basis of the finite element space V h . We emphasize that this is computationally feasible if the basis Φ in the noise approximation W Φ h is chosen as the same, since then b ∼ N (g, M), where g i = (g, φ i,h ) H and M is the mass matrix with respect to the basis Φ, which is usually sparse. Hence, samples of b can be generated from b = g + Gz, where z ∼ N (0, I) and G is a Cholesky factor of M = GG T .
The following estimate of the strong mean-square error between the exact solution u to the fractional order equation (2.8) and the numerical approximation u Q h,k in (2.18) is our main result.
Theorem 2.10. Let the family (V h ) h∈(0,1) of finite element spaces satisfy Assumption 2.6 and assume that the growth of the eigenvalues of the operator L is given by (2.1) for an exponent α with where the values of d ∈ N, r, s > 0, and q > 1 are the same as in Assumption 2.6. Then, for sufficiently small h ∈ (0, h 0 ) and k ∈ (0, k 0 ), the strong L 2 (Ω; H) error between the solution u of (2.8) and the approximation u Q h,k in (2.18) is bounded by (2.20) where the constant C > 0 is independent of h and k.

Partition of the error and error estimates
In order to prove Theorem 2.10 we express the difference between the exact solution u and the approximation u Q h,k as follows: , and partition the strong error in (2.20) accordingly. Here, u N h , u E h , and u Φ h are defined in terms of the truncated white noise in (2.10) and the white noise approximations in (2.14)-(2.15) as the P-almost sure solutions to the following equations: and u N h refers to the truncation with N = N h = dim(V h ) terms in (3.1). Note that by Lemma 2.8 the difference between u E h and u Φ h vanishes identically in L 2 (Ω; H): In the following, we address the three remaining terms separately: the truncation error, the error of the finite element discretization, and the error caused by the quadrature approximation Q β h,k in (2.16) of the discrete fractional inverse L −β h . 3.1. Truncation error. In the lemma below, we bound the strong mean-square error between the exact solution u to the fractional order equation (2.8) and the approximation u N in (3.1) based on the truncated right-hand side g N + W N .
Lemma 3.1. Assume that the eigenvalues of the operator L satisfy (2.1) and that 2αβ > 1. Let u ∈ L 2 (Ω; H) be the solution to (2.8). Then, for any g ∈ H, N ∈ N, there exists a unique solution u N ∈ L 2 (Ω;Ḣ 2β ) to the truncated equation (3.1) and it satisfies where the constant C > 0 depends only on α, β, and the constants in (2.1). In particular, under Assumption 2.6 it holds that u − u N h L2(Ω;H) h d(αβ−1/2) (1 + g H ). (3.5) Proof. Existence and uniqueness of a solution u N in L 2 (Ω;Ḣ 2β ) follows from the fact that for any g ∈ H the truncated right-hand side g N + W N is an element of L 2 (Ω; H) as well as the isomorphism property of L β :Ḣ 2β → H in Lemma 2.1. If 2αβ > 1, we obtain for any N ∈ N: Under Assumption 2.6(i) we have N h ∝ h −d and (3.5) readily follows.

3.2.
Finite element discretization error. The next ingredient in the derivation of (2.20) in Theorem 2.10 is an upper bound for the error caused by introducing the finite element element discretization. More precisely, the solution u N of the truncated problem (3.1) corresponds to the best V N -valued approximation of u ∈ L 2 (Ω; H). Here, the finite-dimensional subspace V N ⊂Ḣ 1 is the linear span of the first N eigenvectors of the operator L. Subsequently, the approximation u E h has been defined in (3.2), which takes values in another finite-dimensional subspace, namely in the finite element space V h ⊂Ḣ 1 .
The purpose of the following lemma is to bound the error between these two approximations when N = dim(V h ).
Lemma 3.2. Let Assumption 2.6 be satisfied. Assume that the eigenvalue growth of the operator L is given by (2.1) for an exponent α with (2.19). Let u E h be the unique element in L 2 (Ω; V h ) satisfying (3.2), i.e., and let u N h ∈ L 2 (Ω; H) denote the solution to the truncated equation (3.1) with N = N h = dim(V h ) terms. Then their difference can be bounded by for sufficiently small h ∈ (0, h 0 ), where the constant C > 0 depends only on α, β, and the constants in (2.1), (2.12), and (2.13).
Proof. In order to show the estimate in (3.6), we first note that Π h g can be expanded in terms of the orthonormal eigenvectors {e j,h } of L h by Thus, we can partition the difference in (3.6) as follows: Since the random variables {ξ j } are independent and standard normally distributed, we obtain for the first term by the Cauchy-Schwarz inequality for sums Assumption 2.6(i), (2.1), and (2.13) complete the estimation of the first term for h ∈ (0, h 0 ), since since dαq ≤ 2s by (2.19). For the second term we find By the mean value theorem, there exists λ j ∈ (λ j , λ j,h ) such that Therefore, we can bound the third term by (2.1) and (2.12) as follows: where we used the relations N h ∝ h −d from Assumption 2.6(i) and αd(q − 1) ≤ r from (2.19) in the last estimate.
3.3. Quadrature approximation error. As the final step for proving the strong convergence result of Theorem 2.10, we investigate the error caused by applying the quadrature Q β h,k in (2.16) instead of the discrete fractional inverse L −β h to the right-hand side Π h g + W Φ h . The difference between these two operators in the norm (2.17) on the space L(V h ) has been bounded by Bonito & Pasciak (2015). The following lemma is a consequence of that result and the distribution of the noise approximation W Φ h . Lemma 3.3. Let Q β h,k : V h → V h be the operator in (2.16) approximating L −β h . Then it holds for sufficiently small k ∈ (0, k 0 ) and any φ h ∈ V h that where the constant C > 0 depends on β and grows linearly in the largest eigenvalue (3.8) Proof. The first assertion is proven in (Bonito & Pasciak, 2015, Lem. 3.4, Thm. 3.5).
Proof of Theorem 2.10 and some remarks. After having bounded the truncation, the discretization, and the quadrature errors in § §3.1-3.3, the strong convergence result of Theorem 2.10 is an immediate consequence.
Proof of Theorem 2.10. As outlined at the beginning of the section, we partition the mean-square error as follows: (II) h min{d(αβ−1/2),r,s} (1 + g H ), and by (3.4) the third term vanishes, (III) = 0. Thus, the assertion is proven.
In §4 we will not only verify the above derived rate of strong convergence by means of numerical experiments, but also investigate weak type errors. The following result is proven similarly to Theorem 2.10.
Corollary 3.4. Let the assumptions of Theorem 2.10 be satisfied. Then there exists a constant C > 0 independent of h ∈ (0, h 0 ) and k ∈ (0, k 0 ) such that the following weak type error estimate between the approximation u Q h,k in (2.18) and the solution u to (2.8) holds: Proof. First we note that the norm on L 2 (Ω; H) attains the following values for u in (2.8), u N h in (3.1), and u Φ h in (3.3): Again, we partition the error,
Remark 3.5. The error estimates in (2.20) and (3.9) imply that the distance of the quadrature nodes k has to be adjusted to the finite element mesh size h. Table 1 shows the calibration between h and k for error studies of strong and weak type as well as the corresponding theoretical convergence rates. For the calibration, we set Remark 3.6. In the non-fractional case β = 1, the discretized problem (3.3) is also non-fractional. Therefore, realizations of its solution can be computed directly and no quadrature is needed. If β ∈ (n, n + 1] for some n ∈ N, one may apply the above described method and error estimates to L := L n+1 and β := β n+1 ∈ (0, 1], since L β = L β . Yet, the finite element theory for the operator L may not be trivial.

An application and numerical examples
In the following numerical experiment we take up the SPDE from (1.2) in §1. More precisely, with the objective of generating computationally efficient approximations of Gaussian Matérn fields on the unit cube D := (0, 1) d in d = 1, 2, 3 spatial dimensions, we consider the following problem: x ∈ D, (4.1a) u(x) = 0, x ∈ ∂D, (4.1b) and study the above presented numerical method generating the approximation u Q h,k in (2.18). As already observed in Example 2.2, the exponent of the eigenvalue growth is given by α = 2/d in this case.
Furthermore, using a finite element discretization on uniform meshes with continuous, piecewise linear basis functions, Assumption 2.6 is satisfied for this problem in all three dimensions with r = s = q = 2, see Example 2.7. The condition in (2.19) of Theorem 2.10 becomes β > d/4. We emphasize that this assumption is meaningful also from the statistical point of view: on all of R d , β > d/4 corresponds to a positive smoothness parameter ν > 0 of the Matérn covariance function (1.1).
Thus, if the quadrature step size k and the finite element mesh width h are calibrated as indicated in Table 1 in §3.4, the theoretical rates of convergence for β ∈ (d/4, 1) are 2β − d/2 for the strong error and min{4β − d, 2} for the weak type error according to Theorem 2.10 and Corollary 3.4. For Problem (4.1) with κ = 0.5, we investigate the empirical convergence rates (i) of the strong error for d = 1, 2, 3 and β = 2d+n 8 , n ∈ {1, . . . , 7 − 2d}; (ii) of the weak type error for d = 1, 2, 3 and β = 2d+1 8 . In each dimension, we use a finite element method in space with continuous, piecewise linear basis functions on uniform meshes with mesh diameter h, mesh nodes x 1 , . . . , x N h and corresponding mass matrix M. We choose k = −1/(β ln h). The numbers of finite element basis functions and the corresponding numbers of quadrature nodes depending on β used in the strong error study are shown in Table 2. For (i), measuring the strong mean-square error between the exact solution u to our model problem (4.1) and the approximation u Q h,k in (2.18), we proceed as follows: First, samples of an overkill approximation of the white noise θ1,...,θ d ) are generated and evaluated on a uniform overkill mesh ofD with N d ok nodes. Here, {ξ θ } are independent standard normally distributed random variables and {e θ } are the eigenfunctions in (2.6). The approximation W ok corresponds to a truncated Karhunen-Loève expansion with N d ok terms. From this, samples of the overkill solution u ok are obtained via where {λ θ } are the eigenvalues from (2.7). For the sake of generating comparable samples of the approximation u Q h,k , we consider the same realizations of W ok and use the load vector b with entries b i = (W ok , φ i,h ) L2(D) instead of b ∼ N (0, M) from Remark 2.9, as ((W, φ i,h ) L2(D) ) N h i=1 ∼ N (0, M) and we treat W ok as the true white d = 1 d = 2 β = 7/8 weak type error In the third graph the observed strong error for β = 7/8 and d = 1, 2, 3 is displayed. The corresponding observed strong convergence rates are shown in Table 3. The final graph shows the observed average weak type errors for three experiments in d = 1, 2, 3 dimensions. For all of them the theoretical rate of convergence is 0.5 and the observed rates are shown in the legend. All errors are shown as functions of the mesh size h used in the computations in a log-log scale.
noise. The resulting approximation is denoted by u h,k . We choose N ok = 2 18 + 1 for d = 1, N ok = 2 12 + 1 for d = 2, and N ok = 5 · 2 6 + 1 for d = 3. For each value of β and in every spatial dimension, we use 50 samples of W ok to generate samples u The results are shown in the first three panels of Figure 1. The data set {(h , err )} is then used to compute the observed rate of convergence r as the least-squares solution to the linear regression ln err = c + r ln h for each combination of (d, β). Table 3. Observed (resp. theoretical) rates of convergence for the strong errors shown in Figure 1   ]| addressed in Corollary 3.4. Note that for this study no sample-wise comparison is needed and, thus, the load vector is sampled from b ∼ N (0, M) as discussed in Remark 2.9. The variance of the exact solution can be computed directly from the known eigenvalues of the differential operator as E[ u 2 L2(D) ] = θ λ −2β θ . The variance of u Q h,k is approximated via Monte Carlo integration by where the number of Monte Carlo samples is N MC = 10 3 and u (i) h,k denotes a realization of the numerical approximation u Q h,k . The observed weak type errors and the observed rates of convergence are displayed in the fourth panel of Figure 1. The theoretical rate of convergence predicted by Corollary 3.4 is 0.5 for all three cases.

Conclusion
We have considered the fractional order equation (2.8) with Gaussian white noise in a Hilbert space setting. We have shown that the fractional operator L β is an isometric isomorphism between theḢ-spaces in (2.4). From this result and the mean-square regularity of the white noise with respect to theḢ-spaces, we have deduced existence and uniqueness of a solution u to (2.8) with a certain regularity.
We have proposed the approximation u Q h,k in (2.18) based on two numerical ingredients: (i) finite-dimensional subspaces V h ofḢ 1 , and (ii) the quadrature approximation (2.16) of the inverse fractional operator. The most advantageous and novel properties of the corresponding numerical scheme are (a) that only solutions to integer order (i.e., local) elliptic equations have to be computed, and (b) that it does not require the knowledge of the eigenfunctions of the differential operator L.
Our main result, Theorem 2.10, shows strong mean-square convergence of the approximation u Q h,k to the exact solution u. If the quadrature step size k and the finite element mesh width h are calibrated appropriately, see Table 1, the resulting strong convergence rate depends only on the fractional order β, the dimension d, the eigenvalue growth α of the operator L, and the approximation properties of the finite element spaces. In order to prove this result, we have partitioned the strong error in three terms: the truncation error, the error caused by the finite element discretization, and the error of the quadrature approximation. We have derived bounds for each of these error terms separately in § §3.1-3.3. By means of similar techniques, we have proven a weak type error estimate in Corollary 3.4.
Finally, in §4 we have applied the proposed numerical method to an explicit problem with relevance for spatial statistics: the solution u to (4.1) can be regarded as an approximation of a Gaussian Matérn field on the unit cube (0, 1) d . The performed numerical experiments with continuous, piecewise linear finite element basis functions in d = 1, 2, 3 spatial dimensions verify the derived theoretical strong and weak type convergence rates, see Figure 1 and Table 3.