Model reduction for fractional elliptic problems using Kato's formula

TL;DR

This paper introduces a stable, efficient model reduction algorithm for solving fractional elliptic PDEs using Kato's formula, significantly reducing computational costs for multiple queries.

Contribution

It develops a novel, stable discretization of the integral representation of fractional Laplacian solutions combined with reduced basis methods for efficient multi-query computations.

Findings

Stable discretization independent of discretization parameters

Gaussian quadrature improves efficiency over previous methods

Order of magnitude reduction in computational cost for multiple queries

Abstract

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.

Equations

$$(-\Delta {)}^{s}u=f,$$

$$\Vert \mathbf{x}{\Vert }_{\mathbf{M}}^2$$

$${\lambda }_{\min }\left(\mathbf{A},\mathbf{B}\right):=\underset{\mathbf{x}\in \text{mathbbm}{R}^N\backslash \{0\}}{\inf }\frac{\Vert \mathbf{x}{\Vert }_{\mathbf{A}}^2}{\Vert \mathbf{x}{\Vert }_{\mathbf{B}}^2}.$$

$${\lambda }_{\min }\left({\mathbf{B}}^{-1/2}\mathbf{A}{\mathbf{B}}^{-1/2}\right)={\lambda }_{\min }\left(\mathbf{A},\mathbf{B}\right),$$

$$L^2(\Omega )$$

$$\begin{array}{rl}-\Delta u=\lambda u,& x\in \Omega ,\\ u(x)=0,& x\in \partial \Omega ,\end{array}$$

$$u\in L^2(\Omega )⟹u(x)=\sum _{n=1}^{\infty }{u}_n{\varphi }_n(x),u_n$$

$$u$$

$${\mathbb{H}}^s(\Omega )$$

$$\begin{array}{rl}(-\Delta {)}^{s}u=f,& x\in \Omega \\ u=0,& x\in \partial \Omega \end{array}$$

$$\displaystyle U\coloneqq\left\{u(s)\;\big{|}\;s\in(0,1)\right\}.$$

$$u(s)$$

$$\begin{array}{rl}-\Delta q(y)+e^{y}q(y)=f,& x\in \Omega \\ q(y)=0,& x\in \partial \Omega .\end{array}$$

$$u(s)=\beta (s){\int }_{-\infty }^{\infty }q(y)e^{(1-s)y}\mathrm{d}y.$$

$$u(s)\approx u_M(s):=\sum _{m=1}^M{\tau }_{m}q(y_m)e^{(1-s)y_m}.$$

$$\begin{array}{rl}-\alpha \Delta w(\alpha ,\beta )+\beta w(\alpha ,\beta )=f,& x\in \Omega \\ w(\alpha ,\beta )=0,& x\in \partial \Omega .\end{array}$$

$$w_{-}(y)$$

$$(-\Delta +e^{-y})w_{-}(y)$$

$$u(s)$$

$$={\beta }_0(s_{-}){\int }_0^{\infty }{w}_{-}\left(\frac{y}{s_{-}}\right)W(y)\mathrm{d}y+{\beta }_0(s_{+}){\int }_0^{\infty }{w}_{+}\left(\frac{y}{s_{+}}\right)W(y)\mathrm{d}y,$$

$${\beta }_0(s)$$

$$u$$

$$=(1-s){\beta }_0(1-s){\int }_0^{\infty }{w}_{-}(y)\exp (-(1-s)y)\mathrm{d}y+s{\beta }_0(s){\int }_0^{\infty }{w}_{+}(y)\exp (-sy)\mathrm{d}y$$

$$={\beta }_0(1-s){\int }_0^{\infty }{w}_{-}\left(\frac{y}{1-s}\right)W(y)\mathrm{d}y+{\beta }_0(s){\int }_0^{\infty }{w}_{+}\left(\frac{y}{s}\right)W(y)\mathrm{d}y,$$

$$\Vert v\Vert$$

$${\tilde{C}}_{\Omega }^2:=\max \left(1,C_{\Omega }^2\right),$$

$$a\left(w,v;\alpha ,\beta \right)$$

$$a\left(w,v;\alpha ,\beta \right):=\alpha ⟨\nabla w,\nabla v⟩+\beta ⟨w,v⟩.$$

$$\underset{s\in (0,1)}{\sup }\Vert u(s)\Vert \le \frac{4{\tilde{C}}_{\Omega }^2}{\pi }\Vert f\Vert .$$

$$a\left(w_{+}(y),v;e^{-y},1\right)$$

$$\Vert f\Vert \Vert w_{+}(y)\Vert$$

$$\ge \left(1+\frac{e^{-y}}{c_{\Omega }^2}\right)\Vert w_{+}(y){\Vert }^2\ge \Vert w_{+}(y){\Vert }^2.$$

$$\underset{y\ge 0}{\sup }\Vert w_{+}(y)\Vert \le \Vert f\Vert .$$

Full text

H. Antil was partially supported by National Science Foundation awards DMS-1818772 and DMS-1521590 and Air Force Office of Scientific Research (AFOSR) under Award NO: FA9550-19-1-0036. Y. Chen’s research was supported by NSF awards DMS-1719698 and AFOSR award FA9550-18-1-0383. E. Cherkaev was supported by NSF awards DMS-0940249 and DMS-1413454. A. Narayan was partially supported by NSF award DMS-1720416, AFOSR award FA9550-15-1-0467, and DARPA EQUiPS contract N660011524053.

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 hysteric and globally-coupled phenomena. For example, fractional differential equations have recently been used to model fluid mechanics, arterial blood flow, cardiac ischemia, and are been used as ingredients in image denoising and image segmentation [36, 30, 10, 24, 2, 3]. Fractional PDEs have shown tremendous potential to model applications in geophysics [40] and manifold learning [4]. Finally, we mention the novel optimal control concepts introduced by fractional equations [7], see also [13, 11, 12]. In this paper we focus on fractional elliptic operators, defined via spectral expansion; a prototypical example of an operator that we use throughout this paper is the fractional Laplacian.

With $\Delta$ the classical Laplacian on a physical domain $\Omega$, we are interested in computing the solution $u$ to the partial differential equation

[TABLE]

with appropriate boundary conditions on $\partial\Omega$, where $s\in(0,1)$ is the fractional order. The precise definition of the fractional operator $(-\Delta)^{s}$ involves the spectral expansion of the classical operator $-\Delta$, which we more formally 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 $(-\Delta)^{s}$ to devise a scheme that computes solutions using the spectral expansion of the associated discretized operators [27, 28, 43, 38]. The disadvantage of this approach is that the procedure is expensive, requiring a full eigendecomposition of a potentially very large matrix. In addition, it is difficult to certify error using this approach.

• •

A second approach uses an extension procedure to write the non-classical $d$-dimensional PDE as a $(d+1)$-dimensional classical PDE [33, 18, 39]. This latter PDE can be solved with existing methods, although some nontrivial tailoring of existing numerical methods is needed [34, 32, 1]. The challenge with this approach is that the spatial dimension is increased, and the extended PDE is degenerate, requiring specialized numerical methods.

• •

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 as a type of Dunford-Taylor integral involving the resolvent of the classical operator see [29, Theorem 2 with simplification $\lambda=0$], see also [16, 8]. This approach discretizes the integral formulation with quadrature and requires several classical PDE solves (equal to the number of quadrature points) in order to compute the solution to the fractional problem. However, this results in an algorithm that can use several queries of an existing PDE solver to compute the solution to the fractional problem. The challenge with this approach is that $\mathcal{O}(100)$ classical 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 classical PDE solves may be accomplished in parallel, but this does not diminish the overall cost. Other operators require coupling of these solves, making parallel approaches more difficult [41].

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 [15], our approach would extend to more general cases as well. Our contributions in this article are as follows:

• •

We provide a rearrangement of the Dunford-Taylor integral considered in [16] 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.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.1. 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 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)\mapsto u$ by at least two orders of magnitude. The offline portion of this algorithm requires approximately as much time as a single $(s,f)\mapsto 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.1.

We remark that while we study PDEs with an operator of the form $(-\Delta)^{s}$, all our results extend to more general fractional elliptic operators. See Remarks 3.1 and 5.1.

This paper is not the first strategy for model reduction for fractional elliptic problems. The authors in [5] provide a model reduction strategy, applied to the second (extension) approach listed above. More recently, the work in [22] employs a reduced basis approach by interpolating operator norms. However, low-rank structure in solution sets to fractional problems has been empirically noted even earlier [42]. 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.

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 [16]. Section 5 utilizes RBM to propose a new model reduction algorithm that computationally accelerates the algorithm from Section 3 and provides a computable error certificate for the model reduction. Finally, section 6 demonstrates our new algorithms on a two-dimensional fractional Laplace problem and compares our algorithm against the predecessor in [16].

2. Notation and setup

Vectors will be denoted in lowercase bold, and matrices in uppercase bold, e.g., $\boldsymbol{x}$ and $\boldsymbol{A}$, respectively. If $\boldsymbol{M}$ is a symmetric positive definite matrix, we define

[TABLE]

and $\|\boldsymbol{x}\|$ is the standard Euclidean norm. The matrix norm $\|\boldsymbol{A}\|$ is the standard induced $\ell^{2}$ norm on matrices. If $\boldsymbol{A}$ is symmetric, then $\lambda_{\mathrm{\min}}\left(\boldsymbol{A}\right)$ denotes the smallest (real) eigenvalue of $\boldsymbol{A}$. If both $\boldsymbol{A}$ and $\boldsymbol{B}$ are symmetric positive definite matrices in $\mathbbm{R}^{N\times N}$, we define the smallest generalized eigenvalue of $(\boldsymbol{A},\boldsymbol{B})$ as

[TABLE]

Note that under these assumptions on $\boldsymbol{A}$ and $\boldsymbol{B}$, the above expression is equal to the smallest $\lambda$ such that $\boldsymbol{A}\boldsymbol{x}=\lambda\boldsymbol{B}\boldsymbol{x}$ has a nontrivial solution $\boldsymbol{x}$, and also we have that

[TABLE]

where $\boldsymbol{B}^{1/2}$ is the symmetric positive definite matrix square root of $\boldsymbol{B}$. Similarly, we use the notation $\lambda_{\mathrm{max}}(\cdot)$ and $\lambda_{\mathrm{max}}(\cdot,\cdot)$ to denote maximum eigenvalues.

Consider a bounded domain $\Omega\subset\mathbbm{R}^{d}$ with Lipschitz boundary $\partial\Omega$; we are mainly concerned with $d\leq 3$. We have

[TABLE]

We will often write $L^{2}=L^{2}(\Omega)$, and we define $\langle\nabla w,\nabla v\rangle=\sum_{j=1}^{d}\langle\frac{\partial}{\partial x_{j}}w,\frac{\partial}{\partial x_{j}}v\rangle$, which induces a definition for the $L^{2}$ norm $\|\nabla v\|$ of vector-valued functions. For brevity will also write $\|v\|=\|v\|_{L^{2}}$ and $\|\nabla v\|=\|\nabla v\|_{L^{2}}$. The standard Laplace eigenvalue problem on $\Omega$ with Dirichlet boundary conditions,

[TABLE]

yields an infinite sequence of eigenvalues $0<\lambda_{1}\leq\lambda_{2}\cdots$ with associated eigenfunctions $\{\phi_{n}\}_{n=1}^{\infty}$. Here, and in all the following, the differential operator $\Delta$ operates on the $x$ variable. The spectral theorem ensures that the eigenfunctions enjoy $L^{2}(\Omega)$-orthogonality and completeness, so that

[TABLE]

where we have further assumed that each $\phi_{n}$ has unit $L^{2}(\Omega)$ norm.

2.1. The fractional Laplace problem

Let $s\in(0,1)$. In this section we describe the spectral definition of the fractional operator $(-\Delta)^{s}$, on bounded domains, supplemented with homogeneous Dirichlet boundary conditions, see [9] 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,

[TABLE]

We likewise use $\lambda_{n}$ to define Sobolev spaces of fractional order. With $s\in(0,1)$:

[TABLE]

where $\mathbb{H}^{s}(\Omega)^{\ast}$ denotes the dual space of $\mathbb{H}^{s}(\Omega)$. With these definitions, $(-\Delta)^{s}:\mathbb{H}^{s}(\Omega)\rightarrow\mathbb{H}^{-s}(\Omega)$. For the relation of $\mathbb{H}^{s}(\Omega)$ to the standard fractional order Sobolev space, see [9].

Given data $f\in L^{2}(\Omega)$, our main goal is to compute the solution $u$ to

[TABLE]

for arbitrary $s\in(0,1)$. 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)\in L^{2}$ denote the solution $u$ to (4) for a fixed value of $s$. We will be interested in developing a computational algorithm for computing the family or manifold of solutions,

[TABLE]

If $f\in L^{2}$, then the solution $u$ to (4) has $\mathbb{H}^{2s}$ membership, so that the natural function space in which to study the manifold $U$ is $\cap_{s\in(0,1)}\mathbb{H}^{2s}$. Hence, all of our investigations will assume $f\in L^{2}$ and study solutions $u(s)$ as elements of $L^{2}$. Note that $f\in L^{2}$ is a stronger requirement than for classical elliptic problems.

In the remainder of this document we will consider the problem (4) with homogeneous Dirichlet boundary conditions. The inhomogeneous boundary case may be handled using the lifting technique in [9], which requires only a small modification of the algorithm proposed here.

2.2. Kato’s integral solution of (4)

The following remarkable result provides an appealing formula for the solution $u$ to (4):

[TABLE]

which is a reformulation of Kato’s formula [29, Theorem 2 with simplification $\lambda=0$]. This representation was first exploited in [16] 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\in\mathbbm{R}$, as the solution to the classical $y$-parameterized PDE,

[TABLE]

This representation reveals that $u$ is actually just an integral of solutions $q$ to classical 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 (4). The straightforward way to compute the solution via (5b) is to approximate the integral with a quadrature rule. This would require computing solutions $q(\cdot;y)$ to the PDE (5a) for many values of the parameter $y$.

More precisely, let $\{y_{m},\tau_{m}\}_{m=1}^{M}$ be a quadrature rule for approximating the integral in (5b). We will give precise choices for this quadrature rule soon. Then we can approximate the solution $u(s)$ as

[TABLE]

One then needs only to compute the ensemble of functions $\{q(y_{m})\}_{m=1}^{M}$, which are solutions to classical PDE’s, in order to approximate the solution to the fractional problem. This is the particular approach adopted in [16], 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 a classical problem; when $M$ is large (which can be required when $s$ is small), this may become computationally prohibitive, see for instance [40]. 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 (5) 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 [16].

3. Fractional Laplace solutions via integral formulation

Recall that given $s\in(0,1)$, we seek to evaluate (5b), which defines the solution $u(s)$ to the fractional Laplace problem (4). 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 [16]; 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 [16] 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 [16] with a Gaussian quadrature rule.111The authors in [16] also propose a “Gaussian” quadrature rule, 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.

3.1. Partitioning of the $y$ variable

To make computations numerically stable, we split our parameterized problem (5a) into regions $y\in(-\infty,0]$ and $y\in[0,\infty)$. To accomplish this introduce a new parameterized PDE for a solution $w$ that is closely related to the solution $q$ from (5a)

[TABLE]

where $(\alpha,\beta)\in(0,\infty)\times[0,\infty)$ is the parameter. In our computational setting, we will only require $(\alpha,\beta)\in(0,1]^{2}$. Comparing (6) with (5a), we see that $q(y)=w(1,e^{y})$. We now define $w_{\pm}(y)$ as two specializations of $w$ that will be used in the following:

[TABLE]

for $x\in\Omega$, with boundary conditions $w_{\pm}(y)=0$ for $x\in\partial\Omega$. We can now formulate the solution to the fractional PDE (4) in terms of these new quantities.

Lemma 3.1**.**

The solution $u(s)$ to (4) is given by

[TABLE]

where $\beta_{0}$, $s_{\pm}$, and $W$ are defined as

[TABLE]

Proof.

Beginning with (5b), we have

[TABLE]

completing the proof. ∎

Given the domain $\Omega$, we define $C_{\Omega}$ as the domain’s Poincaré constant, i.e., the smallest constant such that for every $v\in H_{0}^{1}(\Omega)$,

[TABLE]

The space $H_{0}^{1}(\Omega)$ is the standard Sobolev space of zero-trace $L^{2}(\Omega)$ functions whose gradients are also in $L^{2}$. In what follows, we will also need the following quantity,

[TABLE]

which also depends only on $\Omega$. Fixing $(\alpha,\beta)$, the weak formulation of (6) seeks a solution $w(\alpha,\beta)=w\in H^{1}_{0}\left(\Omega\right)$ as the unique function satisfying the Galerkin formulation,

[TABLE]

with the bilinear form $a(\cdot,\cdot;\alpha,\beta)$ defined as

[TABLE]

With $\alpha>0$ and $\beta\geq 0$, the Poincarè inequality (10) ensures that the coercivity property $a(v,v;y)\geq\alpha\|\nabla v\|^{2}_{L^{2}(\Omega)}\geq k\|v\|_{H^{1}_{0}(\Omega)}$ holds for some $k>0$ uniformly in $\beta$, so that standard Lax-Milgram theory then yields a unique $H_{0}^{1}(\Omega)$ solution.

With this setup, we can demonstrate the utility of a formula like (8) by deriving an $s$-independent $L^{2}$ stability estimate for solutions to (4).

Proposition 3.1**.**

Assume $f\in L^{2}$. Then

[TABLE]

Proof.

The function $w_{+}(y)\in H_{0}^{1}(\Omega)$ is the unique solution to

[TABLE]

Taking $v=w_{+}(y)$ and using the Cauchy-Schwarz and Poincaré inequalities results in

[TABLE]

We thus obtain

[TABLE]

Therefore, taking the $L^{2}$ norm in (8) and using the triangle inequality yields

[TABLE]

where the third inequality uses the fact that $W$ is a probability density on $[0,\infty)$. From the above, we immediately obtain the desired result by noting that

[TABLE]

The proof is complete. ∎

3.2. Spatial discretization

In this section we employ a spatial discretization to the result of Lemma 3.1. We proceed to discretize (13) using a finite element method. Let $T_{\Omega}$ be a conforming triangulation of $\Omega$ 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\geq 1$, we define the finite element space

[TABLE]

where $P_{k}(e)$ is the space of polynomials of degree $k$ or less over the element $e\in T_{\Omega}$. Let $\mathcal{N}=\dim V$. The finite element-discretized version of (12) is the Galerkin formulation seeking $w^{\mathcal{N}}\in V$ satisfying

[TABLE]

Let $w^{\mathcal{N}}(\alpha,\beta)$ be expressed as a linear expansion,

[TABLE]

where $\{\psi_{n}\}_{n=1}^{\mathcal{N}}$ is a basis for $V$, e.g., a basis comprised of compactly supported piecewise polynomials. Collecting the linear degrees of freedom of $w^{\mathcal{N}}(y)\in V$ in the $\mathcal{N}$-dimensional vector $\boldsymbol{w}^{\mathcal{N}}$, then this vector satisfies the linear system,

[TABLE]

for $j,k=1,\ldots,\mathcal{N}$. Above we have defined the $\mathcal{N}\times\mathcal{N}$ $y$-independent stiffness and mass matrices $\boldsymbol{S}$ and $\boldsymbol{M}$, respectively. Both $\boldsymbol{S}$ and $\boldsymbol{M}$ are symmetric and positive-definite.

The matrices $\boldsymbol{S}$ and $\boldsymbol{M}$ can be used to define a “discretized” Poincaré constant $C_{\mathcal{N}}$ by using a standard Rayleigh quotient argument:

[TABLE]

hence leading to the inequalities,

[TABLE]

Our estimates for $x$-discrete quantities will involve $C_{\mathcal{N}}$, but we will sometimes use the above inequality to bound quantities in terms of $C_{\Omega}$. Bounds involving $C_{\Omega}$ emphasize independence of the $x$-discretization. Bounds involving $C_{\mathcal{N}}$ emphasize the explicit computability of the bounds, since $C_{\mathcal{N}}$ is equal to an extremal eigenvalue of finite element matrices, which is computable with iterative eigenvalue solvers.

The maximum generalized eigenvalue of $(\boldsymbol{S},\boldsymbol{M})$ will also play a small role in our estimates. In analogy with (18) we define

[TABLE]

Note that $K_{\mathcal{N}}$ in general tends to 0 as $\mathcal{N}\uparrow\infty$.

From the discretization of $w$ we derive discretizations of $w_{\pm}(y)$ defined in (7). Thus, finite element discretizations for $w_{\pm}$ are specializations of those for $w$. In particular, we define

[TABLE]

Denote by $\boldsymbol{w}_{\pm}^{\mathcal{N}}(y)$ the $\mathcal{N}$-dimensional vectors that are solutions to the linear systems,

[TABLE]

so that, akin to (16), we have

[TABLE]

We can now codify the fact that the solutions $w^{\mathcal{N}}_{\pm}(y)$ are $L^{2}$-stable uniformly in $y$.

Lemma 3.2**.**

Assume $f\in L^{2}(\Omega)$. Then

[TABLE]

Proof.

The result can be obtained by considering the discrete form (21). To begin we relate the $L^{2}$ norm of $f$ to the Euclidean $\ell^{2}$ norm of $\boldsymbol{f}$. Let $P_{V}$ denote the $L^{2}$-orthogonal projector onto $V$. Then:

[TABLE]

Thus we have

[TABLE]

Now since $\boldsymbol{M}$ is symmetric and positive-definite, it has a unique symmetric positive-definite square root $\boldsymbol{M}^{1/2}$. Thus:

[TABLE]

This in turn implies:

[TABLE]

Since $\boldsymbol{M}^{-1/2}\boldsymbol{S}\boldsymbol{M}^{-1/2}$ is symmetric and positive-definite, we have

[TABLE]

Therefore,

[TABLE]

which, when combined with (23) yields (22a). A similar computation for $w_{-}(y)$ yields (22b) by using the definition of $C_{\mathcal{N}}$ in (18). ∎

The result above gives the stability of an algorithm that uses $w^{\mathcal{N}}_{\pm}(y)$ as a spatial discretization. In particular, consider the following semi-discrete approximation of $u(s)$,

[TABLE]

A fully discrete scheme, introduced in the next section, would discretize the $y$ variable. The following result mirrors the stability estimate of Proposition 3.1, showing $s$-uniform $L^{2}$ stability of the semi-discrete solution.

Proposition 3.2**.**

Assume $f\in L^{2}$. Then

[TABLE]

The proof, which we omit, is almost exactly the same as for Proposition 3.1 using the discrete stability estimates (22); the only novelty is the need to exercise the inequality (18).

Having described the spatial discretization, we now proceed to describe a discretization for the $y$-integrals in (8).

3.3. Quadrature for the $y$-integral

We will use an $M$-point $W$-Gaussian quadrature rule to discretize the integrals in (8). The weight function $W(y)$ is a weight function associated to a classical family of orthogonal polynomials: Laguerre polynomials. Let $\{p_{n}\}_{n\geq 0}$ denote the family of Laguerre polynomials, orthonormal under the weight $W$, i.e.,

[TABLE]

where $\delta_{m,n}$ is the Kronecker delta function. Like all orthogonal polynomials, the family $\{p_{n}\}_{n\geq 0}$ satisfies a three-term recurrence formula,

[TABLE]

with $p_{0}\equiv 1$ and $p_{-1}\equiv 0$. The recurrence coefficients $(a_{n},b_{n})$ are explicitly known:

[TABLE]

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},\tau_{m})_{m=1}^{M}$, respectively, and integrates polynomials up to degree $2M-1$ exactly:

[TABLE]

Although $y_{j}$ and $\tau_{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\times M$ symmetric tridiagonal Jacobi matrix,

[TABLE]

consider its associated eigenvalue decomposition,

[TABLE]

where $v_{j,1}^{2}$ is the first component of the vector $\boldsymbol{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 $\sum_{j=1}^{M}\tau_{j}=1$, and furthermore each $\tau_{j}$ is non-negative by (25).

Now let $M_{-}$ and $M_{+}$ be the number of quadrature points used to approximate the “$-$” and “$+$” integrals in (8), respectively. This results in two sets of $W$-Gaussian quadrature rules,

[TABLE]

We apply these rules to the integrals in (24), resulting in the fully discrete approximation,

[TABLE]

We emphasize that the $y$-discretization does not suffer from any numerical instabilities as $s\uparrow 1$ or $s\downarrow 0$: the weights $\tau_{j,\pm}$ are positive, no larger than 1, and independent of $s$, $\beta_{0}(1-s)$ and $\beta_{0}(s)$ are just sinc functions, and $w_{\pm}(y)$ has bounded $L^{2}$ norm for all $y\geq 0$, i.e., for all inputs. The following codifies this stability.

Proposition 3.3**.**

Assume $f\in L^{2}(\Omega)$. Then

[TABLE]

Proof.

The proof is very similar to the proof for proposition (3.1). Take the $\|\cdot\|_{\boldsymbol{M}}$ norm on both sides of (26), use the triangle inequality and (22), and note that the quadrature weights $\tau_{j,\pm}$ all satisfy $0\leq\tau_{j,\pm}\leq 1$ since they are all positive and $W$ is a probability density. Note that (27) 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 [16], assuming we have the ability to compute $\boldsymbol{w}^{\mathcal{N}}_{\pm}$, then the formulation above immediately yields an algorithm to compute $u^{\mathcal{N}}(s)$ via (26).

Remark 3.1*.*

All of our results extend to the case when we solve (4) but replacing $-\Delta$ with a general elliptic operator $\mathcal{E}$ that (i) satisfies the coercivity condition $\left\langle\mathcal{E}v,v\right\rangle\geq\alpha\|v\|_{H_{0}^{1}}$ for $\alpha>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 $\widetilde{C}^{2}_{\Omega}$ with $\widetilde{C}^{2}_{\Omega}/\alpha$. The matrix $\boldsymbol{S}$ should likewise be replaced with the associated matrix defined from the bilinear form of $\mathcal{E}$.

3.4. Algorithm summary

The sections above identify an algorithm for computing $u^{\mathcal{N}}(s)$ in (26). We summarize the procedure in Algorithm 1. The discrete solution adheres to the stability bound in Proposition 3.3. Note that we have not yet described how one should decide on values for $M_{\pm}$. 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_{\pm}$ using (46).

4. Error due to quadrature discretization

The formulation (26) is our numerical approximation to compute solutions to (4). 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 $\widetilde{u}^{\mathcal{N}}(s)$ and $u^{\mathcal{N}}(s)$.

To proceed, we need an auxiliary function that measures the absolute error between a size-$M$ quadrature rule $\left(y_{j},\tau_{j}\right)_{j=1}^{M}$ and the exact integral applied to a particular function:

[TABLE]

We also need to define intervals on the real line enclosing the spectrum of some discretized operators. Recalling the definitions of $C_{\mathcal{N}}$ and $K_{\mathcal{N}}$ in (18) and (20), respectively, define intervals $I$ and $I_{\pm}$ as

[TABLE]

The error committed by the quadrature rule can be understood by studying the quantity,

[TABLE]

The precise statement is as follows.

Proposition 4.1**.**

Assume $f\in L^{2}$. Then for each $s\in(0,1)$,

[TABLE]

and therefore,

[TABLE]

Proof.

The same argument that produces the relations (23) implies that

[TABLE]

so we proceed to study the quantity on the right-hand side. The difference between the “$-$” integral contributions in $\boldsymbol{u}^{\mathcal{N}}(s)-\boldsymbol{\widetilde{u}}^{\mathcal{N}}(s)$ is proportional to

[TABLE]

We can express the solution $\boldsymbol{w}_{-}^{\mathcal{N}}(y)$ as

[TABLE]

where we have defined $\boldsymbol{A}\coloneqq\boldsymbol{M}^{-1/2}\boldsymbol{S}\boldsymbol{M}^{-1/2}$. The matrix $\boldsymbol{A}$ is symmetric and positive definite, and thus has an eigenvalue decomposition

[TABLE]

Then further manipulation of the $\boldsymbol{w}_{-}^{\mathcal{N}}(y)$ expression yields

[TABLE]

where $\boldsymbol{H}(y)$ is a diagonal matrix having entries

[TABLE]

Thus,

[TABLE]

The first inequality is sub-multiplicativity of the $\|\cdot\|$ matrix norm, and the first equality uses the invariance of the same norm under unitary transformations. A similar computation for the “$+$” quantities yields

[TABLE]

The combination of these results implies

[TABLE]

Using the inequality (19) yields the result. ∎

The summation on the right-hand side of (30) can be computed independent of the data $f$, and requires only knowledge of the extremal eigenvalues of the discrete operator, cf. (18) and (20). While we cannot at present provide a theoretical estimate of this error, we numerically investigate the behavior of this error on $M_{\pm}$ in our numerical results section.

Remark 4.1*.*

Comparing (30) with the stability bound (27) suggests that many of the factors in (30) appear due to bounding the error relative to $\|\widetilde{u}^{\mathcal{N}}\|_{L^{2}}$. Thus the supremum over $g$ is the factor that arises due to the quadrature error.

Remark 4.2*.*

The result (30) also shows that the error between $u^{\mathcal{N}}(s)$ and $\widetilde{u}^{\mathcal{N}}(s)$ is stable independent of $s$ since

[TABLE]

uniformly in $a$, $b$, and $M$. Thus,

[TABLE]

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.1 is that the error in the fully discrete approximation (26) that is due to the $y$-quadrature discretization is computable without solving any PDE’s, assuming that the extremal generalized eigenvalues of $(\boldsymbol{S},\boldsymbol{M})$, coded in the quantities $C_{\mathcal{N}}$ and $K_{\mathcal{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:

[TABLE]

Since $C_{\mathcal{N}}$ is bounded above by the analytical Poincaré constant of the domain $C_{\Omega}$, then choosing this $\mathcal{O}(1)$ quantity for $C_{\mathcal{N}}$ is reasonable. The choices above make the intervals $I_{\pm}$ defined in Proposition 4.1 explicit. The finite element discretization from our numerical experiments in Section 6 results in values $C^{2}_{\mathcal{N}}=0.0506$ and $K^{2}_{\mathcal{N}}=2.36\times 10^{-6}$.

We note that $G_{\pm}$ in (29) 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_{\pm}$ by discretizing the intervals $I_{\pm}$ with $200$ logarithmically spaced points. (I.e., $\log I_{\pm}$ is replaced with 200 equispaced points.) This discretization then allows us to compute $G_{\pm}$, and hence allows us to compute approximations to the bound in (30). In Figure 1, we show the behavior of $G_{\pm}$ 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.

To explore this further, we define the smallest value222To avoid computational effects of oscillating errors due to, e.g., even/odd parity of the quadrature rule, we actually compute the smallest value of $M$ so that $M$, $M+1$, $M+2$, and $(M+3)$-point quadrature rules all achieve the stated accuracy requirement. of $M$ needed to assure a given error level $\delta$:

[TABLE]

For $\delta=10^{-2}$, $10^{-4}$, and $10^{-6}$, we display the values of these $\widetilde{M}$ quantities in Figure 2. We see that for small values of $s$, the requisite number of points $\widetilde{M}$ scales like $1/s$. In particular, for small $s$, more effort (quadrature points) is allocated to $\widetilde{M}_{+}$, but for small $1-s$, comparatively more effort is allocated to $\widetilde{M}_{-}$.

Therefore, the number of quadrature queries $M_{-}+M_{+}$ in the fully discrete scheme (26) 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 [8]. Thus, the number of classical 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 (26) is essentially $M_{-}+M_{+}$ queries of finite element solvers for $\boldsymbol{w}_{\pm}^{\mathcal{N}}$. In practice one can require $M_{-}+M_{+}\sim 100$, cf. Figure 2 and earlier work [16, 40], resulting in a substantial computational cost if the cost of computing $\boldsymbol{w}_{\pm}^{\mathcal{N}}$ is high.

We observe that the formulations (7) for $w_{\pm}$ (and also (21) for the discrete counterparts $\boldsymbol{w}_{\pm}^{\mathcal{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_{\pm}$ are $w^{\mathcal{N}}_{\pm}$ defined in (21). The associated truth solution for $u(s)$ is (26). The purpose of RBM procedures is to diminish the cost of evaluating the truth solution.

5.1. Reduced basis methods

Let $\mathcal{L}$ be a (classical) differential operator, and consider the following PDE parameterized by a Euclidean parameter $y\in D\subset\mathbbm{R}^{p}$:

[TABLE]

where $x$ is the spatial variable and $y$ is a parameter. The operator $\mathcal{L}$ is differential in the $x$ variable. For example, the PDE defining $w(y)$ from (6) can be written as (34) with the operator,

[TABLE]

We assume that (34) is well-posed for each $y\in D$. For a fixed $y$, one usually develops an $x$-discretization with $\mathcal{N}\gg 1$ degrees of freedom yielding a solution $w^{\mathcal{N}}$ with membership in an $\mathcal{N}$-dimensional subspace. For us, this is the discretization defined in section 3.2. We assume that $\mathcal{N}$ is large enough so that

[TABLE]

where $\epsilon$ is a user-prescribed tolerance. Thus, the map $y\mapsto w^{\mathcal{N}}(y)\approx w(y)$ requires algorithms whose complexity is dependent on $\mathcal{N}$.333For linear elliptic operators $\mathcal{L}$, this complexity can in principle scale like $\mathcal{O}\left(\mathcal{N}\log\mathcal{N}\right)$, but frequently is $\mathcal{O}(\mathcal{N}^{2})$, or even $\mathcal{O}(\mathcal{N}^{3})$ depending the details of the employed numerical solver. Such an algorithm that performs the operation $y\mapsto w^{\mathcal{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\mapsto w_{N}(y)\approx w^{\mathcal{N}}(y)$, whose complexity behaves like $\mathcal{O}(N^{3})$ or $\mathcal{O}(N^{2})$, where $N\ll\mathcal{N}$. For $N/\mathcal{N}$ sufficiently small, this can result in an emulator $w_{N}$ whose evaluation is substantially cheaper than the truth approximation $w^{\mathcal{N}}$. The RBM emulator takes the form,

[TABLE]

where $\{y_{k}\}_{k=1}^{N}$ 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,

[TABLE]

is “low rank”. The mathematically precise statement of this is that the Kolmogorov $N$-width of the manifold,

[TABLE]

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^{\mathcal{N}}$ when $N/\mathcal{N}$ is very small. We provide empirical evidence in this paper that this condition is true, and show this rigorously in a follow-up paper [6].

• •

The condition that the truth approximation $w^{\mathcal{N}}(y)$ comes with a practically computable a posteriori error estimate $\Delta(y)$, satisfying,

[TABLE]

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^{\mathcal{N}}(y)\|$. This condition ensures that the parameter values $\{y_{n}\}_{n=1}^{N}$ in (35) 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 $\Delta$ can be derived.

• •

The condition that the operator $\mathcal{L}$ and right-hand side $f$ have affine dependence on the parameter $y$. This means that one has the expressions,

[TABLE]

where we have introduced (i) $y$-independent differential operators $\mathcal{L}_{q}$, (ii) $y$-independent functions $f_{q}(x)$, (iii) $x$-independent functions $\gamma_{q}(y)$, and (iv) $x$-independent functions $\sigma_{q}(y)$. More precisely, one requires the weak (variational) form of $\mathcal{L}$ to have such a decomposition. This condition is needed so that evaluation of the RBM emulator map $y\mapsto u^{N}(y)$ can be accomplished using operations that are independent of the truth approximation discretization parameter $\mathcal{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_{\pm}(y)$ that define the fractional solution $u(s)$. We discuss this in the next section.

5.2. RBM formulation

We describe here the RBM procedure for approximating $\boldsymbol{w}_{-}^{\mathcal{N}}$ via a reduced basis emulator; the procedure for $\boldsymbol{w}_{+}^{\mathcal{N}}$ is nearly identical. We recall the discrete truth approximation formulation that defines $\boldsymbol{w}_{-}^{\mathcal{N}}$:

[TABLE]

where $\boldsymbol{S}$, $\boldsymbol{M}$, and $\boldsymbol{f}$ are defined in (17). 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 [35, 26, 37]. 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 described in (35), the RBM method produces emulator $w_{n,-}$ defined as

[TABLE]

where we have made a particular choice of the basis $\phi_{k}$ appearing in (35).444The parameter values $y_{k}$ and coefficients $c_{n,k}$ should be labeled $y_{k,-}$ and $c_{n,k,-}$, respectively, to differentiate them from the analogous quantities resulting from applying RBM to $w^{\mathcal{N}}_{+}$. However, we omit this notational dependence for more clarity in exposition. Also, since the RBM procedure builds $w_{N}$ sequentially by first building $w_{1},w_{2},\ldots,$, we label the RBM dimension as $n$, satisfying $1\leq n\leq N$, in this section. We must specify the parameter values $\{y_{k}\}_{k=1}^{n}$ and the coefficients $\{c_{n,k}\}_{k=1}^{n}$, which is the focus of the following discussion.

5.3. Computing the $c_{n,k}$

We assume that $y_{1},\ldots,y_{n}$ have been chosen and are known, and now seek to define the coefficients $\{c_{n,k}\}_{k=1}^{n}$, equivalently the vector $\boldsymbol{c}_{n}$, whose computation allow evaluation of $y\mapsto w_{n}(y)$. To proceed we define a new matrix $\boldsymbol{U}\in\mathbbm{R}^{\mathcal{N}\times n}$, having entries

[TABLE]

where the vector $\boldsymbol{w}^{\mathcal{N}}$ and its entries $w^{\mathcal{N}}_{j}$ are expansion coefficients for the truth approximation solution, see (16).

Then for each $y$, the coefficients $c_{n,k}$ of the RBM solution are defined by seeking the vector $\boldsymbol{c}_{n}(y)\in\mathbbm{R}^{n}$ satisfying

[TABLE]

Assuming $\boldsymbol{U}$ has linearly independent columns (which is assured by the choice of $y_{k}$ discussed in the next section), then this uniquely defines $\boldsymbol{c}_{n}(y)$ for each $y$, and prescribes the RBM solution $w_{n}$ via (37). One final point of interest is that our truth variational form (36) exhibits affine dependence on the parameter $y$, making it possible to compute $\boldsymbol{c}_{n}$ very efficiently. We may rearrange computations in (38) so that

[TABLE]

where

[TABLE]

so that the quantities $\boldsymbol{B}$, $\boldsymbol{C}$, and $\boldsymbol{g}$, once computed, are all independent of both the truth discretization dimension $\mathcal{N}$ and the parameter $y$. Thus, for each $y$, the coefficients $\boldsymbol{c}_{n}$ (i.e., the RBM solution $w_{n}$) can be computed with complexity that depends only on $n$ and not on $\mathcal{N}$. Since in practice $n\ll\mathcal{N}$ this can result in computational savings, especially if we wish to query $y\mapsto 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 (37). 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

[TABLE]

Unfortunately, this explicit form requires computing the full solution $w^{\mathcal{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 $\Delta_{n}(y)$.

Lemma 5.1**.**

Define the residual vector $\boldsymbol{r}_{n}(y)\in\mathbbm{R}^{\mathcal{N}}$ as

[TABLE]

and the indicator

[TABLE]

Then

[TABLE]

Proof.

We show the result for the “$-$” quantities; a similar proof works for the “$+$” quantities. The residual $\boldsymbol{r}_{n,-}$ satisfies

[TABLE]

so that

[TABLE]

with

[TABLE]

To summarize, we have the estimate

[TABLE]

which is the desired result by using the definition of $C_{\mathcal{N}}$ in (18). ∎

The residual vectors $\boldsymbol{r}_{n,\pm}$ can be efficiently computed for many values of $y$. We illustrate this $\boldsymbol{r}_{n,-}$. We have:

[TABLE]

where $P_{\mathcal{R}(\boldsymbol{A})}:\mathbbm{R}^{\mathcal{N}}\rightarrow\mathbbm{R}^{\mathcal{N}}$ is the $\mathbbm{R}^{\mathcal{N}}$-orthogonal projector onto the column space of a matrix $\boldsymbol{A}$, and

[TABLE]

The orthogonal decomposition in (43) shows that the Pythagorean theorem can be used to compute the Euclidean vector norm $\|\boldsymbol{r}_{n,-}(y)\|_{2}$ in an efficient way for several values of $y$:

• •

$\|P_{\mathcal{R}(\boldsymbol{R}_{n})^{\perp}}\boldsymbol{f}\|_{2}$ is $y$-independent, so that it can be computed once and stored.

• •

$P_{\mathcal{R}(\boldsymbol{R}_{n})}\left(\boldsymbol{f}-\left(\boldsymbol{S}+e^{-y}\boldsymbol{M}\right)\boldsymbol{U}_{n,-}\boldsymbol{c}_{n,-}(y)\right)$ is a vector in a $2n$-dimensional, $y$-independent vector space. Thus, the norm can be computed with only $n$-dependent complexity. The fact that $\boldsymbol{c}_{n}(y)$ appears with linear behavior in this expression ensures that we can rearrange computations so that, for each $\boldsymbol{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 (42) is not efficiently computable for many values of $y$, the left-hand side is efficiently computable for several values of $y$ since $\|\boldsymbol{r}_{n}\|_{2}$ is efficient to compute, and $\lambda_{\mathrm{min}}\left(\boldsymbol{S},\boldsymbol{M}\right)$ 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 (42), and replace the essentially un-computable maximization (39) with the computable maximization

[TABLE]

The above maximization has an objective function that is efficiently computable, and the inequality (42) ensures that the maximization (44) is a weak greedy algorithm. Weak greedy algorithms in turn ensure that the set of chosen parameters $\{y_{1},\ldots,y_{N}\}$ defines an RBM subspace $V_{N}$ in (35) whose best approximation to the truth solution $\boldsymbol{w}^{\mathcal{N}}_{-}$ is comparable to the Kolmogorov $N$-width [14, 23].

We have completed the basic description of the RBM algorithm: the emulators $w_{N,\pm}(y)$ are defined by computing coefficients $\boldsymbol{c}_{N,\pm}(y)$ as described in the previous section, and the parameter values $y_{k}$ are chosen according to (44) by computing the estimators $\Delta_{n,\pm}(y)$ for $n=1,2,\ldots,$. One usually sequentially computes $y_{k}$ until $\sup_{y\geq 0}\Delta_{n,\pm}(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:

Theorem 5.1**.**

Let $N_{\pm}$ denote the RBM dimensions formed for the emulators $w_{N_{\pm},\pm}(y)$. Then

[TABLE]

Proof.

We have

[TABLE]

Taking $L^{2}$ norms of both sides, and using the triangle inequality with (42) yields the result. ∎

We note that the quantities $\Delta_{N_{\pm},\pm}(y)$ are computed during the RBM construction phase, so that these estimators are available. Therefore, (45) provides a computable error bound that can be used to certify error committed by using the RBM algorithm.

Remark 5.1*.*

Just as with Remark 5.1, all our results above extend to a more general elliptic operator $\mathcal{E}$ satisfying the assumptions outlined in Remark 5.1. 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^{\mathcal{N}}_{\pm}(y)$.555In the previous sections we have describe this process only for $w^{\mathcal{N}}_{-}$, but the process for $w^{\mathcal{N}}_{+}$ results in almost the same procedure with only minor differences stemming from the location of the $e^{-y}$ factor. Subsequently, the efficient RBM emulators $w_{N,\pm}(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 (44) 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 (37) uses solution snapshots $w^{\mathcal{N}}_{\pm}(y)$ as basis functions. This is known to generally lead to ill-conditioning in the formulation (38) even for small $n$ (since in practice the columns of $\boldsymbol{U}_{n,\pm}$ are “nearly” linearly dependent). A better prescription is to build the RBM basis functions by orthogonalizing snapshots.

• •

Some naive implementations of the decomposition (43) via quadratic forms leads to numerical roundoff error that results in stagnation of the error indicators $\Delta_{n,\pm}(y)$ near root-machine precision. (E.g., for double precision, stagnation occurs when $\Delta_{n,\pm}(y)$ takes values around $10^{-8}$. More careful computations allow one to overcome this limitation [19, 20, 17, 21].

We again refer to [35, 26, 37] for a more complete description of important but standard RBM algorithm details.

Finally, we note that one can combine the error estimates in (45) and (30) via the triangle inequality to create a computable bound for $\|u^{N}-\widetilde{u}^{\mathcal{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.

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 (4) on $\Omega=[0,1]^{2}$ with homogeneous Dirichlet boundary conditions. We test a total of 3 algorithms:

• •

“SQ” — The sinc quadrature approach from [16].

• •

“GQ” — Algorithm 1 detailed in section 3, utilizing a modified version of the integral formulation in [16] 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 classical PDE solves needed to compute an approximation to $u(s)$. In particular, we make the choices given in [16, 8], which depend on the spatial mesh and fractional order:

[TABLE]

The finite element discretization is accomplished with linear quadrilateral finite elements on a Cartesian tessellation of $\Omega$. 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$.

6.1. Manufactured Solutions on $[0,1]^{2}$

Consider the physical domain $\Omega=[0,1]^{2}$. In this case, an explicit family of eigenfunctions for the Laplacian with homogeneous Dirichlet boundary conditions is available:

[TABLE]

which satisfy

[TABLE]

With this in hand, and using the inverse of the relation (3), 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,

[TABLE]

• •

“Mixed modes” — The function $u$ and data $f$ are, in this case,

[TABLE]

• •

“Square bump” — The data $f$ is the indicator function

[TABLE]

While an analytical solution is available as an infinite sum of eigenfunctions, we instead numerically compute a solution to the above problem on a highly refined mesh and consider this the “exact” solution.

6.2. Spatial convergence

Our first test verifies that we recover spatial convergence in terms of the finite element mesh size. Since are primarily interested in accuracy and not efficiency, this section compares the SQ and GQ methods.

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^{\mathcal{N}}_{\pm}$ 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_{\pm}$ as given by (46). 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.

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 $\Omega$, $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 (44) by discretizing the $y\in[0,\infty)$ domain in a uniform way under a logarithmic map. Precisely: we set $z=e^{-y}$ for $y\in[0,\infty)$ and proceed to discretize $z\in[0,1]$. We take 128 equispaced points in the $z$ variable and map back to $y$-space with $z\mapsto-\log z=y$. We subsequently perform a discrete maximization over this set instead of the continuous optimization (44).

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\gg 1$ quadrature point is essentially the same as that required to evaluate at a single quadrature point.

6.5. RBM accuracy

We now investigate the accuracy delivered by the RBM algorithm in the construction of reduced order models for $w^{\mathcal{N}}_{\pm}$. Our rigorous error certificate for $u(s)$ using the reduced order model is (45), but we here consider a finer estimate using the proof of Theorem 5.1. We define the error estimator,

[TABLE]

where $\Delta_{N,\pm}$ are the computable error indicators defined in (41), and we again choose $M_{\pm}$ as in (46). One can see from the proof of Theorem 5.1 that this quantity bounds the error committed by the RBM procedure. We show the values of this quantity in Figure 6 as a function of $N$, and observe that it decays exponentially.

Finally, we remark that $\Delta_{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\mapsto 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\mapsto 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.

7. Conclusion

We propose a novel model reduction strategy for computing solutions to fractional Laplace PDE’s, in particular (4). Our algorithm builds on the ideas introduced in [16], 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.