Mitigating the Influence of the Boundary on PDE-based Covariance Operators

Gaussian random fields over infinite-dimensional Hilbert spaces require the definition of appropriate covariance operators. The use of elliptic PDE operators to construct covariance operators allows to build on fast PDE solvers for manipulations with the resulting covariance and precision operators. However, PDE operators require a choice of boundary conditions, and this choice can have a strong and usually undesired influence on the Gaussian random field. We propose two techniques that allow to ameliorate these boundary effects for large-scale problems. The first approach combines the elliptic PDE operator with a Robin boundary condition, where a varying Robin coefficient is computed from an optimization problem. The second approach normalizes the pointwise variance by rescaling the covariance operator. These approaches can be used individually or can be combined. We study properties of these approaches, and discuss their computational complexity. The performance of our approaches is studied for random fields defined over simple and complex two- and three-dimensional domains.


1.
Introduction. Gaussian random fields over functions, sometimes referred to as continuously indexed Gaussian random fields, are important in spatial statistical modeling, geostatistics and in inverse problems. They are described through a mean and a covariance operator, usually defined over a Hilbert space. Efficient manipulation of random fields is of critical importance in applications. In particular, one commonly requires the application of the covariance operator and of its inverse, the precision operator, to vectors from the function space. Additionally, computation of realizations from the distribution requires the ability to apply a square root of the covariance operator to vectors.
Constructing covariance operators from elliptic PDE operators, which has recently gained popularity [1][2][3][4][5], allows one to build on available fast PDE solvers for the required manipulations. This leads to a correspondence between domain Green's functions of PDE operators and covariance functions of the Gaussian random fields. On bounded domains, PDE operators require the definition of boundary conditions, which has implications for the resulting covariance operators. Namely, this can lead to increased/decreased correlation and pointwise variance close to the boundary, which is usually undesirable from a statistical perspective. This behavior is illustrated in figure 1 and has also been observed in [2][3][4]. In this work, we present methods to ameliorate these boundary effects. Since we aim at large-scale problems, we are interested in scalable optimal complexity algorithms that avoid dense matrix operations or matrix assembly. Our target is to find domain Green's functions that are as similar as possible to the free-space Green's functions of the precision operator, which are Matérn covariance functions. We present two methods towards achieving this objective.
The first method combines the PDE operator with a homogeneous Robin boundary condition βu + ∂u dn = 0, with a varying Robin coefficient β = β(x). This coefficient function is derived as solution to an optimization problem that aims at making the difference between the domain and the free-space Green's functions small. Our approach exploits the definition of the domain Green's function and uses the fact that explicit expressions for the free-space Green's functions are available or can easily be computed numerically. For one-dimensional domains, β can be chosen such that the effect of boundary conditions vanishes completely. For twoand three-dimensional domains, β can be chosen to minimize boundary effects in an averaged sense. The approach only requires computation of inner products and is thus feasible for large-scale problems.
The second method we propose amounts to a rescaling of the covariance operator C that is constructed from elliptic PDE operators. It can be combined with the approach discussed above. This rescaled operator has constant pointwise variance (a property that C above does not have). The idea is most easily understood in finite dimensions: For a covariance matrix Σ, with diagonal D ij := Σ ij δ ij , the rescaled matrix Σ = D − 1 2 ΣD − 1 2 is also a covariance matrix, and it has a constant unit diagonal.
Related Work. In spatial statistics, the use of covariance operators is motivated by the need for fast computations [6,7]. The connection between inverse elliptic operators and Gaussian fields was originally established in [8]. Building on this connection and results in [9], the authors of [3] show that discretizing an inverse elliptic covariance operator is valid, from a statistical perspective. This results in a (discretely-indexed) Gaussian field with a sparse precision operator due to the locality of differential operators. This sparsity allows for fast application of the precision. Fast application of the covariance operator is possible building on fast elliptic PDE solvers.
A parallel approach aiming at Bayesian inverse problems was established by Stuart [1]. Contrary to [3], the author's motivation is to develop the theory of Bayesian inference in function spaces. The advocated approach (which we follow) is that all algorithms should be presented and studied in function spaces, i.e., in infinite dimensions. Taking this approach, the author is lead to the use of "Laplacian-like" precision operators [1, Assumption 2.9], which are used to define Gaussian priors for Bayesian inverse problems. In some respect, [1] and [3,6] draw similar conclusions. They argue that using covariance operators is superior to using covariance functions, both from a theoretical as well as from a computational perspective.
The role of PDE-operator boundary conditions if the domain Ω ⊆ R d is bounded was already observed to cause variance inflation near the boundary in [3]. To avoid this effect, domain extension was proposed in [3,4,6]. As an alternative and related to one of the methods we propose in this paper, in [4] the authors propose to use a Robin boundary condition of the form βu + ∂u ∂n = 0. They conduct numerical experiments to empirically find a constant, boundary effect mitigating coefficient β 2). The magnitude of the left covariance function exceeds the gray scale used to show the covariance between the centers and the points of the domain. The discrepancy between the covariance is due to the use of Neumann boundary conditions for the differential operator.
in the Robin condition for a two-dimensional circular domain. In [10] the authors suggest sampling values on the boundary according to the correct distribution and then using these values as Dirichlet data for the domain. This approach is technical in higher dimensions and it requires assembled matrices.
Contributions. The main contributions of this work are as follows: (1) The proposed methods mitigate boundary effects arising in continuously indexed Gaussian random fields when elliptic PDE operators are used to construct covariance operators. They are computationally feasible, do not require assembled matrices, nor the extension of the computational domain. (2) We present simple and fast algorithms for the approximation of the quantities used in our methods. Once these upfront computations, which depend on the domain and the PDE operator, are available, all remaining computations (covariance and precision application, and computation of samples) are as efficient as in the original method. (3) We perform a comprehensive numerical study of the proposed methods on simple as well as complex geometries, such as the Antarctica domain from [5].
Limitations. We also remark limitations of our methods. (1) To compute the optimal Robin coefficient, an integration over the domain must be performed for each point on the boundary, where the Robin coefficient is needed. 1 However, this integration can be accelerated by realizing that in many cases the Robin coefficient varies smoothly, and hence one may use interpolation between adjacent points. Moreover, the integrands decay rapidly and thus the integration can be restricted to a small part of the domain. (2) Computation of the integrals in the Robin method can be challenging due to the singularity of the integrands. As a remedy, we discuss approximations that allow computation of these integrals at an accuracy that suffices for our purposes. (3) For the variance normalization method, we require knowledge of the pointwise variance over the domain. Fortunately, this field is often smooth and thus interpolation from a few points to the entire domain is possible. Additionally, one can leverage potential symmetries in the geometry to speed up the computation of the pointwise variance. (4) The upfront computations our methods require depend on the PDE operator used to define the covariance operator. If one uses a hierarchical method in which the PDE operator varies, the proposed approach can become computationally expensive.

Preliminaries.
Let Ω ⊆ R d , d = 1, 2, 3 be a bounded open domain with piecewise smooth boundary ∂Ω. Throughout this paper, we are concerned with Gaussian measures over spaces of functions defined on Ω. We first recall definitions of Gaussian measures and Gaussian fields-see, e.g., [11,12] for details.

Gaussian measures.
Let µ a measure on a separable Banach space X and u ∼ µ. We say µ is Gaussian if ∀ ∈ X * , if there exist m real and σ nonnegative such that (u) ∼ N (m , σ ) is Gaussian. Consider X = C(Ω), the space of continuous functions on Ω with the sup-norm, so that X ⊆ L 2 (Ω). Taking this view, one can specify a Gaussian measure µ on X by first taking a mean m ∈ X and a (linear) self-adjoint positive definite trace class covariance operator C : L 2 (Ω) → L 2 (Ω). Since samples N (m, C) are continuous for the choices of C we consider below, X has full measure and by [11,Ex. 3.39] we have a Gaussian measure on C(Ω). We still write N (m, C) for the corresponding Gaussian measure on X. If h ∈ L 2 (Ω) is discontinuous, then C −1 h is empty and h, C −1 h = C − 1 2 h = ∞, informally making the likelihood of observing h zero. Thus, the Gaussian measure gives full measure to X.

2.2.
Gaussian random fields. For our purposes, a Gaussian random field is a random function u : Ω → R such that for all finite sets {x i } n i=1 ⊆ Ω, the random vector (u(x 1 ), ..., u(x n )) T is a multivariate normal. For simplicity, here we consider a centered field, i.e., m(x) := E[u(x)] ≡ 0. The corresponding covariance function c : Ω × Ω → R is defined as c(x, y) := E[u(x)u(y)]. A covariance function can also be used as a kernel for an integral operator. The resulting operator is given by If c is positive-definite [?], C is a valid covariance operator in the sense of section 2.1. This connection motivates considering C to be an inverse elliptic operator, making c the Green's function of that operator. In this case, writing c(x, y) = (Cδ x )(y) is well-defined from a PDE perspective and we use this identity below. Now, the connection with Gaussian measures is straightforward-a Gaussian random field defines a Gaussian measure on C(Ω).

Inverse elliptic covariance operators.
On Ω, consider the elliptic differential operator (1) A := −γ∆ + α with constants γ, α > 0. The domain on which A is defined depends on the choice of boundary conditions. We will discuss different domains Dom(A) and the implied properties for covariance operators derived from A. We assume that Ω is such that A is a Laplacian-like operator in the sense of [1, Assumption 2.9] when equipped with homogeneous Dirichlet, Neumann or Robin boundary conditions. The operator A −p is a valid covariance operator for p > d/2, with samples that are s-Hölder continuous for all s < min{1, p − d/2}. The covariance function of the free-space operator has a characteristic length of 8(p − d/2 γ/α meaning that at that distance away from a source x, the covariance decays to 0.1 of its maximal value (attained at x) [3]. Specifically, A −1 is a covariance operator for d = 1 and A −2 is a covariance operator for d = 1, 2, 3 [1]. The boundary conditions of A 2 are inherited from the boundary conditions of A, which we denote by BC(·) = 0, i.e., u = A −2 f is defined as the solution of the mixed system Av = f in Ω, This implies that A −1 is a square root of A −2 . This choice of boundary conditions renders sampling from a centered Gaussian with covariance operator A −2 , N (0, A −2 ), straightforward. Namely, samples are generated as u ∼ A −1 W where W is white noise, and this can be interpreted in infinite dimensions-see [3,8].
The key property is, as hinted in section 2.2, that the Green's function, G p , of A p with appropriate boundary conditions is the covariance function of a Gaussian measure with covariance operator The covariance function G p (·, ·), however, depends strongly on the boundary condition of A, as can be seen in figure 1.

2.4.
Causes of boundary effects. The reason for these boundary effects can be understood from either PDE theory or from probability theory. To illustrate the PDE perspective, consider the covariance operator A −1 on Ω := [0, 1] with homogeneous Dirichlet boundary conditions, and x ∈ Ω, y ∈ ∂Ω. Then that G(x, y) = 0, since G(x, ·) has to satisfy the boundary condition. By continuity, Green's function is small near the boundary, even if y is only close to the boundary. So for a Gaussian field u ∼ N (0, C) and x, y ∈ Ω near the boundary, Cov(u(x), u(y)) is smaller than what it would have been without the boundary condition. To illustrate the probabilistic perspective, consider the same operator A and domain and γ = 1, only with homogeneous Neumann boundary conditions. Loosely speaking, the Green's function is the amount of time a particle spends near y, given that it started its walk at x, if it is killed at rate κ : . Then the Green's function is interpreted as the amount of time a branching particle spends at y, had it started at x and if it is killed at the same rate κ. The Neumann boundary means the particle reflects off the boundary upon hitting it. 2 Since the particle is reflected off the boundary, it spends more time near it, making G(x, y) large near the boundary. Thus, the opposite happens -the covariance is larger near the boundary than what it would be without the boundary. These boundary effects (G(x, ·) is too big or too small near the boundary) can be undesirable from a statistical modeling point of view. In the next section, we review approaches based on extending the domain, and in sections 4 and 5 we present two novel methods to mitigate these boundary effects.
3. Extending the domain. The presented problem has a seemingly appealing solution-considering an extended open domain Ω ⊃ Ω with sufficiently regular boundary ∂Ω , which is far enough from Ω that boundary effects arising from ∂Ω are negligible in Ω. In this section, we present variants of this approach and discuss challenges that arise for large-scale problems.
Recall that we are particularly interested in scalable algorithms for the application of the (discretized) covariance operator, its inverse and its square root to vectors. Before discussing concrete methods that are based on domain extension, some comments are in order. First, extending the domain may result in undesired correlations between parts of the domain. An extreme example would be a domain which consists of two disjoint subdomains. In such a case, a connected domain Ω that encompasses these subdomains inevitably introduces correlations between them. Second, creating an extended domain Ω comes at a cost, both in terms of development time and computing time. For instance, it might require to extend a given mesh for Ω to a mesh for the extended domain Ω , and to manage the increased overall number of unknowns of the problem.
Let us start with introducing some notation. We consider the covariance and the precision operators A −2 and A 2 , respectively. Here, A is an elliptic PDE operator defined over Ω , which incorporates, for instance, homogeneous Neumann or Dirichlet boundary conditions at ∂Ω . Assume we are given a discretization (e.g., based on finite elements or finite differences) for functions defined over Ω, which we extend to a discretization of functions defined over Ω . We denote the number of degrees of freedom of the discretization for functions defined over Ω by n, and assume that the corresponding unknowns are ordered such that the first n 1 unknowns correspond to points that are inside or on the boundary of Ω. The remaining n 2 = n − n 1 unknowns correspond to points in Ω c , the domain extension. This implies the following block structure of the covariance and precision matrices Σ , Q ∈ R n×n , respectively.
In this setting, Σ 11 ∈ R n1×n1 can be used as covariance matrix for unknowns corresponding to points inside Ω. Note that the matrices Σ , Q in (3) might not be available in assembled form, and we might only be able to apply them to vectors. The application of the blocks to vectors can then be computed efficiently by appropriate padding of vectors with zeros, followed by truncation. To be precise, we denote by P 1 : R n → R n1 and P 2 : R n → R n2 the (Boolean) operators that restrict vectors to their first n 1 and to their last n 2 components, respectively. The corresponding adjoint operators P * 1 : R n1 → R n and P * 2 : R n2 → R n are padding-by-zero operators. For instance, for v ∈ R n1 , we can efficiently compute Σ 11 v as P 1 Σ P * 1 v. The precision for unknowns corresponding to points in Ω is found as the Schur complement, [4] Hence, fast application of the precision Σ −1 11 to vectors requires that we can apply Q −1 22 efficiently. The ability to do this depends on the specific choice of the discretization, and we discuss some special cases next. Additionally, we discuss options for applying the square root of the covariance operator, Σ 1/2 11 , to vectors, as is required for computing sample realizations from Gaussian distributions with covariance matrix Σ 11 .
3.1. Domain extension from [4]. First, we summarize the approach proposed in [4], where the authors use finite differences (for simple geometries) or finite elements (for more complicated geometries) to discretize elliptic operators defined on Ω . They assume that the matrices (3) are available in assembled form, which allows them to apply Q −1 22 using standard solvers for positive definite sparse matrices that are available in assembled form. For computing samples, which requires a square root of Σ 11 or its inverse, they assume a factorization Q = L T L, which is reasonable as we assumed that the precision is the product of two elliptic PDE operators. Computing samples from the distribution is then carried out by defining are the first n 1 columns of L and L 2 = LP * 2 ∈ R n×n2 are the last n 2 columns. With a short calculation, one verifies thatL T 1L1 = Σ −1 11 . Thus, we can obtain samples from N (0, Σ 11 ) by solvingL T 1 u = z, where z ∼ N (0, I 1 ), with the n 1 × n 1 -identity matrix I 1 . However,L 1 might not be sparse as Q −1 22 is in general dense. Thus, having to solve a linear system with coefficient matrixL 1 might not be feasible for large n.

3.2.
Modfication of approach from [4]. A modification of the approach taken in [4] is to start by defining how samples from N (0, Σ 11 ) are generated. Namely, let u ∼ P 1 L −1 z, with z a finite element approximation to white noise. Then, a short argument yields that u ∼ N (0, Σ 11 ) and using a decomposition analogous to (4), we may recover the corresponding precision. Provided systems with L can be solved efficiently, 3 this provides a method to compute samples and to apply the covariance matrix to vectors that does not require assembled matrices. However, now the bottleneck is in the need to apply the precision which requires the block Q −1 22 that is usually not available unless one assembles the matrix Q . One possibility is to use an iterative method, such as the conjugate gradient method, to solve systems with Q 22 . However, unless an efficient preconditioner for this solve is available, this can require large numbers of iterations, as can be seen in the following section.
3.3. Domain extension and Fourier bases. So far, we have considered approaches based on local discretizations for the elliptic operator A. As an alternative, we may consider using a global basis, such as a discretization based on the (non-uniform) fast Fourier transform (FFT), which allows fast application of A and A −1 . This requires the extended domain Ω to be a box, and enables fast application of Σ and Q without requiring these matrices in assembled form. However, similar as above, we do not have access to Q −1 22 because we cannot extract and invert the submatrix Q 22 . As discussed above, one option would be to solve systems with Q 22 iteratively. However, this requires a large number of iterations, making the method inefficient in practice.
3.4. Practical aspects. Each of the approaches discussed above has limitations for large-scale problems. Sampling using the method suggested in [4] requires assembling and solving a dense system. While our modification provides a fast method to compute samples, it either requires matrix assembly or an iterative Krylov method for applying the precision operator. Using a global Fourier basis on a rectangular domain extension requires an iterative method for applying the precision operator as well. We found this to take a large number of iterations that has to be performed each time the precision is applied. The methods we proposed in the next section require some upfront computation to estimate an optimal Robin coefficient or the pointwise variance field. After this step, all computations with the covariance operator can be performed efficiently and without requiring assembled matrices.

4.
Optimal Robin boundary conditions. In this section we aim at finding Robin boundary conditions that mitigate the undesirable boundary effects shown in figure 1. We derive a coefficient in the Robin condition such that the Green's functions (which are also the covariance functions) of the domain are close to the free space Green's functions.
For x ∈ Ω, we denote the free-space Green's functions for A p , centered at x, by Φ p (x, ·), p = 1, 2. Their explicit expressions depend on the dimension d of the domain, see appendix B. For p > d/2 these free-space Green's function are known as Matérn covariance functions.
The corresponding domain Green's functions with Robin boundary condition are denoted by G p (x, ·), p = 1, 2, and they satisfy: in Ω, (5a) where δ x is the Dirac-delta function centered at x ∈ Ω, and β : ∂Ω → R ≥0 is a non-negative function defined for all boundary points y ∈ ∂Ω. Following [13], we refer to the difference between the free-space and the domain Green's functions, ·) for p = 1, 2, as the correctors. These correctors satisfy the following equations: From (7), it can be seen that if the right hand sides in (7b) and (7d) were to vanish everywhere on ∂Ω, the correctors φ x p ≡ 0 and thus Φ p (x, ·) = G p (x, ·) for p = 1, 2. If these were to vanish for all x ∈ Ω, then Φ p = G p for p = 1, 2. In the remainder of this section, we present an optimization problem for the Robin coefficient β(y), y ∈ ∂Ω that aim at making the boundary right hand sides in (7b) and (7d) small, and thus Φ p ≈ G p . 4.1. One-dimensional case. In one dimension, both A and A 2 with appropriate boundary condition are valid precision operators [1]. For A, we only have to consider the system (7a) and (7b). The one-dimensional free-space Green's function appearing in (7b), is Φ 1 (x, y) = exp(−κ|x−y|) 2κγ , with κ = α/γ. It can be verified that for β := κ, the right hand side of (7b) vanishes. Thus, G 1 = Φ 1 , i.e., the domain and the free-space Green's functions coincide.
If one considers A 2 as covariance, this choice of β does not guarantee that G 2 = Φ 2 . While for β = κ, the right hand sides in (7b) and (7c) vanish, the right hand side in the boundary condition (7d) does not. Thus, for A 2 one should choose a different value for β following the approach presented in section 4.2 below.

Higher dimensions.
We consider the precision operator A 2 and propose an optimization approach for deriving an optimal Robin coefficient. As discussed above, we would like to make β(y)Φ p (x, y) + ∂Φp(x,y) ∂n , p = 1, 2 (the right hand sides of (7b) and (7d)) vanish. For a fixed x ∈ Ω and y ∈ ∂Ω, we may choose, as a compromise, β = β(y) to be the average of the roots of these terms. Note that both terms are linear in β with positive slopes Φ p (x, y), p = 1, 2. Recall that a convex parabola attains its minimum value at the mean of its two roots. Thus, the parabola (in the variable β) attains its minimum in the average of its roots. Thus, for a fixed x ∈ Ω, y ∈ ∂Ω, the minimum of the parabola may serve as a compromise between the two competing terms. However, this compromise is made for a single x ∈ Ω. In order to take into account all x ∈ Ω, we average, recovering the following optimization problem: (8) β(y) := arg min This quadratic minimization problem can be solved easily for β, noting that the constraint β ≥ 0 can be enforced on the solution. This leads to the following expression for β(y): β(y) := max(0,β(y)), Note that the integrals occurring in (9) are finite for all dimensions d = 1, 2, 3. Computingβ(y) requires the computation of two integrals over Ω. From the explicit expressions (16) and (17) forβ(y), one can verify thatβ(y) > 0, if Ω is convex.

Practical aspects.
Numerical evaluation of the singular integrals in (9) is a challenging task. We have used two practical approaches for computing Robin coefficients in the context of finite element discretizations. The first approach approximates the fundamental solutions with piecewise constants, found by evaluating the fundamental solutions and their derivatives at element centers. This avoids singularities and the integrals in (9) for the resulting  piecewise constant functions can be computed exactly. Robin boundary coefficients computed using this approach are shown in (c) and (d) in figure 2. As the mesh is refined, the Robin coefficients converge to the numerically accurate Robin coefficient, which is obtained from adaptive quadrature [14]. Our second approach is motivated by the derivation of β as presented in section 4, but for the discretized problem. We consider discrete approximations Φ h 1 and Φ h 2 of the free-space Green's functions Φ 1 and Φ 2 , and aim at solving the optimization problem (8) with these discrete Green's functions rather than their continuous counterparts. Our motivation is that Φ 1 and Φ 2 cannot be represented in finite dimensions and thus the discrete domain Green's functions can never be good approximations to the continuous free-space Green's functions. The best we can hope for is that the numerically computed domain Green's functions approximate discrete free-space Green's functions Φ h 1 and Φ h 2 . Unless for uniform discretizations, Φ h 1 (x, ·) and Φ h 2 (x, ·) depend on the discretization mesh in a neighborhood of x and thus would have to be computed for every x. To avoid this, and using the radial symmetry of Green's functions, we compute a one-dimensional finite element approximation to the free-space Green's function as illustrated next for d = 2-an analogous approach can be taken for d = 3. Recall that for a radially symmetric function v, we can use polar coordinates (r, φ) for the Laplacian operator to write: Hence, using a Dirac-delta δ 0 , we find the weak form: where the last equality follows from integration by parts and radial symmetry. Now, we solve for Φ 1 as a function of the radius r using the finite element method in one dimension. The space discretization length scale h and the polynomial order for this one-dimensional finite element calculation should be representative of their higher-dimensional counterparts, such that the resulting discrete free-space Green's functions can be used to compute the optimal Robin coefficient functions for the discrete problem. We truncate the integration over the radius to [0, R], with R sufficiently large such that the Neumann boundary condition imposed at r = R has negligible effect. To compute Φ 2 as a function of r, we solve the discretized problem twice. The usual finite element quadrature can now be used for computing the Robin coefficients, since the numerically computed free-space Green's functions are finite element functions (or interpolations of radially symmetric one-dimensional finite element functions to a two/three-dimensional mesh). The results are shown in (a) and (b) in figure 2. Moreover, Robin coefficients computed with these discrete free-space Green's functions are (close-to) optimal for a discrete version of the optimization problem (8), which is particularly relevant for coarser discretizations, i.e., in the pre-asymptotic regime.

5.
Normalizing the variance. The approach presented in this section can be used to mitigate boundary effects in covariance operators derived from elliptic PDEs with Neumann or Robin boundary conditions. In particular, it can be used in combination with the optimal Robin coefficient approach introduced in the previous section. Recall, that the correlation between two (real valued) random variables X, Y is defined as Var(X)Var(Y ) .
Now, let us consider a Gaussian random field, u, which is defined over Ω and has the covariance function G 2 with Robin or Neumann boundary conditions, and a Gaussian random field, v, defined over R d with covariance function Φ 2 . Then, A key property of v is that where σ 2 is a constant given explicitly in (14). This means that the covariance of the field v coincides with its correlation (up to a multiplicative constant). This is a desirable property from a modeling point of view, as it means that v(x) and v(y) vary on the same scale. Said differently, it is as likely to observe v(x) at a certain distance from its mean E[v(x)] as it is to observe v(y) at the same distance from its mean E[v(y)]. This is not the case, however, for u. A property similar to (10) does not hold for Var(u(x)) = G 2 (x, x). This, as will be seen in the numerical simulations, is a significant part of the boundary effect illustrated in figure 1. The idea of the approach proposed in this section is to modify the covariance operator A −2 so that its Green's functions satisfy (10) (with the constant σ 2 ).
Before presenting our method in function space, we consider its simpler analogue in R n . Consider a (symmetric positive definite) covariance matrix Σ ∈ R n×n with non-constant diagonal and define a diagonal matrix D by D ii = σ −1 Σ ii , with σ > 0. Let Λ := D − 1 2 ΣD − 1 2 and v ∼ N (0, Λ). Then, (10) holds for v in the sense that The covariance operator modification presented below is the infinite-dimensional analogue to the computation of Λ.
Consider A as in section 2, equipped with a homogeneous Robin boundary condition βu + ∂u ∂n = 0 with β : ∂Ω → R ≥0 bounded. Note that this includes a homogeneous Neumann boundary condition for β ≡ 0. We define g(x) := σ/ G 2 (x, x), the infinite-dimensional analogue of the matrix D − 1 2 defined above. Note that G 2 (x, x) is the pointwise variance field of N (0, A −2 ) and σ 2 is the pointwise variance of the free-space covariance function defined in (14) in the appendix. In the appendix (proposition 1) we show that g is bounded away from zero and infinity and that it is differentiable. Thus, C := gA −2 g is a valid covariance operator and has constant pointwise variance σ 2 (proposition 2). Also, u ∼ N (0, C) are characterized by u ∼ gv, where v ∼ N (0, A −2 ), which is a consequence of [12,Proposition 1.18].
Note that this transformation can be interpreted probabilistically using particles that follow a Brownian motion. For simplicity, we assume d = 1 such that A −1 with Neumann boundary conditions is a valid covariance operator. Then, the time a particle starting at x ∈ Ω spends in a set A ⊂ Ω before being killed (killing occurs independently at a rate κ 2 = α/γ) is A G 1 (x, y) dy. Multiplying A by g changes both the Laplacian part of the operator (responsible for Brownian motion) and the κ 2 (responsible for killing of particles). Multiplying the Laplacian by g corresponds to a time change. This does not change the distribution of Brownian paths (without killing), but changes the particle velocities along the paths. If one only changes the traveling speed, one changes the measure on paths, because the rate of killing stays the same. If the killing rate is changed by the same factor, one obtains the same distribution of paths but particles are sped up where the pointwise variance was too large and slowed down where it was too small. 5.1. Practical aspects. Note that this method requires knowledge of the pointwise variance of the covariance operator A −2 with some choice of boundary conditions. This can be an expensive computation, but there are several options to approximate the pointwise variance field.
One option is to calculate the pointwise variance through samples. For the finite element method, this involves applying a square root of the mass matrix M to vectors [2]. Since this can be a difficult task, we suggest an alternative. Denote by K the symmetric finite element discretization of A. Then, the covariance matrix is K −1 M K −1 [2], and pointwise variances of the finite element function are known to be the diagonal entries of the covariance matrix. If we set Z ∼ N (0, I) and let Thus, we may estimate the pointwise variance as follows. Draw {Z k } N k=1 iid as above, set Additionally, often symmetry properties of the domain Ω can be used to speed up the computation of the pointwise variance (as, e.g., in [2]), or an approximation for the pointwise variance field, which is typically smooth, can be obtained through interpolation with a small number of points.
The problem of estimating the diagonal of a matrix inverse has been studied extensively in the literature. For fast estimation methods for diagonals of Green's functions we refer to [15,16]. Alternatively, low-rank matrix approximation of the discretized covariance operator can be used to approximate the diagonal. The problem is considered for a wider class of matrices in [17,18] using stochastic estimation. A method based on applying an inverse of a sparse matrix to carefully chosen vectors is proposed in [19]. 6. Numerical Experiments. In this section, we study the ability of the methods proposed in sections 4 and 5 to mitigate boundary effects in two and threedimensional numerical examples. For comparison, we also present results obtained with homogeneous Neumann boundary conditions as used in [2,5], with homogeneous Dirichlet conditions, and with the constant Robin coefficient as suggested in [4]. We use the finite element library FEniCS [20] for our numerical tests, 4 and rely on linear finite elements in our computations. Unless otherwise specified, we compute the Robin boundary coefficient using the numerically computed approximate L 2 -projection of Green's functions discussed in section 4.3. For the parallelogram and Antarctica meshes we calculate the pointwise variance at every discretization point x ∈ Ω directly as (A −2 δ x )(x). For the cube mesh we do so using our stochastic estimator derived in section 5.1 with 10, 000 samples, which we find to result in reasonable approximations.  6.1. Parallelogram example. We first illustrate our methods on a two-dimensional domain that is more challenging than the square domain in figure 1. The results shown in figure 3 show cross sections through covariance functions centered at a particularly challenging point close to a corner of the domain. We use A = −∆ + 121 as the square root of the precision operator. We discretized the unit square by 128 2 points and then transformed it to the parallelogram using a linear transformation such that its vertices become (0, 0), (cos θ − , sin θ − ), (cos θ + , sin θ + ) and (cos θ − + cos θ + , sin θ − + sin θ + ), where θ ± = π 4 ± π 8 . As can be seen in figure 3, using Robin boundary conditions results in a significant improvement over Neumann boundary conditions. In this problem, the constant Robin coefficient β = √ α/1.42 from [4] and the variable Robin coefficient perform similarly. Moreover, figure 3 also shows that the variance normalization method results in constant pointwise covariance variances, but that the resulting covariance functions differ from the free-space covariance functions. Combining varying Robin boundary conditions with variance normalization leads to the best results.
6.2. Antarctica domain example. We also show Green's functions on the Antarctica domain used for Bayesian inference in [5]. We use A = −∆+α, α = 10 −5 as the square root of the precision operator, and measure distances in kilometres (Antarctica extends laterally between 3000 and 6000 kilometers). Note that in [5], the authors used α = 10 −6 , which leads to stronger point correlation. 5 We used a finite element discretization with 27,749 cells. Figure 4 shows a comparison of two domain Green's functions, one centered far and one close to the boundary. The differences between the covariance functions on the left of the domain (which is West Antarctica) is due to the different boundary conditions. As for the previous example, we find that using Robin boundary conditions largely mitigates undesired boundary effects.
We also show pointwise standard deviation (i.e., the square root of the pointwise variance) fields in figure 5. Since the free-space Green's functions are independent of the boundary, deviation from a constant pointwise standard deviation is an indicator for the strength of undesired boundary effects. We only show standard deviation fields for variants of Robin boundary conditions as the variance normalization methods discussed in section 5 ensures constant standard deviation for the resulting operator. We find that using Dirichlet or Neumann boundary conditions can have a significant effect also on the pointwise standard deviations fields. These boundary effects are significantly diminished for the cases with Robin boundary conditions.   If there was no positive lower bound for G 2 (x, x), then due to the compactness ofΩ, there was x ∈Ω such that G 2 (x, x) = 0 (we extend Green's functions to ∂Ω). However, (11) implies that G 1 (x, ·) = 0 almost everywhere. From the probabilistic interpretation of these Green's functions as (density of) time spent at a point this can only happen if a particle is immediately killed at x. This is not possible for an interior point and can only be possible for a boundary point with a homogeneous Dirichlet boundary condition, which we exclude. We may conclude that no such sequence exists and that there is some lower bound c > 0 for which G 2 (x, x) > c > 0.
We know G 2 (x, x) is the pointwise variance of a u ∼ N (0, A −2 ) and by the Karhunen-Loève expansion u(x) = k∈K λ 1/2 k φ k (x)ξ k with ξ k ∼ N (0, 1) iid and {λ k , φ k } k∈K eigenpairs of the (trace-class) operator A −2 . Then Since A is a Laplacian-like operator according to [1, Assumption 2.9] we have a uniform bound on φ k ∞ and since A −2 is trace-class k∈K λ k < ∞. Using these facts in (12) gives a uniform upper bound G 2 (x, x) < C < ∞, as desired.
Definition A.1. Let C be defined via [Cu](x) := g(x)[A −2 (gu)](x) with g(x) as in proposition 1. In the following, we write C = gA −2 g, with the understanding that A −2 operates on the product of all functions to its right.
We now show that C is a valid covariance operator with constant pointwise variance.
Proof. First observe that: Since its diagonal is constant, C is trace class with trace T r(C) = E u∼N (0,C) u 2 2 = σ 2 |Ω| (the first equality follows from the Karhunen-Loève expansion). Let u ∈ L 2 (Ω). By proposition 1, ug ∈ L 2 (Ω). Then positive definiteness follows from the fact that A −2 is positive definite. A straightforward calculation shows that C is self-adjoint in the L 2 (Ω) inner product. Using g > 0 from proposition 1 and the fact that A is invertible on Dom(A), it is easy to verify that The last statement follows from [12, Proposition 1.18].
Appendix B. Explicit expressions for varying Robin coefficients. Here, we give explicit expressions for the Robin coefficient function β from section 4. Let us first recall the expressions for the free-space Green's functions for A and A 2 . For an elliptic differential operator L, the free-space Green's function is defined (informally) as the solution to the equation where L operates in y.