{'fraction_non_alphanumeric': 0.04259359234506767, 'fraction_numerical': 0.01493009661745816, 'mean_word_length': 4.068692206076618, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A probabilistic finite element method based on random meshes: Error estimators and Bayesian inverse problems

We present a novel probabilistic finite element method (FEM) for the solution and uncertainty quantification of elliptic partial differential equations based on random meshes, which we call random mesh FEM (RM-FEM). Our methodology allows to introduce a probability measure on standard piecewise linear FEM. We present a posteriori error estimators based uniquely on probabilistic information. A series of numerical experiments illustrates the potential of the RM-FEM for error estimation and validates our analysis. We furthermore demonstrate how employing the RM-FEM enhances the quality of the solution of Bayesian inverse problems, thus allowing a better quantification of numerical errors in pipelines of computations.


Introduction
In recent years, there has been a growing interest in developing and analyzing probabilistic counterparts of traditional numerical methods spanning most areas of computational mathematics. This gave rise to the field of Probabilistic Numerics (PN), whose founding principles and aims are summarized in the review papers [15,23,31]. All methods belonging to the field of PN share the idea of introducing a probability measure on the solution of traditional numerical methods. The underlying rationale is to quantify the uncertainty due to numerical errors in a probabilistic manner, rather than with standard error estimates. Indeed, a probability measure over approximate solutions can be readily pushed through a pipeline of computations, thus justifying the need of probabilistic methods especially when the solution of the problem at hand is employed as the input of a subsequent analysis. A typical example of computational pipelines for which probabilistic methods are successfully employed is given by Bayesian inverse problems, where introducing a probability measure on the forward model allows for a better quantification of the uncertainty in the inversion procedure.

Literature Review
Several contributions to the field of PN concern differential equations. For ordinary differential equations (ODEs), the methodologies can be roughly split in two different areas. In [11,25,26,[38][39][40]45] the authors present a series of schemes which rely in different measure on Bayesian filtering techniques. These methodologies proceed by updating Gaussian measures over the numerical solution with filtering formulae and evaluations of the right-hand side of the ODE, which are interpreted as observations. While being not involved computationally, analyzing the convergence properties of this class of methods is not always possible, and one can only marginally rely on standard techniques for this purpose. A valuable effort in this sense can be found in [26], where the authors show rates of convergence of the mean of the Gaussian measure towards the exact solution. A different approach is presented in the series of works [2,16,28,43,44], where the authors propose probabilistic schemes which are based on perturbing randomly the approximate solution and on letting evolve these perturbations through the dynamics of the ODE. In this manner, it is possible to obtain empirical probability measures over the otherwise deterministic numerical solution. A random perturbation can be applied directly to the state, as it was presented and analysed for one-step methods in [16,28], with a particular focus on implicit schemes in [43] and for multistep methods in [44]. Another approach, which was presented in [2], consists in perturbing the scheme itself by randomizing the time steps of a Runge-Kutta method. This allows to maintain certain geometric properties of the deterministic scheme in its probabilistic counterpart, such as the conservation of invariants or the symplecticity.
There has been a keen interest from the PN community on developing probabilistic numerical solvers for partial differential equations (PDEs), too [11, 13, 14, 16, 21, 30, 32-34, 36, 37]. In [13], the authors present a meshless Bayesian method for PDEs, which they then apply to inverse problems in [14], and in particular to a challenging time-dependent instance drawn from an engineering application in [30]. Their methodology consists of placing a Gaussian prior on the space of solutions, thus updating it with evaluations of the right-hand side, which are interpreted as noisy observations. A similar idea has been presented in [11], where the main focus are time-dependent problems, and in [36,37], where the method is recast in the framework of machine learning algorithms. In [33,34], a probabilistic approach involving gamblets is applied to the solution of PDEs with rough coefficients and by multigrid schemes, with a particular interest to reducing the complexity of implicit algorithms for time-dependent problems [34]. Moreover, in [32] the author presents a Bayesian reinterpretation of the theory of homogenization for PDEs, which can be seen as a contribution to the field of PN. To our knowledge, the only perturbation-based finite element (FE) probabilistic scheme for PDEs is presented in [16], where the authors randomize FE bases by adding random fields endowed with appropriate boundary conditions, thus obtaining an empirical measure over the space of solutions. By tuning the covariance of these random fields, they obtain a consistent characterization of the numerical error, which can then be employed to solve Bayesian inverse problems and to quantify the uncertainty over their numerical solution.

Our Contributions
In this work we present a probabilistic finite element method (FEM) which is based on a randomization of the mesh, and which we call RM-FEM. The idea underlying our method stems from both [16], where the authors propose a probabilistic FEM based on random perturbations, and from [2], in which the first instance of randomizing the discretization instead of the solution itself is presented. In the context of ODEs, a careful randomization of the time step in Runge-Kutta methods allows to maintain certain convergence and geometric properties, either path-wise or in the mean-square sense. In a similar fashion, creating a probability measure on the space of solutions by randomizing the mesh has the advantage that each sample is a FEM solution itself, and therefore a projection of the exact solution on some random finite-dimensional space.
Keeping in mind the fundamental goal of PN, we consider the problem of employing probabilistic methods to quantify numerical errors in the context of PDEs. Indeed, in [2,16,26] and other works concerning ODEs, the authors show that the probabilistic solution converges to the true solution with the same rate as the deterministic method, which represents a consistency result. No work so far shows that PN methods can be readily employed for an a posteriori estimation of the error. Some forms of adaptivity for nonlinear ODEs based on probabilistic information can be found in [6,10,39], where the arguments are based on heuristics but are not rigorously analyzed. In this work, we construct and present a posteriori error estimators which can be readily employed for mesh adaptation in elliptic PDEs. Our estimators are entirely based on probabilistic information, are simple to compute and do not entail considerable computational cost. We present an analysis in the one-dimensional case that shows that our error estimators based on the RM-FEM are equivalent to a classical estimator by Babuška and Rheinboldt [5], which employs the jumps of the derivative of the solution at the nodes to quantify the numerical errors. Our one-dimensional theoretical analysis is complemented by a series of numerical experiments confirming the validity of our theory in higher dimensions.
As stated above, probabilistic numerical methods are especially appealing when employed in pipeline of computations such as Bayesian inverse problems. In particular, employing deterministic methods for approximating forward maps leads to overly confident posterior measures, which can be corrected by appropriate probabilistic approximations. Similarly to [2,16], we show in this paper how the RM-FEM can be employed to construct empirical distributions over the forward problem and compute a random posterior measure, solution to the inverse problem in the Bayesian sense. The solution is consistent asymptotically with respect to the mesh spacing, but its quality is enhanced if the latter is relatively large, i.e., if the forward model is approximated cheaply.

Outline
The outline of the paper is as follows. In Section 2 we state the problem of interest, introduce the RM-FEM and the main assumptions and notation required by our analysis. We then present the two main applications of the RM-FEM, i.e., a posteriori error estimators and Bayesian inverse problems, in Sections 3 and 4, respectively. For both applications, a series of numerical experiments in the one and two-dimensional cases illustrate the usefulness and efficiency of the RM-FEM. In Section 5 we present a rigorous a priori and a posteriori error analysis. Finally, in Section 6 we draw our conclusions.

Notation
Let d = 1, 2, 3 and D ⊂ R d be an open bounded domain with sufficiently smooth boundary ∂D.
For v ∈ R d , we denote by v 2 the Euclidean norm on R d . We denote by L 2 (D) the space of square integrable functions, by (·, ·) the natural L 2 (D) inner product, and by H p (D) the Sobolev space of functions with p weak derivatives in L 2 (D). Moreover, we denote by H 1 0 (D) the space of functions in H 1 (D) vanishing on ∂D in the sense of traces, by H −1 (D) the dual of H 1 0 (D) and by ·, · the natural pairing between H −1 (D) and H 1 0 (D). We equip the space H 1 0 (D) with the norm v H 1 0 (D) = ∇v L 2 (D) , i.e. the H 1 (D) seminorm. For an event space Ω, with a σ-algebra A and a probability measure P , we let the triple (Ω, A, P ) denote a probability space. For an event A ∈ A, we say that A occurs almost surely (a.s.) if P (A) = 1. For n ∈ N we call random variables the measurable functions X : Ω → R n , and denote by L 2 (Ω) the space of square integrable random variables, with associated inner product. Denoting by B(R n ) the Borel σ-algebra on R n , we say that a probability measure µ X on the measurable space (R n , B(R n )) satisfying µ X (B) = P (X −1 (B)) for all B ∈ B(R n ) is the measure induced by X, or equivalently the distribution of X. For a set of random variables {X i } n i=1 which are independent and identically distributed, we say they are i.i.d., and denoting by µ their common induced measure on (R n , B(R n )), we write

Problem and Method Presentation
Let κ ∈ L ∞ (D, R d×d ), f ∈ H −1 (D) and u be the weak solution of the partial differential equation i.e., the function u ∈ V ≡ H 1 0 (D) satisfying for all functions v ∈ V . We assume there exist positive constants κ andκ such that for all ξ ∈ R d κ ξ where · 2 is the Euclidean norm on R d , so that there exist constants m, M > 0 such that for all The Lax-Milgram theorem then guarantees that the problem (2.2) is well-posed.
Let N be a positive integer and let T h = N i=1 K i be a partition of D, where for all i = 1, . . . , N , the element K i ⊂ D is a segment, triangle or tetrahedron for d = 1, 2, 3 respectively. We denote by h i = diam(K i ) the radius of the smallest ball containing K i , and by h = max i h i the maximum radius, indexing the mesh T h . We denote by V h the set of all vertices of the elements of T h , and in particular as V I h ⊂ V h the set of vertices which do not lie on the boundary of D, and by Moreover, we denote by N I the number of internal vertices, i.e., N I = V I h . We assume the partition to be conforming, i.e., if two elements have non-empty intersection, than the latter consists of a point (for d = 1), of either a vertex or a side (for d = 2), and of either a vertex, a segment or a face (for d = 3). We then denote by V h ⊂ V , dim(V h ) < ∞ the space of continuous piecewise linear finite elements on T h , i.e., where P 1 is the space of linear functions. Let us remark that imposing u h = 0 on ∂D yields dim(V h ) = N I . The FEM proceeds by finding u h ∈ V h such that for all v h ∈ V h , which is equivalent to solving the linear system Au = f , where and where {ϕ j } N −1 j=1 are the Legendre basis functions defined on the internal vertices of T h . The assumptions on κ guarantee that A is symmetric positive definite, and in turn that u is uniquely defined and the problem (2.3) is well-posed.
We now introduce the random-mesh finite element method (RM-FEM), which is based on a random perturbation of the mesh T h obtained by moving the internal vertices. First, we here detail how we build perturbed meshes and which kind of random perturbations we consider to be admissible.
be a sequence of random variables and let us define the set of internal points We then define the set of perturbed vertices as , the vertices on the boundary are left unchanged. The perturbed mesh is then simply T h = N i=1 K i , where each element K i has the same vertices as its corresponding element K i in the original mesh, modulo the random perturbation (2.4). In other words, we compute the internal points of the perturbed mesh following (2.4), and keep the connectivity structure of the original mesh T h . Clearly, the mesh so defined is not conforming for any sequence of random variable α, for which we therefore introduce an assumption.
Assumption 2.1. The sequence of random variables α is such that (i) its components α i admit densities F αi with respect to the Lebesgue measure on R d , which satisfy supp(F αi ) ⊂ B ri , where B ri ⊂ R d is the ball centered in the origin and of radius r i > 0, and which are radial, i.e., F αi (x) = F αi ( x 2 ), (ii) the perturbed mesh T h is conforming a.s.
Let us remark that the assumption (i) actually implies for all p ≥ 1 the assumption (ii) a.s., provided the radii r i are chosen small enough. We assume in (i) the densities F αi to be radial functions so that the random perturbations do not have a privileged direction.
we callh i the minimum element size for the two intervals sharing the point x i as a vertex, i.e.,h i := min{h i , h i+1 }. Then, a choice of random variables satisfying Assumption 2.1 is given by where for a set D ∈ R d we denote by U(D) the uniform distribution over D. With this choice, indeed, we have that x i < x i−1 a.s., and therefore the perturbed mesh is conforming. In the two-dimensional case, we introduce for i = 1, . . . , N I the notation Analogously to the one-dimensional case, we writeh i := min j:Kj ∈∆i h j . In this case, it is possible to verify that choosing for all i = 1, . . . , N I then α satisfies Assumption 2.1. We verify this graphically in Fig. 1, where we show a realization of the perturbed mesh based on a generic Delaunay mesh and on a structured mesh on D = (0, 1) 2 along with the sets where the perturbed points are constrained to belong a.s. We notice that for p > 1 the magnitude of the perturbations clearly tends to vanish. Finally, we remark that similar admissible perturbations can be introduced in higher dimensions.
Having defined the perturbed mesh, we now proceed with describing the RM-FEM. Let V h be the space of continuous piecewise linear finite elements on T h . Let moreover { ϕ i } N I i=1 be the Legendre basis functions defined on the internal vertices of T h and I : C 0 (D) ∩ V → V h be the Lagrange interpolation operator onto V h , i.e., for a function v ∈ C 0 (D) ∩ V and for x ∈ D we define We are then interested in the two functions belonging to the finite element space V h whose definition we give below.
Remark 2.5. Clearly, either for any fixed p ≥ 1 and h → 0 or for any fixed h < 1 and p → ∞, the functions u h , Iu h and u h tend to coincide. We visualize this for u h and u h in Fig. 2, where we simply fix κ = 1 and the right-hand side f such that u = sin(2πx) in (2.1), and consider the effects of increasing p and decreasing h. For this simple problem, we notice that for p = 2 and N = 20 the FEM solution u h and the RM-FEM solution u h are almost indistinguishable. Remark 2.6. All the quantities distinguished by a tilde (e.g., T h , V h , I) are random variables with values in appropriate spaces. For example u h is a random function u h : Ω × D → R, such that Ω × D (ω, x) → u h (ω, x). For economy of notation, in the following we drop the argument ω from all random variables.
Remark 2.7. The coefficient p in (2.4) has the same role as the coefficient identified by the same symbol in both [2,16], i.e., it controls the variability of the probabilistic solutions by tuning the variability of the noise which is applied to the method. Before proceeding with the two main applications of the RM-FEM, i.e., a posteriori error estimators and Bayesian inverse problems, we state an a priori error estimate, which suggests how to balance the sources of error due to numerical discretization and to the randomization of the method, respectively. Theorem 2.9. Let the solution u of (2.1) be u ∈ H 2 (D). Then, it holds for a constant C > 0 independent of h. Moreover, if p = 1 in (2.4) the numerical and discretization errors are balanced with respect to h, i.e., it holds This results indicates that one should fix p = 1 in (2.4) in order to obtain a family of probabilistic solutions whose statistical properties should reflect the true error. This is crucial when the RM-FEM is employed in a pipeline of computations such as Bayesian inverse problems, which will be presented in detail in Section 4. The proof of Theorem 2.9 is elementary and discussed in Section 5.1.

A Posteriori Error Estimators based on the RM-FEM
The first and foremost application of the RM-FEM is deriving a posteriori error estimators which are entirely based on the statistical information carried on by the mesh perturbation. We say that a quantity E h is an a posteriori error estimator if it gives an error estimate on the numerical approximation and is computable only by knowledge of the numerical solution. Moreover, if there exist constants C up and C low independent of h and of u such that we say that the a posteriori error estimator is reliable and efficient, respectively. Indeed, the upper bound above guarantees that when the estimator is small, so is the numerical error. The lower bound, instead, gives an insurance on the quality of the estimator, as it shows that the estimation of the error is not exceedingly pessimistic. There exist in the literature a huge number of a posteriori error estimators, and we refer the reader to the surveys given e.g. in [4,47]. Most a posteriori error estimators are expressed in the form where the η K are local quantities depending on the solution and the data on the element K and its neighbors. For example, in the two-dimensional case a valid a posteriori error estimator is given by the expression of its local components where · is the jump operator and ν K denotes the unitary vector normal to the boundary of K (see e.g. [46,Section 3] or [4, Chapter 2]). Other a posteriori error estimators are based on recovered gradients, which are employed as surrogates of the gradient of the exact solution to estimate the error. A notable member of these methodologies is the Zienkiewicz-Zhu (ZZ) patch recovery technique [49,50], which is proved to be superconvergent on special meshes, and which is in practice widely employed on any mesh.
It has been heuristically noted for ODEs in [6,10,39] that information on the variability of a probabilistic solution can be employed to estimate the error and thus adapt the numerical discretization. Indeed, building probabilistic solution to otherwise deterministic problems should pursue the goal of quantifying numerical errors through uncertainty. Guided by this observation, we now introduce two probabilistic error estimators for elliptic PDEs.

Definition 3.1. Let Iu h be the RM-FEM interpolant defined in Definition 2.3 and for each
K ∈ T h , let us denote by K ∈ T h its corresponding element in T h . We define the first RM-FEM a posteriori error estimator as .
Moreover, we define the second RM-FEM a posteriori error estimator as in the definition of η K,1 and η K,2 are necessary to obtain well-calibrated error estimators. This is made clearer in the one-dimensional case by the analysis presented in Section 5.2. For higher dimensions, they can be partially explained with the ansatz (5.2), especially for the first estimator E h,1 , and they appear in practice to be the correct scaling. Remark 3.3. Computing the estimator E h,1 is more involved than the estimator E h,2 . Indeed, for the latter it is sufficient to compute the interpolant Iu h and the gradients over each element of u h and of the interpolant. For E h,1 , instead, one has to compute on each element K the quantity By construction, each element K overlaps with the elements corresponding to its neighbors in the original mesh in a non-trivial manner, and if d > 1 one has to rely to the construction of a "supermesh" (see e.g. [18,19]) such that on each of its elements the quantity ∇(u h − Iu h ) is constant. A super-mesh has to be built for each realization of the perturbed mesh T h , which could therefore be expensive.
In this article, we show in the one-dimensional case that the estimators given in Definition 3.1 are reliable and efficient in the sense of (3.1). In the statement of our theoretical result, which is given below, we make use of a quantity Λ ∈ R which is of higher order in most practical scenarios and which is defined as where for each K j the real constant C j will be specified in the analysis of Section 5.2 (see e.g. [5, Equation (8.7)]). Moreover, we consider one-dimensional meshes which are λ-quasi-uniform, i.e., we assume there exists a constant λ ∈ (1, ∞) such that it holds uniformly in h. Finally, we consider perturbations satisfying such that |ᾱ 1 | ≤ 1/2 a.s. These perturbations are indeed the same as the ones presented in Example 2.2, but without the assumption of {ᾱ i } N i=1 to be uniformly distributed, which is not necessary in the following. In practice, a uniform distribution is nevertheless advisable, as it is still general enough and satisfies the radial assumption of Assumption 2.1(i). We moreover introduce the following technical assumption on the perturbation. Assumption 3.4. Let the family of meshes T h be λ-quasi-uniform, let p be the coefficient introduced in (2.4) and assume that for all h and p there exists C > 0 such that Remark 3.5. We note that Assumption 3.4 holds for p > 1 and h sufficiently small, and is therefore not restrictive in practice.
We can now state the main result involving a posteriori error estimators.
, up to higher order terms in h and under Assumption 3.4 for k = 1. If additionally κ ∈ C 2 (D) and f ∈ C 1 (D), then there exist constants C low and C up independent of h and of the solution u such up to higher order terms in h and under Assumption 3.4 for k = 1.
Let us notice that the estimators given in Definition 3.1 involve the computation of an expectation with respect to the random perturbations of the mesh, and therefore a Monte Carlo simulation is realizations of the estimator E h,k , obtained with independent perturbations of the mesh. Then, in practice we compute Remark 3.7. It could be suggested that the application of Monte Carlo techniques increases dramatically the simulation time. We argue that in practice the computational overhead is not relevant, mainly for three reasons. First, it has been proved in [2] that the variance of Monte Carlo estimators drawn from probabilistic numerical methods decreases with respect to the discretization size h. Hence, the number of simulations M does not need to be large, nor increasing if h → 0, to guarantee a good quality of the estimator. The same arguments hold for the RM-FEM, too. Second, the Monte Carlo estimation is completely parallelizable, thus reducing the cost by a factor equal to the number of available computing units. Finally, the computation of the RM-FEM interpolant Iu h is not computationally involved, neither when it is repeated N MC times.

Numerical Experiments
We now present numerical experiments on one and two-dimensional test cases to demonstrate the validity of our a posteriori error estimators. In particular, we are interested in determining whether the probabilistic error estimators introduced in Definition 3.1 are indeed reliable estimators for the numerical error in the FEM, and in employing these estimators for local refinements of the mesh. Setting a tolerance γ > 0, our goal is building a mesh T h such that Indeed, in this case and thus (3.4) holds. Let us remark that C 2 up is not known a priori in practice, and therefore we just decide to employ the condition (3.5) fixing C up = 1 in our experiments. We therefore adapt the mesh by computing the local contributions and comparing them with γ loc , thus locally refining the mesh if the condition (3.5) is not met, and coarsening if the local estimators are excessively small with respect to γ loc .
In the following we employ for both the one and the two-dimensional cases the uniform distributions given in Example 2.2 for the random perturbations of the points. In light of Lemma 5.3 and Lemma 5.4, we decide to correct the estimators by normalizing them with respect to the random perturbations. In particular, in the following, the estimators are normalized as

Figure 4:
Exact solution u 1 for the experiment of Section 3.1.2. Both the contour and the three-dimensional view highlight the steep gradient that features u 1 .

One-Dimensional Case
We first consider d = 1 and the two-point boundary value problem (5.4) with κ and the exact solution u given by where we fix a = 15 and b = 50, and where we choose the right-hand side f so that u is indeed the solution. As a goal, we set the tolerance γ = 10 −2 in (3.4) and stop the algorithm when condition (3.5) is met by all elements of the mesh. We consider the RM-FEM implemented with uniform random variables as in Example 2.2 and fix p = 3 in (2.4). Moreover, we consider N MC = 20 realizations of the probabilistic mesh to approximate the error estimator as in (3.3). We then compute both the error estimators given in Definition 3.1 and employ E h,1 for adapting the mesh by refinement and coarsening, guided by the condition (3.5). The adaptivity algorithm is initialized with a mesh T h built on N = 30 elements of equal size and proceeds by refinement and coarsening. Results, given in Fig. 3, confirm the validity of our probabilistic error estimators. In particular, we remark that the local error estimators succeed in identifying the regions where the mesh has to be refined, thus getting a solution with an approximately equal distribution of the error over the domain. Both probabilistic estimators, moreover, succeed in bounding the global error until the tolerance is reached, with the estimator E h,2 which appears to be more efficient than E h,1 .

Two-Dimensional Case
We now present two numerical experiments conducted in the two-dimensional case. For the adaptivity algorithm, we start from a coarse mesh and apply regular local refinements if the condition (3.5) is not met by the local error estimator η K,2 . In the two-dimensional case we do not apply coarsening to the mesh.
We first consider D = (0, 1) 2 , the conductivity κ = 1, so that (2.1) reduces to −∆u = f with homogeneous Dirichlet boundary conditions. Moreover, we choose the right-hand side f such that where β > 0. The solution has a steep transition around the line {y = 4 √ 2/5 − x}, whose steepness is proportional to the parameter β. In Fig. 4, we show the exact solution for β = 20, which we fix for this experiment. We initialize the adaptivity procedure with a mesh with maximum element size h = 1/5 and proceed with adaptation until a tolerance γ = 0.1. In Fig. 5 we show the convergence of E h,2 with respect to the convergence of the true error, as well as the the effictivity index for this experiment. We can see that the estimator indeed captures the error globally. In Fig. 6, we show the behavior of the local contributions η K,2 with respect to the true error on each element, as well as the mesh adaptation. We can see that the error estimator succeeds in identifying the region where gradients are the steepest and proposes a mesh which appears adapted to this problem.
We then consider the L-shaped domain with the re-entrant corner on the origin, i.e. D = (−1, 1) 2 \ (−1, 0) 2 . We set κ = 1, f = 0 and fix a inhomogeneous Dirichlet boundary conditions u = g on ∂D, with g chosen such that the exact solution satisfies where (r, ϑ) ∈ R + × (0, 2π] are the polar coordinates in R 2 . The exact solution of this problem is given in Fig. 7. Let us remark that the gradient of the exact solution is singular at the re-entrant corner, and we expect the mesh to be refined consequently at the singularity. For this experiment, we fix the tolerance γ = 0.03, and initialize the mesh to have a maximum element size of h = 1/3. Results, given in Fig. 8 and Fig. 9, show on the one hand that the estimator reproduces well the behavior of the global error during adaptation, and on the other hand that the mesh is progressively refined at the singularity as expected.

The RM-FEM for Bayesian Inverse Problems
Probabilistic numerical methods have been demonstrated to be particularly effective in the context of Bayesian inverse problems [2,11,13,14,16,29,30]. where Σ ∈ R m×m is a non-singular covariance matrix on R m . Given an observation y * = (u(x * 1 ), . . . , u(x * m )) associated to an unknown value ϑ * ∈ X and corrupted by observational noise β ∈ R m the inverse problem can then be stated as: find ϑ * ∈ X given observations y * = G(ϑ * ) + β. (4. 2) The randomness and the mismatch between the dimensionality of the unknown and of the observation make problem (4.2) ill-posed. Regularization can be achieved in the Bayesian framework (see e.g. [20,41]) by introducing probability measures on the unknown, which summarize prior knowledge and the information provided by data. We briefly introduce the Bayesian paradigm in the remainder of this section. First, we restrict ourselves to the space H = C 0 (D) ∩ V , which is a valid subspace of admissible values for ϑ, i.e., H ⊂ X. We then introduce a prior measure µ 0 on H, encoding all knowledge on the unknown ϑ before observations are obtained. In particular, we consider a Gaussian prior measure µ 0 = N (0, Γ 0 ) on H, where Γ 0 is a positive semi-definite covariance operator on H, and such that µ 0 (H) = 1, so that any sample from µ 0 is in H, a.s. We set the mean of the prior measure to zero without loss of generality. A broader class of prior measures could be employed, such as Besov or heavy-tailed measures [20,42], but we restrict ourselves to the Gaussian case for simplicity. Finally, we can obtain a measure µ y on X encoding all the knowledge on ϑ given the prior and the observations. We call this measure the posterior, and formally compute with Bayes' formula its Radon-Nykodim derivative with respect to the prior as dµ y dµ 0 (ϑ) = 1 Z y exp(−Φ(ϑ; y)), and where for any y ∈ R m the potential Φ(·; y) : X → R is given due to the Gaussian assumption on the noise β ∼ N (0, Σ) by For economy of notation, in the following we drop the dependence of µ y and Z y on the data, and simply denote these quantities µ and Z. In order for (4.3) to be well-defined, the posterior measure needs to be absolutely continuous with respect to the prior. This is ensured under some conditions on the forward map G, in particular, Lipschitz continuity and some controlled growth condition. Precisely, G can be shown to satisfy [ where the potential Φ h is given by  and where the normalization constant Z h is defined equivalently to (4.4). A natural question arising from this setting is whether the approximate posterior µ h converges to the true posterior µ in the limit h → 0. This is indeed true, as it holds d H (µ, µ h ) → 0 for h → 0, where the Hellinger distance d H (·, ·) is defined as For a proof of this result, see [41,Theorem 4.6], where the statement above is made more precise by the introduction of convergence rates.
It has been demonstrated heuristically that the approximate posterior measure µ h can be overly confident on the parameter if h > 0 is a finite value and if the observation model is precise, i.e., when the covariance Σ of the observational noise is small [2,11,13,14,16,29,30]. In particular, this is accentuated when h is big relatively to the forward problem, or in other words when we employ a poor approximation of the forward map G. It is therefore useful in applications to have a cheap surrogate which can be evaluated quickly, without renouncing to a complete uncertainty quantification of the solution to the inverse problem. Probabilistic numerical methods can be employed for this purpose. Let S h : X → V h be the solution operator mapping S h : ϑ → u h , where u h is the RM-FEM solution of Definition 2.4. In particular, in this context it is advisable to fix p = 1 in (2.4), so that the random deviations of the probabilistic numerical solution are of the same order of magnitude than the error itself by Theorem 2.9. Coherently to the definitions above, we define the random forward map G h : X → R m as G h = O • S h and the approximate random posterior measure µ h on X as where the potential Φ h and the normalization constant Z h are defined as above. Let us remark that the posterior µ h is a random measure, as it depends on the random variable α : Ω → R N I governing the random perturbations of the mesh. To be more precise, the posterior µ h is a random variable µ h : Ω → M(H), where M(H) is the space of probability measures over the space H.
Employing the tools of [29] and due to Theorem 2.9, it is possible to prove convergence results for µ h towards the true posterior µ for h → 0 in the Hellinger metrics.
Remark 4.1. There exist other approaches to factor the effects of discretization into Bayesian inverse problems. In particular, numerical error can be treated as modelling discrepancies. Under the assumption that the error is independent of the observational noise and of the inferred parameter ϑ, a viable alternative to probabilistic methods is employing the techniques introduced in [8,9], and further applied and analysed e.g. in [1,3]. We argue that while assuming numerical errors to be independent of the observational noise is reasonable, their independence from the inferred parameter itself is not, at least for the inverse problem (4.2).

Implementation Details
We now detail how one can solve the inverse problems above in practice.
h . Then, we define and ]. Still, the random variable ϑ is infinite-dimensional, and we need to define a finite-dimensional approximation in order to obtain a practical procedure to generate the above samples and thus solve the inverse problem. We recur to the Karhunen-Loeve expansion (KL). Denoting by {(λ i , ϕ i )} i≥1 the ordered eigenvalues/eigenfunctions of the prior covariance Γ 0 , a function ϑ ∼ µ 0 is given by the convergent sum ∼ N (0, 1). We then let N KL be a positive integer and truncate the sum above as thus obtaining a function ϑ ∈ H which is approximately sampled from µ 0 . Due to the superquadratic decay of the eigenvalues of Γ 0 = −∆ −α for α > d/2, disregarding the tail of the sum causes a negligible error in case N KL is chosen appropriately large. Our inversion problem is therefore shifted to computing the posterior distribution on a finite-dimensional parameter, comprising the coefficients of the expansion (4.8). We define the mapping K : R N KL → H, K : ξ → ϑ by (4.8).
The prior measure on ξ is µ 0 = N (0, I), with I being the identity matrix of dimension N KL × N KL , and we denote by π 0 the density of µ 0 with respect to the Lebesgue measure. The density π h of the posterior on the coefficients ξ given the observations is then where Φ h is defined in (4.5) and Z h is the normalization constant The same procedure can be applied seamlessly to the probabilistic case, thus obtaining a random posterior density π h over R N KL for the coefficient ξ.
The last detail missing is how to produce samples in order to obtain the approximations (4.6) and (4.7). Being the normalizations constant unknown, we employ Markov chain Monte Carlo techniques (MCMC) (see e.g. [24,Chapter 3] or [27,Chapter 6]), which proceed by generating an ergodic Markov chain whose invariant density is the desired posterior. Successive states of the aforementioned Markov chain then serve as samples from the posterior density. We choose to employ the Metropolis-Hastings (MH) algorithm, which we here briefly detail. The Markov chain is built employing a symmetric proposal q : for all x, y ∈ R N KL and such that for any fixed x ∈ R N KL the function q(x, ·) is a probability density, and with an acceptance-rejection strategy. In particular, given an initial guess ξ (1) , the algorithm proceeds for i = 2, . . . , N MC as (ii) Set ξ (i) =ξ (i) with probability α, and ξ (i) = ξ (i−1) with probability 1 − α, where Let us remark that the normalization constant Z h does not need to be known to run the algorithm, since we only compute ratios of posterior densities. Moreover, we note that the proposal distribution is the only tunable element of the MH algorithm. The easiest choice, at least for implementation, would be to fix q(x, ·) = N (x, σ 2 I) for some user-prescribed variance σ 2 . Unfortunately, the quality of the resulting Markov chain is not robust with respect to σ. In particular, if σ is too small, the probability to accept is too large and the Markov chain fails to effectively explore the posterior. At the other end of the spectrum, if σ is too large the probability of accepting a new sample reduces drastically, and the Markov chain presents a sticky behaviour. We therefore decide to employ the robust adaptive Metropolis algorithm (RAM) (see [48] for details), in which the proposal is q(x, ·) = N (x, Σ q ), where the covariance Σ q is adapted on the fly to obtain a user-specified final acceptance ratio, i.e. the ratio between the accepted and the total number samples, which should be roughly 25% (see e.g. [48]). Another viable option for the implementation of MCMC could have been the preconditioned Crank-Nicolson MCMC (pCN-MCMC) of [17,22], which is tailored for high-dimensional inverse problems.
For the probabilistic case, we perform a run of the MH, implemented with RAM proposal, with N in MC iterations for each one of the N out MC realizations of the random mesh, thus obtaining the approximation (4.7).
Remark 4.2. The total number of samples is given by in the probabilistic case by N in MC · N out MC . One could argue that, for a fair comparison between the probabilistic and the deterministic case in terms of computational cost, one should choose N MC ≈ N in MC · N out MC . In fact, since the "outer" Monte Carlo simulation can be performed in parallel, the correct scaling is N MC ≈ N in MC . Moreover, due to Remark 3.7, the number of random meshes does not need to be chosen excessively large.

Mean Eµ[κ]
Confidence Interval Figure 10: Numerical results for κ * 1 in Section 4.2.1. In all plots, the solid line represents the true conductivity, the dashed line is the posterior mean, and the shaded grey area is a confidence interval. In the first row, results are obtained by approximating the forward map with FEM, and in the second with the RM-FEM.

Numerical Experiments
In this section we present numerical experiments highlighting the beneficial effects of adopting the probabilistic framework of RM-FEM in the context of Bayesian inverse problems.

One-Dimensional Case
We first consider D = (0, 1) and solve the inverse problem presented above for two different true diffusion fields κ * . In both cases, we consider the prior on H to be given by N (0, Γ 0 ), with Γ −1 0 = −d 2 /dx 2 with homogeneous boundary conditions, so that the Bayesian inverse problem is well-posed. First, we consider κ * 1 = exp(ϑ * 1 ), where the log-conductivity ϑ * 1 ∈ H is given by are the first four ordered eigenpairs of Γ 0 . Second, we consider ϑ * 2 ∈ X ∩ H C , so that the true conductivity does not belong to the domain in which we solve the inverse problem, but it is still admissible for (4.1) to be well-posed. In particular, we consider the discontinuous conductivity We sample with the MH algorithm from the posterior distributions µ h and µ h , with N MC = 2 · 10 5 for µ h and with N out MC = 50 and N in MC = 2 · 10 5 for µ h . Knowing for the first conductivity κ * 1 that the true conductivity is fully determined by four coefficients, we fix the truncation index N KL = 4 in the Karhunen-Loève expansion (4.8). For the second conductivity κ * 2 , we fix N KL = 9. We then approximate the mean and pointwise standard deviation with (4.6) and (4.7) for the deterministic and probabilistic posteriors, respectively. Moreover, we arbitrarily fix a pointwise confidence interval at twice the standard deviation away from the mean. Numerical results are given in Fig. 10 and Fig. 11. Results highlight that for a coarse approximation, specifically for N = 10, the posterior distribution µ h is overly confident on the result. Indeed, the posterior mean fails to capture precisely the true conductivity in both the continuous and discontinuous case, and the confidence interval is extremely sharply concentrated around the mean. Conversely, the distribution µ h based on the probabilistic forward model accounts better for the uncertainty due to numerical discretization. Increasing the number of elements N , the mean computed under µ h and µ h tends to approximate better the true conductivity field. In particular, for N = 40 the posteriors µ h and µ h are already practically undistinguishable and are close to the true field. Moreover, let us remark that while the width of the confidence interval seems independent of N for µ h , it shrinks coherently to the discretization for µ h . Finally, we note that for κ * 2 even for larger values of N the posterior µ h seems to capture with its uncertainty local errors in the solution of the inverse problem. Indeed, the posterior mean is particularly off the true field on the left side of the domain, where the confidence interval is wider with respect to areas where the solution is more accurate.

Two-Dimensional Case
We consider now a two dimensional example on the domain D = (0, 1) 2 . We fix a Gaussian prior µ 0 on H for the log-conductivity ϑ chosen as µ 0 = N (0, Γ 0 ), where Γ 0 = −∆ −1.3 with homogeneous boundary conditions, so that the inverse problem is well-posed. We fix N KL = 6 and let the true conductivity κ * = exp(ϑ * ) in (4.1) be given by are the first six ordered eigenpairs of Γ 0 , and where ξ * i = (−1) i+1 · 10 for i = 1, 2, . . . , 6. Let us remark that ϑ * ∈ H. The right-hand side in (4.1) is chosen as f (x, y) = 8π 2 sin(2πx) sin(2πy). Synthetic observations are obtained by evaluating a reference solution on m = 50 random locations sampled from U(D) and then corrupted by an observational noise distributed as N (0, 10 −6 I). We then approximate the forward map in the inverse problem with the FEM and the RM-FEM. We choose a structured mesh T h as the one in Example 2.2 (or the second row of Fig. 1). In particular, in this case we let h denote the constant length of the short side of the triangular elements, i.e., the inverse of the number of subdivisions of each side of D.
In particular, we consider h = 0.1 · 2 −i , i = 0, 1, . . . , 3. The RM-FEM is implemented with p = 1 in (2.4) as per Theorem 2.9, and with an uniform choice for the random perturbations as the one described in Example 2.2.
Employing the notation introduced in Section 4.1, we then sample from the posterior distributions µ h and µ h employing the RAM method, considering only the first N KL = 6 coefficients in the KL expansion (4.8). In particular, we consider N MC = 10 5 samples for the deterministic case, and for the probabilistic case we generate N in MC = 10 5 samples for N out MC = 24 parallel chains, each corresponding to a realization of the random mesh in the RM-FEM. We then compute for each value of h the mean and standard deviation of the field computed under µ h (resp. µ h ) and denote their Monte Carlo approximations as κ µ h and σ κ µ h (resp. κ µ h , σ κ µ h ). Moreover, we consider the statistics of the 6-dimensional coefficient σ of the KL expansion, and denote by ξ µ h and σ ξ µ h the Monte Carlo approximation of mean and standard deviation computed under µ h (resp. ξ µ h , σ ξ µ h ). We show in Fig. 12 the posterior mean κ µ h for three values of h, compared to the truth κ * , and remark that the mean approximation is sensibly better for smaller values of h. The mean value under the probabilistic posterior µ h is not shown, as it is essentially equal to the deterministic case. The beneficial effect of employing the RM-FEM-based posterior distribution µ h , with respect to the FEM-based posterior µ h , consists of the approximate equalities which indicate that the error on the conductivity field, or on the coefficients of its KL expansion, are well represented by the uncertainty in the posterior distribution. This is shown in Fig. 12, where we notice that under µ h the standard deviation is practically independent of h and small with respect to the error on the solution of the inverse problem. Conversely, under µ h we have that the posterior standard deviation converges accordingly to the error, both for the L 2 -norm of the error on the mean and for the coefficients of the KL expansion.

Error Analysis for the RM-FEM
In this section, we present our a priori and our a posteriori error analysis for the RM-FEM. Let us remark that while the a priori error analysis is carried on for a general space dimension d and the adaptive algorithm has been shown to be efficient in higher dimensions (see Section 3.1), we present a rigorous a posteriori error analysis only in case d = 1. Conversely, in the a priori analysis we fix the coefficient p = 1 in (2.4), whereas in the a posteriori analysis we consider general perturbations, i.e., general coefficients p ≥ 1 in the same equality.

A Priori Error Estimates
We first prove the a priori error estimate given in Theorem 2.9. The convergence properties of the FEM for the elliptic problem (2.1) are well-established. In particular, without any additional assumptions on the exact solution, i.e., when u ∈ V , it holds u − u h V → 0 for h → 0. Under the more restrictive assumption u ∈ H 2 (D) ∩ V , we have a linear convergence rate, i.e.
for a constant C > 0, which is independent of h and u [7,12,35]. It is desirable that the RM-FEM is endowed with the same property. Moreover, we wish the error due to randomization to be balanced with the error due to the FEM discretization, which is shown in the proof of Theorem 2.9 below.
Proof of Theorem 2.9. Since (5.1) holds independently of the mesh, we have for a constant C independent of h and u and of the coefficient p in (2.4). Hence, by the triangle inequality we have for p = 1 , which shows the desired result.
Let us remark that we have shown above that the probabilistic solution converges with the same rate with respect to h in case p = 1, but we have not considered the case p > 1, for which the probabilistic term may be of higher order. Indeed, a preliminary theoretical and numerical investigation leads us to conjecture that so that, at least in the mean-square sense, the error due to randomization should converge faster than the error due to discretization if p > 1.

A Posteriori Error Analysis in the One-Dimensional Case
In this section we prove our main result for the a posteriori error estimator of the RM-FEM given in Definition 3.1, namely Theorem 3.6. Our goal is to prove in the one-dimensional case that the probabilistic a posteriori error estimators are reliable and efficient, i.e., that there exist positive constants C low and C up independent of h and u such that where κ ∈ L ∞ (D) satisfies κ(x) ≥ κ almost everywhere in D, and where we assume f ∈ L 1 (D). We recall that the notation for one-dimensional problems has been introduced and discussed in Example 2.2 and at the end of Section 3. Additionally, we introduce here for a function w which is piecewise constant on each K i ∈ T h the jump operator Our strategy for proving that the error estimator introduced in Definition 3.1 satisfies (5.3) relies on showing it is equivalent to known valid estimators. In particular, we consider the following estimator, defined in [5, Definition 6.3].
Definition 5.1. Let κ be the diffusion coefficient of (5.4) satisfy κ ∈ C 0 (D) and κ ≥ κ > 0. We define the error estimator Clearly, the quantity E h is computable up to quadrature error due to the approximation of the local estimators η j . Let us finally introduce more precisely the higher-order quantity Λ appearing in (3.2), i.e., for a constant C independent of h and of the solution u. If moreover the family of meshes T h is λ-quasi-uniform and if κ ∈ C 2 (D) and f ∈ C 1 (D) then, up to higher order terms, it holds for constants C low , C up independent of h and u.
We recall that in the one dimensional case the probabliistic error estimators for the RM-FEM are given by .
Our strategy to prove Theorem 3.6 relies on showing that the deterministic estimator E h of Definition 5.1, as well as its probabilistic counterparts E h,1 and E h,2 above are all equivalent to the quantity i.e., the sum of all squared jumps of the derivatives on the internal nodes. We first prove the equivalence for E h,1 .

Lemma 5.3. Let Assumption 2.1 hold. Then, if the mesh is λ-quasi-uniform it holds
where E h,1 is given in Definition 3.1.
Proof. Let K i , i = 1, . . . , N , be a generic element of the perturbed mesh and let us compute the derivative of the interpolant on K i , which is given by where an exact Taylor expansion allows to compute Hence, it holds which we can rewrite rearranging terms as It is clear then that the expression above depends on the signs of the variables α i−1 and α i . For simplicity of notation, we therefore introduce the events A In the trivial case α i−1 > 0 and α i < 0, i.e., if A i−1,i occurs, the equality (5.6) simplifies to

Since in this case
where we denote We thus have an expression for E I i due to (5.7). We now turn to we have by the law of total expectation Let us remark that adding and subtracting u h Ki yields The same computations employed for I i allow to conclude that which, replaced into (5.12) gives the final expression for E I i−1,i . Similarly, for I i+1,i we have We now reassemble the quantity I i + I i−1,i + I i+1,i by grouping terms with regards to their conditioning on the sign of (α i−1 , α i ). In particular, some algebraic simplifications yield We now can compute the estimator E h,1 by summing its local contributions, as in where J 1 and J 2 are given by Let us consider J 1 and J 2 separately. Rearranging the sum, noticing that under Assumption 2.1(i) it holds E [α i | α i > 0] = − E [α i | α i < 0] = E |α i | and recalling that α i = (h i h −1 ) pᾱ i , we obtain Now, let us remark that if the mesh is λ-quasi-uniform, it holds We now turn to J 2 . Clearly, we have J 2 ≤ 0, which implies the desired upper bound together with (5.13). For the lower bound, we remark that in both cases A i−1,i occurs, we have x i − x i−1 ≥ h i . Hence, simplifying the conditioning in the first and second terms, we obtain We now consider ξ i given in (5.11) and use (a + b) 2 ≤ 2(a 2 + b 2 ) for a = α i−1 u h xi−1 and b = α i u h xi to obtain Therefore Rewriting the sum and replacing the definition of α i yields where we applied the law of total expectation on the second line. Finally, we haveh i ≤ h and p ≥ 1, which yield Combining this with (5.13) then yields the desired lower bound and thus concludes the proof.
Let us remark that the coefficient appearing in the lower bound of Lemma 5.3 is positive if Assumption 3.4 holds. We now prove the equivalence of the estimator E h,2 given in Definition 3.1 with J (u h ).

Lemma 5.4. Let Assumption 2.1 hold and let the mesh T h be λ-quasi uniform. Then, it holds
where E h,2 is given in Definition 3.1.
We now remark that if the mesh T h is λ-quasi-uniform and under the assumptions on κ it holds for k ∈ {0, 1} m 2 (1 + λ) 2 which, in turn, implies We now focus on the upper bound. Reassembling the global error estimator E h , we have where we recallh j = min{h j , h j+1 }, so that h j + h j+1 ≤ (1 + λ)h j . We conclude the proof proceeding similarly for the lower bound as in Lemma 5.4.
We can finally prove Theorem 3.6 and conclude the error analysis.
Proof of Theorem 3.6. Let us first consider E h,1 . Under Assumption 3.4, we have for the lower bound of Lemma 5.3 for a constant C > 0. Hence, due to Lemma 5.5 we have that there exists a constant C such that

Conclusion
We have introduced a novel probabilistic methods for PDEs based on the FEM and random meshes, the RM-FEM. We demonstrated how our methodology can be successful when employed in pipelines of computations, such as Bayesian inverse problems. We also show a rigorous use of probabilistic methods for a posteriori error estimators, often speculated in the field. Extending such analysis to PN methods for ODEs would be of interest, thus creating a link between the guiding principles of PN and more classical theories. Generalizing the analysis of the RM-FEM to higher-dimensional PDEs as well as for parabolic or hyperbolic problems represent also interesting future work.