{'fraction_non_alphanumeric': 0.04559594251901944, 'fraction_numerical': 0.013119188503803888, 'mean_word_length': 4.324122412241224, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Model reduction for fractional elliptic problems using Kato's formula

We propose a novel numerical algorithm utilizing model reduction for computing solutions to stationary partial differential equations involving the spectral fractional Laplacian. Our approach utilizes a known characterization of the solution in terms of an integral of solutions to classical elliptic problems. We reformulate this integral into an expression whose continuous and discrete formulations are stable; the discrete formulations are stable independent of all discretization parameters. We subsequently apply the reduced basis method to accomplish model order reduction for the integrand. Our choice of quadrature in discretization of the integral is a global Gaussian quadrature rule that we observe is more efficient than previously proposed quadrature rules. Finally, the model reduction approach enables one to compute solutions to multi-query fractional Laplace problems with order of magnitude less cost than a traditional solver.


1.
Introduction. Differential equations involving fractional derivative powers have gained in popularity in recent years. These non-classical differential equations have shown potential to model nonlocal and time-delay effects, making them good candidates for modeling hysteretic and globally-coupled phenomena. For example, fractional differential equations have recently been used to model fluid mechanics, arterial blood flow, cardiac ischemia, and geophysics; as components of optimal control problems; and as ingredients in image denoising and image segmentation [38,30,22,2,5]. In particular, optimal control problems frequently require optimization of an objective function that involves the solution to a partial differential equation (PDE). Algorithms for solving the optimization problem would require many evaluations of the objective, hence of the PDE. Thus, for control problems it is necessary to ensure that the computational solution of the PDE can be obtained as quickly as possible. In this paper we focus on the acceleration of computational algorithms for such PDEs, in particular those involving fractional elliptic operators, which are defined via spectral expansions; a prototypical example of an operator that we use throughout this paper is the fractional Laplacian.
With ∆ the "standard" (classical) Laplacian operator, we are interested in computing the solution u to the partial differential equation valid on some physical domain Ω (which is open, bounded, and Lipschitz) with appropriate boundary conditions on ∂Ω, where s ∈ (0, 1) is the fractional order. The precise definition of the fractional operator (−∆) s involves the spectral expansion of the local operator −∆ (with appropriate associated boundary conditions), which we more precisely describe in Section 2. There are already several numerical algorithms for computing the solution to such an equation: • Perhaps the most conceptually straightforward idea is to use the spectral expansion definition of (−∆) s to devise a scheme that computes solutions using the spectral expansion of the associated discretized operators [27,28,45,40]. The disadvantage of this approach is that the procedure is expensive, requiring a full eigendecomposition of a potentially very large matrix. • A second approach uses an extension procedure to write the nonlocal ddimensional PDE as a (d + 1)-dimensional local PDE [34,13,41]. This latter PDE can be solved with existing methods, although some nontrivial tailoring of existing numerical methods is needed [35]. The challenge with this approach is that the spatial dimension is increased, and the extended PDE is degenerate, requiring specialized numerical methods. More recent efforts have significantly reduced the computational cost of this method. These methods exploit the tensor product structure that occurs as a result of the additional dimension [33,1] but rely on discretizations that depend on the exponent s, i.e., essentially assume a non-zero lower bound for s. • A final approach that we use as the starting point for the method proposed in this paper is an integral operator approach, which writes the solution (i.e., the inverse of the operator) as a type of Dunford-Taylor integral involving the resolvent of the local operator: See [29, Theorem 2 with λ = 0]. Numerical approaches for solving nonlocal PDEs using this formula discretize the integral in (1) with quadrature and require many local ("standard") PDE solves (equal to the number of quadrature points) in order to compute the solution to the fractional problem [11,6,9]. However, this results in an algorithm that can require, depending on the problem, many or hundreds of queries of an existing PDE solver to compute the solution to the fractional problem. The challenge with this approach is that often many local PDE solves may be necessary to ensure accuracy for a single solution of the fractional PDE, making this approach quite expensive compared to traditional solvers. For certain operators these local PDE solves may be accomplished in parallel, but this does not diminish the overall cost. Other operators require coupling of these solves, e.g., in cases when a lifting strategy is used to solve non-homogenous nonlocal Poisson equations [43], making parallel approaches more difficult. In this paper, we develop a novel model reduction algorithm for the third approach listed above to substantially alleviate the cost of traditional PDE solves. We concentrate on this approach for the spectral definition of the fractional Laplacian in this paper, but due to a similar integral formulation for the integral fractional Laplacian [10], our approach would extend to more general cases as well.
This paper is not the first strategy for model reduction for fractional elliptic problems. The authors in [4] provide a model reduction strategy, applied to the second (extension) approach listed above. More recently, the work in [18,19] employs a reduced basis approach by interpolating operator norms. Reduced order methods for this problem using the integral operator approach are also investigated in [9]. However, low-rank structure in solution sets to fractional problems has been empirically noted even earlier [44]. For problems involving nonlocal integral kernels, the authors in [25] also proposed a reduced basis approach, but use a different strategy to perform model reduction. Our contributions in this article are as follows: • We provide a rearrangement of the Dunford-Taylor integral considered in [11] that improves numerical stability. We first show that the analytical solution has s-independent L 2 stability bounds. This stability extends to the numerical discretization, independent of all discretization parameters. See Lemma 3.1 and Proposition 3. • Our approach to discretize the Dunford-Taylor integral is a novel application of a global Gaussian quadrature rule. Our numerical results suggest that our quadrature choice is more efficient than previously proposed choices, cf. Figure  4. We cannot provide an analytical error bound in terms of the number of quadrature points, but we do provide a rigorous, computable error certificate; see Proposition 4. Previous work has required a number of quadrature points proportional to max(1/s, 1/(1 − s)) in order to obtain a specified level of accuracy. Our empirical results suggest that our approach also suffers from this limitation, see Figure 2. • We employ the reduced basis method (RBM) to effect model reduction which, after a single offline computational investment, can accelerate subsequent computations of (s, f ) → u by at least two orders of magnitude. The offline portion of this algorithm requires approximately as much time as a single (s, f ) → u solve using the traditional Dunford-Taylor approach; see Algorithm 3.4. • We provide a rigorous a posteriori error estimate for our solution computed via model reduction. This error estimate is computed as a by-product of the offline investment, and is therefore directly available; see Theorem 5.2. Regarding alternative model reduction approaches to solve this problem: After initial dissemination of the first draft of this manuscript, the authors in [9] provided a complementary approach. For s 0 some fixed constant, [9] provides rigorous convergence analysis for energy norms for s ≥ s 0 > 0, with estimates depending on 1/s 0 . The space interpolation-based approach in [18,19] likewise provides convergence the L 2 norm for the discretized setting. Although comparably strong convergence analysis is missing in our work (i.e., we do not yet have rigorous convergence analysis available), we do accomplish different goals: We establish L 2 -stability of the fractional problem uniformly as s ↓ 0 both at the discrete and continuous levels. In addition, we provide rigorous, computable error certification for the fractional problem (and not just for the auxilliary problem defined by the Dunford-Taylor integral); these results distinguish our study from the similar approach in [9]. In comparison to extension approaches, the additional reduced order modeling based approach considered in [4] requires solving a PDE in an additional (extended) dimension, and due to this requirement is less efficient, even if combined with efficiency gains in full-order models achieved in [33] and related work.
We remark that while we study PDEs with an operator of the form (−∆) s , all our results extend to more general fractional elliptic operators. See Remarks 1 and 4.
This paper is structured as follows. Section 2 lays out our notation and describes the problem. Section 3 describes a new algorithm for expressing and computing the Dunford-Taylor solution that was first proposed in [11]; this new algorithm is a novel utilization of a Gaussian quadrature rule. The accuracy of this quadrature rule is empirically investigated in Section 4. Section 5 is our main algorithmic section that introduces an RBM algorithm for model reduction that computationally accelerates the fractional PDE algorithm from Section 3, and provides a computable error certificate for the model reduction. Finally, section 6 includes numerical results that demonstrate our new algorithms on a two-dimensional fractional Laplace problem and compares our algorithm against the predecessor in [11].
2. Notation and setup. Vectors will be denoted in lowercase bold, and matrices in uppercase bold, e.g., x and A, respectively. If M is a symmetric positive definite matrix, we define , and x is the standard Euclidean norm. The matrix norm A is the standard induced 2 norm on matrices. If A is symmetric, then λ min (A) denotes the smallest (real) eigenvalue of A. If both A and B are symmetric positive definite matrices in R N ×N , we define the smallest generalized eigenvalue of (A, B) as Note that under these assumptions on A and B, the above expression is equal to the smallest λ such that Ax = λBx has a nontrivial solution x, and also we have that where B 1/2 is the symmetric positive definite matrix square root of B. Similarly, we use the notation λ max (·) and λ max (·, ·) to denote maximum eigenvalues.
Consider a bounded domain Ω ⊂ R d with Lipschitz boundary ∂Ω; for spatial domains we are mainly concerned with d ≤ 3 but this restriction is not necessary. We have We will often write L 2 = L 2 (Ω), and we define ∇w, ∇v = d j=1 ∂ ∂xj w, ∂ ∂xj v , which induces a definition for the L 2 norm ∇v L 2 of vector-valued functions. For brevity will also write v = v L 2 and ∇v = ∇v L 2 . The standard Laplace eigenvalue problem on Ω with Dirichlet boundary conditions, −∆u = λu, x ∈ Ω, yields an infinite sequence of eigenvalues 0 < λ 1 ≤ λ 2 · · · with associated eigenfunctions {φ n } ∞ n=1 . Here, and in all the following, the differential operator ∆ operates on the x variable. The spectral theorem ensures that the eigenfunctions enjoy L 2 (Ω)-orthogonality and completeness, so that where we have further assumed that each φ n has unit L 2 (Ω) norm.
2.1. The fractional Laplace problem. In this section we provide the definition of the fractional operator (−∆) s , and summarize some regularity properties of solutions. Let s ∈ (0, 1). In this section we describe the spectral definition of the fractional operator (−∆) s , on bounded domains, supplemented with homogeneous Dirichlet boundary conditions, see [7] for the inhomogeneous case. Related definitions of similar nonlocal or fractional operators can be found in the literature [31]. If u has an expansion in eigenfunctions, then a formal definition for application of the fractional operator is, We likewise use λ n to define Sobolev spaces of fractional order. With s ∈ (0, 1): where H s (Ω) * denotes the dual space of H s (Ω). With these definitions, (−∆) s : H s (Ω) → H −s (Ω). For the relation of H s (Ω) to the standard fractional order Sobolev space H s (Ω), see [32, Chapter 1, Section 9] or [21] for a more modern survey. Given data f ∈ L 2 (Ω), our main goal is to compute the solution u to for arbitrary s ∈ (0, 1). We emphasize that only for s > 1/2 can the boundary conditions in (5) be understood in the classical trace sense. This is due to the fact that, for smooth enough boundaries, H s (Ω) = H s 0 (Ω) for s > 1/2. For s < 1/2, H s (Ω) = H s 0 (Ω) = H s (Ω) and in this case the boundary conditions are not meaningful in the classical trace sense. One possible interpretation of the boundary conditions in this case is using the integration-by-parts formula, as described in [7, (Ω) i.e., the Lions-Magenes space, in this case, it follows from the definition of the space that the functions need to have certain decay at the boundary.
One can assume f ∈ H −s , yielding a solution u in H s . Note therefore that f ∈ L 2 is a stronger requirement than necessary for this fractional problem, and is even stronger than usual for local elliptic problems. We make this stronger assumption to be assured, for fixed f , that u ∈ L 2 uniformly for s ∈ (0, 1). Note that this is a reasonable assumption in the optimal control context where controls are frequently assumed to be smoother than the explicit image of a differential operator. Notationally, we omit showing explicit dependence of u on the spatial variable x, and only show dependence on the fractional order s, which is a parameter. This convention will be used in the remainder of this paper when considering solutions to parameterized PDE's: notational dependence on parameters will be explicit, but that on the spatial variable x will be implicit. Therefore, we let u(s) ∈ L 2 denote the solution u to (5) for a fixed value of s. We will be interested in developing a computational algorithm for computing the family or manifold of solutions, If f ∈ L 2 , then the solution u to (5) has H 2s membership. In order to ensure that u is also in the classical Sobolev space H 2s , we need to assume the domain to be more regular, for example, quasi-convex. Recall, that quasi-convexity ensures that [23,Theorem 10.4]. See also [24] for a related discussion on solution regularity. Our goal is to study the manifold U so that the natural function space is ∩ s∈(0,1) H 2s . Therefore, all of our investigations will assume f ∈ L 2 and study solutions u(s) as elements of L 2 .
In the remainder of this document we will consider the problem (5) with homogeneous Dirichlet boundary conditions. The inhomogeneous boundary case may be handled using the lifting technique in [7], which requires only a small modification of the algorithm proposed here. While the previous discussion has centered on the definition of the nonlocal operator and the associated PDE, our main analytical and computational tool will be an equivalent formulation via Kato's formula.

2.2.
Kato's integral solution of (5). The following remarkable result provides an appealing formula for the solution u to (5): which is a reformulation of Kato's formula (1) from [29, Theorem 2 with λ = 0]. This representation was first exploited in [11] for designing numerical algorithms, and is derived via a special kind of Dunford-Taylor integral. To write the above more explicitly, define q(y), for fixed y ∈ R, as the solution to the local y-parameterized Then u is given by This representation reveals that u is actually just an integral of solutions q to local Laplace-type problems. A solution method employing a discretization of the above formula then only requires solves of (classical) local PDE's in order to solve the nonlocal problem (5). The straightforward way to compute the solution via (6b) is to approximate the integral with a quadrature rule. This would require computing solutions q(·; y) to the PDE (6a) for many values of the parameter y.
More precisely, let {y m , τ m } M m=1 be a quadrature rule for approximating the integral in (6b). We will give precise choices for this quadrature rule soon. Then we can approximate the solution u(s) as One then needs only to compute the ensemble of functions {q(y m )} M m=1 , which are solutions to local PDEs, in order to approximate the solution to the fractional problem. This is the particular approach adopted in [11], wherein a sinc quadrature rule is adopted, and associated error bounds are derived.
The observation we make in this paper is that the approach above requires approximately M times the work of solving a local problem; when M is large (which can be required when s is small), this may become computationally prohibitive, see for instance [42]. However, there are by-now standard model reduction approaches that allow one to efficiently compute solutions to parameterized PDE's when the number of queries M is very large; one such method that is directly applicable (to some extent) here is the reduced basis method, which we exploit in Section 5. The next section expresses (6) in a more computationally robust formulation, and proposes a new kind of quadrature for y-discretization. We observe in Section 6 that our new quadrature approach is much more efficient than the most efficient strategy considered in [11].
3. Fractional Laplace solutions via integral formulation. Recall that given s ∈ (0, 1), we seek to evaluate (6b), which defines the solution u(s) to the fractional Laplace problem (5). In this section we describe our algorithm for doing so, which discretizes the y variable using quadrature. The main components of this algorithm come in two stages: first we describe a partitioning formulation for the y variable, followed by a quadrature discretization of the y integral.
The approach described in this section augments the approach presented in [11]; our improvements include s-independent stability in both the continuous and discrete case. For small values of s and large y-quadrature rule size, the algorithm in [11] results in discrete operators whose norm becomes very large, which can be problematic for numerical implementation. In our reformulation, the discrete operators are uniformly bounded in s (for any quadrature rule). Second, we replace the sinc quadrature in [11] with a Gaussian quadrature rule. (The authors in [11] also propose using Gaussian quadrature rules, but theirs is a composite rule, whereas ours is a global rule and is designed differently.) Our results in section 6 indicate that this quadrature rule is substantially more efficient than sinc quadrature.
The sections below perform the following tasks: • Section 3.1 rewrites the integral formulation (6b) in a form that we actually discretize. Lemma 3.1 is this rewritten expression. This rewritten expression allows one to easily observe L 2 stability results for solutions. This is Proposition 1. These stability results can naturally also be observed from (6b), to analyze the algorithm it is more convenient to use the rewritten form. • Section 3.2 imposes a Galerkin-based spatial discretization on certain local auxilliary PDEs, which translates into a spatial discretization for the fractional solution u. Proposition 2 establishes that the discrete solution inherits the stability of the continuous problem.
• Section 3.3 introduces the particular global quadrature rule that we employ to discretize the y variable, and establishes a stability result for this fully discrete scheme in Proposition 3.
To accomplish this, we introduce a new parameterized PDE for a solution w that is closely related (but not identical) to q from (6a) where (α, β) ∈ (0, ∞) × [0, ∞) is the parameter. In our computational setting, we will only require (α, β) ∈ (0, 1] 2 . Comparing (7) with (6a), we see that q(y) = w(1, e y ). We now define w ± (y) as two specializations of w that will be used in the following: for x ∈ Ω, with boundary conditions w ± (y) = 0 for x ∈ ∂Ω. We note that for y ∈ [0, ∞), we have q(y) = e −y w + (y). We can formulate the solution to (5) in terms of these new quantities.
where β 0 , s ± , and W are defined as Proof. Beginning with (6b), we have completing the proof.
Given the domain Ω, we define C Ω as the domain's Poincaré constant, i.e., the smallest constant such that for every Note that C Ω = 1/λ 2 1 , where λ 1 is the minimal eigenvalue of the Laplacian from (2). The space H 1 0 (Ω) is the standard Sobolev space of zero-trace L 2 (Ω) functions whose gradients are also in L 2 . In what follows, we will also need the following quantity, which also depends only on Ω. Fixing (α, β), the weak formulation of (7) seeks a solution w(α, β) = w ∈ H 1 0 (Ω) as the unique function satisfying the Galerkin formulation, with the bilinear form a(·, ·; α, β) defined as With α > 0 and β ≥ 0, the Poincarè inequality (11) ensures that the coercivity holds for some c > 0 uniformly in β, so that standard Lax-Milgram theory then yields a unique H 1 0 (Ω) solution. With this setup, we can demonstrate the utility of a formula like (9) by deriving an s-independent L 2 stability estimate for solutions to (5).
Taking v = w + (y) and using the Cauchy-Schwarz and Poincaré inequalities results in f w + (y) ≥ f, w + (y) = e −y ∇w + (y), ∇w + (y) + w + (y), w + (y) We thus obtain Therefore, taking the L 2 norm in (9) and using the triangle inequality yields , where the third inequality uses the fact that W is a probability density on [0, ∞). From the above, we immediately obtain the desired result by noting that The proof is complete.
Note that formulas like (9) are not the only way to show such stability results, but Proposition 1 is meant to illustrate one possible use of these relations. Our main utility of such results will be to design a numerical scheme.
3.2. Spatial discretization. In this section we apply a spatial discretization to the result of Lemma 3.1. We proceed to discretize (14) using a finite element method. Let T Ω be a conforming triangulation of Ω with K elements. We assume each element in the triangulation is isoparametrically equivalent to a standard canonical triangle/tetrahedron. For a fixed polynomial degree k ≥ 1, we define the finite element space where P k (e) is the space of polynomials of degree k or less over the element e ∈ T Ω . Let N = dim V . The finite element-discretized version of (13) is the Galerkin Let w N (α, β) be expressed as a linear expansion, where {ψ n } N n=1 is a basis for V , e.g., a basis comprised of compactly supported piecewise polynomials. Collecting the linear degrees of freedom of w N (y) ∈ V in the N -dimensional vector w N , then this vector satisfies the linear system, where the N × N matrix A(α, β) has entries, for j, k = 1, . . . , N . Above we have defined the N × N y-independent stiffness and mass matrices S and M , respectively. Both S and M are symmetric and positive-definite. The matrices S and M can be used to define a "discretized" Poincaré constant C N by using a standard Rayleigh quotient argument: hence leading to the inequalities, Our estimates for x-discrete quantities will involve C N , but we will sometimes use the above inequality to bound quantities in terms of C Ω . Bounds involving C Ω emphasize independence of the x-discretization. Bounds involving C N emphasize the explicit computability of the bounds, since C N is equal to an extremal eigenvalue of finite element matrices, which is computable with iterative eigenvalue solvers.
The maximum generalized eigenvalue of (S, M ) will also play a small role in our estimates. In analogy with (19) we define Note that K N in general tends to 0 as N ↑ ∞.
From the discretization of w we derive discretizations of w ± (y) defined in (8). Thus, finite element discretizations for w ± are specializations of those for w. In particular, we define so that, akin to (17), we have We can now codify the fact that the solutions w N ± (y) are L 2 -stable uniformly in y.
Proof. The result can be obtained by considering the discrete form (22). To begin we relate the L 2 norm of f to the Euclidean 2 norm of f . Let P V denote the L 2 -orthogonal projector onto V . Then: Thus we have Now since M is symmetric and positive-definite, it has a unique symmetric positivedefinite square root M 1/2 . Thus: This in turn implies: Since M −1/2 SM −1/2 is symmetric and positive-definite, we have λ min e −y M −1/2 SM −1/2 + I ≥ λ min (I) = 1. Therefore, which, when combined with (24) yields (23a). A similar computation for w − (y) yields (23b) by using the definition of C N in (19).
The result above gives the stability of an algorithm that uses w N ± (y) as a spatial discretization. In particular, consider the following semi-discrete approximation of u(s), A fully discrete scheme, introduced in the next section, would discretize the y variable. The following result mirrors the stability estimate of Proposition 1, showing s-uniform L 2 stability of the semi-discrete solution.
The proof, which we omit, is almost exactly the same as for Proposition 1 using the discrete stability estimates (23); the only novelty is the need to exercise the inequality (19).
Having described the spatial discretization, we now proceed to describe a discretization for the y-integrals in (9).

3.3.
Quadrature for the y-integral. We will use an M -point W -Gaussian quadrature rule to discretize the integrals in (9). The weight function W (y) is a weight function associated with a classical family of orthogonal polynomials: Laguerre polynomials. Let {p n } n≥0 denote the family of Laguerre polynomials, orthonormal under the weight W , i.e., where δ m,n is the Kronecker delta function. Like all orthogonal polynomials, the family {p n } n≥0 satisfies a three-term recurrence formula, yp n (y) = b n+1 p n+1 (y) + a n+1 p n (y) + b n p n−1 (y), n ≥ 0, with p 0 ≡ 1 and p −1 ≡ 0. The recurrence coefficients (a n , b n ) are explicitly known: b 0 = 1, b n = n a n = 2n − 1, n ≥ 1.
Among the properties of orthogonal polynomial families is the existence of a unique M -point quadrature rule with optimal polynomial exactness, the Gaussian quadrature rule. This rule has abscissae and weights, (y m , τ m ) M m=1 , respectively, and integrates polynomials up to degree 2M − 1 exactly: p ∈ span{1, y, . . . , y 2M −1 }.
Although y j and τ j depend on the value of M , we omit this explicit notational dependence. This rule can also be easily computed with knowledge of the recurrence coefficients. Defining the M × M symmetric tridiagonal Jacobi matrix, consider its associated eigenvalue decomposition, where V is unitary since J M is symmetric. The Gaussian quadrature rule can be computed from these quantities. In particular, where v j,1 is the first component of the vector v j . To compute an M -point quadrature rule for W thus requires a size-M eigenvalue computation. Since W is a probability density function, then likewise M j=1 τ j = 1, and furthermore each τ j is non-negative by (26).
The cost of computing the above M -point quadrature rule is dominated by the cost of computing the spectrum of symmetric tridiagonal matrix. For computing only the eigenvalues, many algorithms can compute this with O(M 2 ) cost [36]. . We apply these rules to the integrals in (25), resulting in the fully discrete approximation, We emphasize that the y-discretization does not suffer from any numerical instabilities as s ↑ 1 or s ↓ 0: the weights τ j,± are positive, no larger than 1, and independent of s, β 0 (1 − s) and β 0 (s) are just sinc functions, and w ± (y) has bounded L 2 norm for all y ≥ 0, i.e., for all inputs. The following codifies this stability.
Proof. The proof is very similar to the proof for proposition (1). Take the · M norm on both sides of (27), use the triangle inequality and (23), and note that the quadrature weights τ j,± all satisfy 0 ≤ τ j,± ≤ 1 since they are all positive and W is a probability density. Note that (28) holds for any quadrature rule for the y variable if augmented by a multiplicative constant equal to the quadrature condition number (sum of absolute value of weights).
Just as in [11], assuming we have the ability to compute w N ± , then the formulation above immediately yields an algorithm to compute u N (s) via (27).
Remark 1. All of our results extend to the case when we solve (5) but replacing −∆ with a general elliptic operator E that (i) satisfies the coercivity condition Ev, v ≥ α v H 1 0 for α > 0 and (ii) can be associated with a symmetric variational bilinear form. In this more general case, we need only replace all the instances of C 2 Ω with C 2 Ω /α. The matrix S should likewise be replaced with the associated matrix defined from the bilinear form of E.
3.4. Algorithm summary. The sections above identify an algorithm for computing u N (s) in (27). We summarize the procedure in Algorithm 1. The discrete solution adheres to the stability bound in Proposition 3. Note that we have not yet described how one should decide on values for M ± . We provide a concrete computational strategy for accomplishing this in the next section. However, one of our goals in our numerical results section is to compare our algorithm to existing ones, which determine M ± using (48).

4:
for j ← 1 to M − do 5: (22)  Compute u N (s) from (27). 9: return u N (s) 4. Error due to quadrature discretization. The formulation (27) is our numerical approximation to compute solutions to (5). This formulation is a discretization over both the x and y variables (via a finite element formulation and a quadrature rule, respectively). To understand the error that the y discretization contributes, we analyze the discrepancy between u N (s) and u N (s). To proceed, we need an auxiliary function that measures the absolute error between a size-M quadrature rule (y j , τ j ) We also need to define intervals on the real line enclosing the spectrum of some discretized operators. Recalling the definitions of C N and K N in (19) and (21), respectively, define intervals I and I ± as The error committed by the quadrature rule can be understood by studying the quantity, The precise statement is as follows.
Proposition 4. Assume f ∈ L 2 . Then for each s ∈ (0, 1), and therefore, Proof. The same argument that produces the relations (24) implies that so we proceed to study the quantity on the right-hand side. The difference between the "−" integral contributions in u N (s) − u N (s) is proportional to We can express the solution w N − (y) as Then further manipulation of the w N − (y) expression yields where H(y) is a diagonal matrix having entries (H(y)) j,j = 1 t j + e −y = The first inequality is sub-multiplicativity of the · matrix norm, and the first equality uses the invariance of the same norm under unitary transformations. A similar computation for the "+" quantities yields The combination of these results implies Using the inequality (20) yields the result.
The summation on the right-hand side of (31) can be computed independent of the data f , and requires only knowledge of the extremal eigenvalues of the discrete operator, cf. (19) and (21). While we cannot at present provide a theoretical estimate of this error, we numerically investigate the behavior of this error on M ± in our numerical results section.
Remark 2. Comparing (31) with the stability bound (28) suggests that many of the factors in (31) appear due to bounding the error relative to u N L 2 . Thus the supremum over g is the factor that arises due to the quadrature error.
This again suggests that, independent of all discretization parameters, our numerical algorithm is stable. However, a rigorous convergence analysis for our quadrature rule is not yet available.

4.1.
Empirical behavior of quadrature error. The main result from Proposition 4 is that the error in the fully discrete approximation (27) that is due to the y-quadrature discretization is computable without solving any PDE's, assuming that the extremal generalized eigenvalues of (S, M ), coded in the quantities C N and K N , are known. In particular, this implies that most details of the spatial discretization need not be utilized to understand the quadrature error; we only require extremal eigenvalues of discretized operators. We empirically investigate the accuracy of the quadrature rule in this section. Throughout our tests, we will use the following values: Since C N is bounded above by the analytical Poincaré constant of the domain C Ω , then choosing this O(1) quantity for C N is reasonable. The choices above make the intervals I ± defined in Proposition 4 explicit. The finite element discretization from our numerical experiments in Section 6 results in values C 2 N = 0.0506 and K 2 N = 2.36 × 10 −6 . We note that G ± in (30) can be numerically approximated for each (M, s) by replacing the supremum over a with the maximum over a discrete mesh. We compute the supremums in G ± by discretizing the intervals I ± with 200 logarithmically spaced points. (I.e., log I ± is replaced with 200 equispaced points.) This discretization then allows us to compute G ± , and hence allows us to compute approximations to the bound in (31). In Figure 1, we show the behavior of G ± as a function of (M, s). We note that ensuring small values of G + requires more quadrature points when s is close to 0. In contrast, controlling G − requires more quadrature points when s is close to 1. However, the behavior of G + for small s + = s is more restrictive than the behavior of G − for small s − = 1 − s. Thus, we expect that G + is the term that requires more computational investment to guarantee a certain error level. Figure 1. Values of log 10 G ± defined in (30) as a function of (M, s). We show s-dependence as log 10 (1/s) since the error behavior for small s is the most restrictive. We observe that, for fixed M , G ± has large values when s ± is small.
To explore this further, we define the smallest value 2 of M needed to assure a given error level δ: For δ = 10 −2 , 10 −4 , and 10 −6 , we display the values of these M quantities in Figure  2. We see that for small values of s, the requisite number of points M scales like 1/s. In particular, for small s, more effort (quadrature points) is allocated to M + , but for small 1 − s, comparatively more effort is allocated to M − . Therefore, the   (33) for various values of the tolerance δ. We show s-dependence as log 10 (1/s) since the error behavior for small s is most restrictive. For visual reference, a 1/s curve is also plotted. We see that for small values of s ± , the corresponding value of M ± is large.
number of quadrature queries M − + M + in the fully discrete scheme (27) can be quite large when the fractional order s is very small. This general observation, including the 1/s-type behavior shown in Figure 2, is consistent with earlier work [6]. Thus, the number of local PDE solutions M = M − + M + needed to compute an accurate solution is large. This motivates a need to make these solves more efficient; we achieve this in the next section via model reduction.

5.
Model reduction for the integral formulation. This section proposes an augmentation of the algorithm in the previous section. The cost of computing the fully discrete solution (27) is essentially M − + M + queries of finite element solvers for w N ± . In practice one can require M − + M + ∼ 100, cf. Figure 2 and earlier work [11,42], resulting in a substantial computational cost if the cost of computing w N ± is high.
We observe that the formulations (8) for w ± (and also (22) for the discrete counterparts w N ± ) are quintessential examples of parameterized PDE's where RBM algorithms are used to accelerate solution queries. Thus, RBM can be used to ameliorate the cost of performing M − + M + queries of these PDE's. In RBM terminology, an available expensive discrete solution is called a truth solution. Thus, our truth solutions for the auxiliary PDE problem for w ± are w N ± defined in (22). The associated truth solution for u(s) is (27). The purpose of RBM procedures is to diminish the cost of evaluating the truth solution.
5.1. Reduced basis methods. Let L be a (local) differential operator, and consider the following PDE parameterized by a Euclidean parameter y ∈ D ⊂ R p : where x is the spatial variable and y is a parameter. The operator L is differential in the x variable. For example, the PDE defining w(y) from (7) can be written as (35) with the operator, We assume that (35) is well-posed for each y ∈ D. For a fixed y, one usually develops an x-discretization with N 1 degrees of freedom yielding a solution w N with membership in an N -dimensional subspace. For us, this is the discretization defined in section 3.2. We assume that N is large enough so that where is a user-prescribed tolerance. Thus, the map y → w N (y) ≈ w(y) requires algorithms whose complexity is dependent on N . 3 Such an algorithm that performs the operation y → w N (y) is called a truth approximation or solver.
The reduced basis method (RBM) is a thematic collection of model reduction strategies for parameterized PDEs that compute an emulator y → w N (y) ≈ w N (y), whose complexity behaves like O(N 3 ) or O(N 2 ), where N N . For N/N sufficiently small, this can result in an emulator w N whose evaluation is substantially cheaper than the truth approximation w N . The RBM emulator takes the form, w N (y) := N k=1 c N,k (y)φ k , φ k ∈ V N := span w N (y 1 ), . . . , w N (y N ) ⊂ V, (37) where {y k } N k=1 are particular parameter values that are chosen during the RBM construction procedure. The success of RBM algorithms rely on three main components: • The condition that the manifold of solutions, is "low rank". The mathematically precise statement of this is that the Kolmogorov N -width of the manifold, decays quickly with N . "Quickly" ideally means exponentially, but high algebraic rates of decay are also suitable. This condition ensures an RBM emulator w N can achieve L 2 -proximity to the truth approximation w N when N/N is very small. We provide empirical evidence in this paper that this condition is true. Note that for our particular problem (36), exponential decay of the N width is known if we can assume that α is bounded away from 0 [17]. However Kato's formula requires α to take values arbitrarily close to 0. In a follow-up paper in preparation, we present N width estimates uniformly for α ∈ (0, 1] [3]. • The condition that the truth approximation w N (y) comes with a practically computable a posteriori error estimate ∆(y), satisfying, This usually comes in the form of a posteriori finite element estimates, and in practice in the algorithm are actually used to measure w N (y)−w N (y) . This condition ensures that the parameter values {y n } N n=1 in (37) can be chosen in a computationally tractable manner. In our case the PDE's we consider are linear so that efficient residual-based error indicators ∆ can be derived.
• The condition that the operator L and right-hand side f have affine dependence on the parameter y. This means that one has the expressions, where we have introduced (i) y-independent differential operators L q , (ii) yindependent functions f q (x), (iii) x-independent functions γ q (y), and (iv) xindependent functions σ q (y). More precisely, one requires the weak (variational) form of L to have such a decomposition. This condition is needed so that evaluation of the RBM emulator map y → u N (y) can be accomplished using operations that are independent of the truth approximation discretization parameter N . We will briefly justify this for our situation in the next section. Our next goal is to apply the RBM algorithm to the truth discretizations of w ± (y) that define the fractional solution u(s). We discuss this in the next section.

RBM formulation.
We describe here the RBM procedure for approximating w N − via a reduced basis emulator; the procedure for w N + is nearly identical. We recall the discrete truth approximation formulation that defines w N − : where S, M , and f are defined in (18). Detailed exposition of application of the RBM algorithm to this parameterized PDE (and to much more general cases) can already be found in existing textbook literature [37,26,39]. Here we give a only brief synopsis of the major steps in the algorithm for completeness, but refer to the previously-mentioned references for details and motivating explanation of the algorithm. In particular, in what follows we describe the algorithm using vectors and matrices instead of more common functional-analytic mathematical statements; this choice is made for simplicity of exposition since the algorithm itself is not new and we instead focus on the application of the algorithm. As indicated by relation (37), an RBM algorithm produces an emulator w n,− given by where we have made a particular choice of the basis φ k appearing in (37). 4 Also, since the RBM procedure builds w N sequentially by first building w 1 , w 2 , . . . , we label the RBM dimension as n, satisfying 1 ≤ n ≤ N , in this section. We must specify the parameter values {y k } n k=1 and the coefficients {c n,k } n k=1 , which is the focus of the following discussion.

5.3.
Computing the c n,k . We assume that y 1 , . . . , y n have been chosen and are known, and now seek to define the coefficients {c n,k } n k=1 , equivalently the vector c n , whose computation allow evaluation of y → w n (y). To proceed we define a new matrix U ∈ R N ×n , having entries U n,− = w N (y 1 ) w N (y 2 ) · · · w N (y n ) , (U n,− ) j,k = w N j (y k ), where the vector w N and its entries w N j are expansion coefficients for the truth approximation solution, see (17).
Then for each y, the coefficients c n,k of the RBM solution are defined by seeking the vector c n (y) ∈ R n satisfying Assuming U has linearly independent columns (which is assured by the choice of y k discussed in the next section), then this uniquely defines c n (y) for each y, and prescribes the RBM solution w n via (39). One final point of interest is that our truth variational form (38) exhibits affine dependence on the parameter y, making it possible to compute c n very efficiently. We may rearrange computations in (40) so that so that the quantities B, C, and g, once computed, are all independent of both the truth discretization dimension N and the parameter y. Thus, for each y, the coefficients c n (i.e., the RBM solution w n ) can be computed with complexity that depends only on n and not on N . Since in practice n N this can result in computational savings, especially if we wish to query y → w n − (y) numerous times. This is one of the major attractions of model reduction with RBM.

5.4.
Choosing y k . The ingredient we are left to provide to complete our description of the RBM algorithm is the choice of parameter values y k in (39). Given an RBM approximation w n , we focus on the choice of y n+1 . We accomplish this via the standard greedy procedure in RBM algorithms. Ideally, this choice is given by Unfortunately, this explicit form requires computing the full solution w N (y) at all parameter values y, which RBM seeks to avoid. To circumvent this restriction, the standard strategy is to resort to residual-based error indicators. The following lemma identifies one such computable residual-based error indicator ∆ n (y).
and the indicator Then ∆ n,± (y) ≥ w N ± (y) − w n,± (y) (44) Proof. We show the result for the "−" quantities; a similar proof works for the "+" quantities. The residual r n,− satisfies S + e −y M w N − (y) − w n,− (y) = r n,− , so that To summarize, we have the estimate which is the desired result by using the definition of C N in (19).
The residual vectors r n,± can be efficiently computed for many values of y. We illustrate this r n,− . We have: where P R(A) : R N → R N is the R N -orthogonal projector onto the column space of a matrix A, and The orthogonal decomposition in (45) shows that the Pythagorean theorem can be used to compute the Euclidean vector norm r n,− (y) 2 in an efficient way for several values of y: • P R(Rn) ⊥ f 2 is y-independent, so that it can be computed once and stored.
• P R(Rn) (f − (S + e −y M ) U n,− c n,− (y)) is a vector in a 2n-dimensional, yindependent vector space. Thus, the norm can be computed with only ndependent complexity. The fact that c n (y) appears with linear behavior in this expression ensures that we can rearrange computations so that, for each c n (y), the norm of this quantity can be computed using complexity that is dependent only on n.
In summary, while the right-hand side of (44) is not efficiently computable for many values of y, the left-hand side is efficiently computable for several values of y since r n 2 is efficient to compute, and λ min (S, M ) does not depend on y and can be computed either directly or iteratively with generalized eigenvalue solvers once and subsequently stored. Standard greedy algorithms for RBM methods require a computable quantity satisfying (44), and replace the essentially un-computable maximization (41) with the computable maximization y n+1 = argmax y≥0 ∆ n (y). (46) The above maximization has an objective function that is efficiently computable, and the inequality (44) ensures that the maximization (46) is a weak greedy algorithm. Weak greedy algorithms in turn ensure that the set of chosen parameters {y 1 , . . . , y N } defines an RBM subspace V N in (37) whose best approximation to the truth solution w N − is comparable to the Kolmogorov N -width [8,20]. We have completed the basic description of the RBM algorithm: the emulators w N,± (y) are defined by computing coefficients c N,± (y) as described in the previous section, and the parameter values y k are chosen according to (46) by computing the estimators ∆ n,± (y) for n = 1, 2, . . . ,. One usually sequentially computes y k until sup y≥0 ∆ n,± (y) is smaller than some specified tolerance, so that one can rigorously certify the error committed by the RBM emulators. This tolerance condition is usually how the terminal RBM dimension N is computationally specified.
One final observation we make is a major theoretical result of this paper: Proof. We have Taking L 2 norms of both sides, and using the triangle inequality with (44) yields the result.
We note that the quantities ∆ N±,± (y) are computed during the RBM construction phase, so that these estimators are available. Therefore, (47) provides a computable error bound that can be used to certify error committed by using the RBM algorithm.
Remark 4. Just as with Remark 4, all our results above extend to a more general elliptic operator E satisfying the assumptions outlined in Remark 4. In this case, all of our formulas in this section carry over. 5.5. Algorithm summary. The full algorithm of this section first uses the RBM algorithm to perform model reduction on the parameterized PDE solutions w N ± (y). 5 Subsequently, the efficient RBM emulators w N,± (y) are used in the GQ algorithm from Section 3.4. We describe the full algorithm in Algorithm 2.
The algorithm we have described in this section is a skeleton version of a modern RBM algorithm. We summarize various improvements that should be implemented in order for an RBM algorithm to be efficient and accurate: • One does not usually solve (46) by maximizing over the parameter continuum, and instead maximizes over a discrete set. For bounded parameter domains, it is common to use a uniform grid and subsequently adaptively (e.g., dyadically) refine the grid to ensure that no local maxima are skipped. Over the unbounded domain, we employ a logarithmic map; see section 6.4 for details. • Our RBM ansatz (39) uses solution snapshots w N ± (y) as basis functions. This is known to generally lead to ill-conditioning in the formulation (40) even for small n (since in practice the columns of U n,± are "nearly" linearly dependent). A better prescription is to build the RBM basis functions by orthogonalizing snapshots.
• Some naive implementations of the decomposition (45) via quadratic forms leads to numerical roundoff error that results in stagnation of the error indicators ∆ n,± (y) near root-machine precision. (E.g., for double precision, stagnation occurs when ∆ n,± (y) takes values around 10 −8 . More careful computations allow one to overcome this limitation [14,15,12,16]. We again refer to [37,26,39] for a more complete description of important but standard RBM algorithm details.
Finally, we note that one can combine the error estimates in (47) and (31) via the triangle inequality to create a computable bound for u N − u N . This error would estimate the error committed by the two novel innovations of this paper: our Gaussian quadrature approach and the model reduction procedure. 5 In the previous sections we have describe this process only for w N − , but the process for w N + results in almost the same procedure with only minor differences stemming from the location of the e −y factor. Algorithm 2 RBM algorithm: Produces solution to the fractional Laplace problem (5). The function OfflineFracLapRBM needs to be completed only once and has complexity dependent on N = dim V . Afterwards, OnlineFracLapRBM can be called for arbitrary values of s ∈ (0, 1) with a computational cost dependent only on dim V N = N N .

6:
Compute and store w N − (y n ).

11:
Call FracLapGQ(s), Algorithm 1, replacing w N ± (·) with RBM emulators w N±,± (·). 12: return u N (s). 6. Numerical examples. In these experiments, we compare the effectiveness of our improvements to current methods. We first describe the setup of the test problems that we consider in our simulations. We solve (5) on Ω = [0, 1] 2 with homogeneous Dirichlet boundary conditions. We test a total of 3 algorithms: • "SQ" -The sinc quadrature approach from [11].
• "GQ" -Algorithm 1 detailed in section 3, utilizing a modified version of the integral formulation in [11] along with Gaussian quadrature. • "RBM" -The approach detailed in section 5.5, leveraging the reduced basis method to accelerate the GQ algorithm.
To compare the three methods, we will use the same number of y-quadrature points M = M − + M + in each approach. This number represents the total number of local PDE solves needed to compute an approximation to u(s). In particular, we make the choices given in [11,6], which depend on the spatial mesh and fractional order: The finite element discretization is accomplished with linear quadrilateral finite elements on a Cartesian tessellation of Ω. The one-dimensional grids that define this Cartesian tessellation are isotropic with respect to the two dimensions, and are defined as equidistant meshes with 2 K points. We will use various values of K.
With this in hand, and using the inverse of the relation (4), we can easily construct explicit solutions for testing using eigenfunction expansions. We explore the effectiveness of our algorithms through three manufactured solutions: • "Sine" -The function u and data f are, in this case, f (x) = sin(πx 1 ) sin(πx 2 ) • "Mixed modes" -The function u and data f are, in this case, • "Square bump" -The data f is the indicator function While an analytical solution is available as an infinite sum of eigenfunctions, we instead numerically compute a solution on a grid formed by 2 9 + 1 equally spaced nodes in both directions and consider this the "exact" solution. In comparison, the finest mesh used to evaluate the accuracy of our solver has 2 6 + 1 nodes in both directions.
6.2. Spatial convergence. Our first test verifies that we recover spatial convergence in terms of the finite element mesh size. Since we are primarily interested in accuracy and not efficiency, this section compares the SQ and GQ methods, with the results summarized in Figure 3.
We can see that the proposed GQ algorithm performs slightly better than the existing SQ algorithm for the same number of quadrature points M . We also remark that the GQ implementation allows us to generate solutions for small fractional parameters with less numerical difficulty. With the SQ approach, the difficulty arises when application of the quadrature rule results in large values of a term involving e y appearing in operators that must be inverted. 6.3. Quadrature rule efficiency. In this section we compute errors committed by the GQ and SQ algorithms for different values of the quadrature rule size M . The purpose of this test is to understand the efficiency of the quadrature rule, i.e., the number of solutions of w N ± required. Figure 4 illustrates errors for the three test cases as a function of the total number of quadrature nodes. The test in this section does not fix M ± as given by (48). Instead, given a number of quadrature points M (the abscissa in Figure 4), we generate M -point quadrature rules for the Gaussian quadrature and sinc approaches. Thus, the quadrature rule for each M is generated anew.
The results indicate that the GQ algorithm converges faster than the SQ with respect to the number of PDE solves. This shows that the GQ algorithm appears to be far more efficient than sinc quadrature for computing solutions to these fractional problems. Figure 3. Convergence of SQ (solid blue) and GQ methods (red dashed) as the spatial mesh is refined. Each used a dyadic mesh along the spatial variable with increasing resolution. A parameter value of s = 0.2 was used and similar results were seen for value of s between 0.1 and 0.9. We used the number of quadrature points for the integral suggested by current methods [6]. 6.4. RBM offline efficiency. This section investigates the RBM algorithm. For now, we restrict our attention to one-time queries of u(s), i.e., to situations when, given Ω, f , and s, we seek to compute only u(s) for this given s. For the GQ algorithm this involves a single run of the routine FracLapGQ in Algorithm 1. For the RBM algorithm, this entails a single run of the OfflineFracLapRBM routine in Algorithm 2, followed by a single run of OnlineFracLapRB routine.
For the RBM algorithm, we solve (46) by discretizing the y ∈ [0, ∞) domain in a uniform way under a logarithmic map. Precisely: we set z = e −y for y ∈ [0, ∞) and proceed to discretize z ∈ (0, 1]. We take 128 equispaced points in the z variable and map back to y-space with z → − log z = y. We subsequently perform a discrete maximization over this set instead of the continuous optimization (46).
In figure 5 we compare the GQ algorithm to the accelerated RBM algorithm, including the offline construction time. We see that the initial investment of the RBM algorithm in the offline phase is substantial, accumulating to the time required for the direct GQ method with 150-200 quadrature points. However, we see that after this initial offline investment, subsequent evaluations of the RBM surrogate are extremely efficient, so that the effort required to evaluate M 1 quadrature point is essentially the same as that required to evaluate at a single quadrature point. Figure 5. "Offline" (i.e., one-time) computational investment for a single solve of (5) with a fixed value of s = 0.2. These experiments compare both the direct (dotted and dashed) and reduced basis methods (solid) using the gaussian quadrature. Each used a dyadic mesh with 7 levels. Similar results where seen for value of s between 0.1 and 0.9. 6.5. RBM accuracy. We now investigate the accuracy delivered by the RBM algorithm in the construction of reduced order models for w N ± . Our rigorous error certificate for u(s) using the reduced order model is (47), but we here consider a finer estimate using the proof of Theorem 5.2. We define the error estimator, where ∆ N,± are the computable error indicators defined in (43), and we again choose M ± as in (48). One can see from the proof of Theorem 5.2 that this quantity bounds the error committed by the RBM procedure. We plot ∆ N in Figure 6 as a function of N , and observe that it decays exponentially.
Finally, we remark that ∆ N only certifies the error committed by the model reduction RBM algorithm; the error committed by the y-quadrature rule is not certified by this quantity. 6.6. RBM accuracy and efficiency. Finally, we explore the accuracy afforded by the RBM procedure for computing s → u(s), and also verify the computational efficiency of the procedure. In Figure 7 left and center, we demonstrate that the error committed by the RBM algorithm is stable, even for relatively small values of the parameter 0.01 < s < 0.1. Furthermore, the right pane of this figure demonstrates that if we wish to repeatedly query the map s → u(s) for several values of s, the RBM algorithm is undeniably more efficient by an order of magnitude even for just one query, and by three orders of magnitude if 1000 queries are needed. This is an especially salient point for optimization or control applications involving forward models with fractional PDEs of varying fractional order s, where repeated evaluation of such PDE solutions is needed. Figure 7. Accuracy of the RBM algorithm over a range of values of s (left and center). The Sine example is plotted in a red dot-dashed, the Mixed Modes in a blue solid line, and the Square Bump case in black crosses. In the right pane we show the cumulative computational time required by the GQ algorithm (blue) versus the RBM algorithm (red). Each query refers to an evaluation of the map s → u(s). In particular this cumulative time for the RBM solver includes the one-time offline cost required by OfflineFracLapRBM in Algorithm 2.

Conclusion.
We propose a novel model reduction strategy for computing solutions to fractional Laplace PDE's, in particular (5). Our algorithm builds on the ideas introduced in [11], improving accuracy and stability, and accelerating that algorithm considerably. Our model reduction strategy hinges on the fact that the solution to the fractional problem can be written in terms of classical, local elliptic PDE's, for which RBM-based model reduction is known to be efficient.
We provide novel stability bounds for both the continuous and discrete problems, and our numerical experiments suggest that our Gaussian quadrature approach is more efficient than alternative quadrature methods. All of our algorithmic and theoretical results apply to solutions to differential equations involving fractional powers of general elliptic operators. A rigorous proof of the convergence for our quadrature rule is the subject of ongoing study.