A simple solver for the fractional Laplacian in multiple dimensions

A simple solver for the fractional Laplacian in multiple dimensions

Victor Minden Institute for Computational and Mathematical Engineering, Stanford University, Stanford, CA 94305. Current address: Center for Computational Biology, Flatiron Institute, Simons Foundation, New York, NY 10017 (vminden@flatironinstitute.org). Funding: U.S. Department of Energy Advanced Scientific Computing Research program (grant number DE-FC02-13ER26134/DE-SC0009409).    Lexing Ying Department of Mathematics and Institute for Computational and Mathematical Engineering, Stanford University, Stanford, CA 94305 (lexing@stanford.edu).Funding: National Science Foundation (grant number DMS-1521830) and U.S. Department of Energy Advanced Scientific Computing Research program (grant number DE-FC02-13ER26134/DE-SC0009409).

We present a simple solver of the fractional Laplacian based on the hypersingular integral representation. Through singularity subtraction, we obtain a regularized integrand that is amenable to the trapezoidal rule with equispaced nodes. The resulting quadrature scheme gives a discrete operator on a regular grid that is translation-invariant and thus can be applied quickly with the fast Fourier transform. For discretizations of problems related to space-fractional diffusion on bounded domains, we observe that the underlying linear system can be efficiently solved via preconditioned Krylov methods with a preconditioner based on the finite-difference (non-fractional) Laplacian. We show numerical results illustrating the error of our simple scheme as well the efficiency of our preconditioning approach, both for the elliptic (steady-state) fractional diffusion problem and the time-dependent problem.

1 Introduction

Fractional powers of the Laplacian operator arise naturally in the study of anomalous diffusion, where the fractional operator plays an analogous role to that of the standard Laplacian for ordinary diffusion (see, e.g., the review articles by Metzler and Klafter [34, 35] and Vázquez [44]). By replacing Brownian motion of particles with Lévy flights [28], whose increments are drawn from the -stable Lévy distribution for , we obtain a fractional diffusion equation (or fractional kinetic equation) in terms of the fractional Laplacian operator of order [41], defined for sufficiently nice functions via the Cauchy principal value integral


with known normalizing constant [25].

For a bounded domain with complement , we consider fractional diffusion with homogeneous extended Dirichlet conditions given in terms of (1) by


Also of interest is the related elliptic problem


Somewhat unintuitively, the nonlocality of (1) implies that the solutions of (2) and (3) depend on data prescribed everywhere outside [11, 39, 9], though other definitions of the fractional Laplacian on a bounded domain are also in common use [44]. Further, a more general formulation of fractional diffusion involves augmenting (2) by incorporating fractional time derivatives of Caputo or Riemann-Liouville type. We focus in this work on the case of space-fractional diffusion and do not discuss the discretization of time-fractional differential operators, though the latter is of independent interest [31, 26, 49, 50].

1.1 Contribution

The contribution of this paper is a simple discretization scheme for (2) and (3) on Cartesian grids, and an efficient algorithm for solving the resulting linear systems. The discretization generalizes easily to domains that can be represented as occluded Cartesian grids, though we focus on regular domains in this work.

Our approach is based on using a Taylor expansion around each point to replace the singular integrand in Eq. 1 with a sufficiently smooth function of on all of via singularrity subtraction. The resulting integral can be easily discretized using the trapezoidal rule on a regular grid of points, leading to a translation-invariant linear operator that can be applied at a cost of using the fast Fourier transform (FFT). The resulting discrete linear system approximating (3) can then be efficiently solved using standard Krylov methods. As , the resulting linear systems can exhibit the ill-conditioning characteristic of discretizations of the Laplacian operator on a regular grid. To circumvent this, we develop an efficient preconditioning strategy based on the fact that our discrete fractional Laplacian operator may be written as the sum of a standard finite-difference Laplacian and another matrix with mostly small entries.

When the solution to is sufficiently smooth, standard results on convergence of the trapezoidal rule and finite-difference operators imply that the error of our approach for computing the fractional Laplacian at a point goes to zero as , where is the linear spacing between grid points, which we show in Section 2. In general, however, the solution to the fractional Laplace problem on bounded domains is only times continuously differentiable [40], leading to a natural deterioration of the rate of convergence of our simple approach to an observed rate of .

1.2 Related work

A discretization scheme similar to that presented here appears in Pozrikidis [38], though without discussion of accuracy or the importance of windowing for singularity subtraction. Huang and Oberman [22, 23] derive a scheme for the one-dimensional case based on singularity subtraction and finite-difference approximation, but do not tackle the multidimensional case (see also Tian and Du [42], Gao et al. [16], and Duo, Van Wyk, and Zhang[10]). Chen et al. [7] consider a multidimensional discretization and fast preconditioners based on multigrid, but their scheme uses the so-called “coordinate fractional Laplacian” that takes a tensor product of one-dimensional operators and is not equivalent to (1) (see also related finite-difference approaches with different operators [51, 45, 33, 32]).

Other work on efficient solution of fractional Laplacian systems using fast preconditioned iterative methods includes Wang and Du [46], Pang and Sun [37], Fu and Wang [15], and Fu, Ng, and Wang [14], though only in one spatial dimension.

Notable schemes for discretizing the fractional Laplacian based on different ideas include work based on the Caffarelli-Silvestre extension [6, 36, 20], finite-element-based approaches [2, 5, 4, 43], and work on spectral approaches [48, 29, 3]. General references for fractional Laplacians on bounded domains include, e.g., Ros-Oton [39], D’Elia and Gunzburger [9], Felsinger [11], and Lischke et al. [27].

2 Spatial discretization of the fractional Laplacian

To begin, we outline our scheme for discretization of (1) in the one-dimensional case where the function vanishes outside of some interval. Following that, we give more details in our discussion of the multidimensional case.

2.1 Singularity subtraction in one dimension

Concretely, consider the task of approximating the integral


where for . For , this integral is hypersingular due to the high-order pole at , which generally leads to large inaccuracies when simple quadrature schemes are applied directly to (4). Therefore, we proceed by regularizing the integrand to remove the singularity and obtain an integral for which simple quadratures are accurate.

Assuming that the function is sufficiently smooth, we may write a Taylor series expansion about the point to obtain


where the smooth remainder as . For brevity, we have grouped terms that are odd about into , as they will not play an explicit role in what follows.

Our regularization strategy is singularity subtraction based on adding and subtracting a calculable integral that matches terms in the Taylor series. Suppose is a sufficiently smooth windowing function with compact support such that and . Then we may write


where we define (I) to be the first integral and (II) to be the second. By construction, (I) is no longer hypersingular, as we see from (5) that the integrand can be equivalently written

By our smoothness assumptions on and , this integrand is and continuously differentiable as with a second derivative that is integrable. This implies that the standard trapezoidal rule would exhibit second-order convergence when applied to (I); see Cruz-Uribe and Neugebauer [8]. Of course, this requires knowledge of and in general, which we do not assume. In the context of discretization of the integral using a uniform grid, however, the situation simplifies.

2.2 The first integral in one dimension

Consider discretizing (I) using the trapezoidal rule on a one-dimensional lattice and take to be one of the lattice points. Without loss of generality, we may shift the domain such that . This discretization yields the second-order accurate approximation


where and are constants independent of and and the last sum in the second line is identically zero due to oddness considerations. We note that the sum remaining on the final line is over a finite range, as is compactly supported. Since is assumed to be smooth, we replace with the finite-difference approximation

which gives our final approximation for (I),


2.3 The second integral in one dimension and final quadrature

Having established a method for approximating the integral (I) in (6), we turn to (II). Again using oddness considerations, we see that the contribution from vanishes such that


where the constant given by

is well-defined (since is compactly supported) and we again take for convenience. We once again replace the second derivative with its finite-difference approximation to obtain Combining this with our quadrature for (I) gives our approximation for ,

which applies equally well not only to but in general to for any grid point , i.e.,


This is our final quadrature for the fractional Laplacian in one dimension.

2.4 Singularity subtraction in higher dimensions

We turn now to the multidimensional integral, i.e., (1) with or . Once again we will assume that the function is compactly supported and sufficiently smooth, as we will make explicit. Our basic strategy is the same as in one dimension.

Lemma 1.

Suppose that and let be a windowing function symmetric about such that as . Let the third-order Taylor approximation of about the point be given in multi-index notation by


where the remainder is given in explicit form as

Then, defining the function


we have that and as for , , and .


It is clear that

By inspection, the order of differentiability of is limited by that of and of . Given the explicit form of , it is at least in as a function of , whereas by assumption. Further, for , since the first summand is and the second summand is at least . Explicit term-by-term differentiation of with the product rule concludes the proof. ∎

By subtracting off the windowed multivariate Taylor series we obtain an integral that is no longer hypersingular. In particular, we write


where we define (Id) to be the first integral and (IId) to be the second.

2.5 The first integral in higher dimensions

To numerically approximate (Id) we use a quadrature rule on a uniform lattice . We assume the lattice is constructed such that the point coincides with with some lattice point , which we take to be without loss of generality.

Replacing the integral with a weighted sum over the lattice, we obtain


which we note does not include a term for . This corresponds to the standard trapezoidal rule for and the punctured trapezoidal rule for , though more involved quadrature corrections may be used (see, e.g., Marin, Runborg and Tornberg [30]). Assuming is symmetric about the origin, we see that for many values of the multi-index the corresponding summand vanishes due to oddness considerations. Taking these symmetries into account, we let denote the first coordinate of and observe that

which we plug back into our quadrature scheme to obtain


with correspondingly defined constants

Theorem 1.

Suppose the same setup as Lemma 1 with , , and such that and . Assume further and are compactly supported with for all . Then the above approximation for (Id) is second-order accurate. That is,

with and as in (12).


The described approximation is numerically equivalent to the (punctured) trapezoidal rule, so this amouts to bounding the error of the trapezoidal rule applied in dimensions with integrand , where is as in Lemma 1 with . Letting be such that both and for , we proceed by breaking the integral into three contributions: one for the subdomain “near” the singularity, one for the “mid-range” subdomain , and one for the “far” subdomain . We write

each piece of which we analyze separately.

Near the singularity, we see due to symmetry considerations that

where under our assumptions the integrands are both and is bounded. Since , we see , which implies that the corresponding (punctured) trapezoidal rule approximation to the integral is , since we gain a factor of due to the quadrature weights. Therefore, the contribution to the error from the integral over the near subdomain is , since .

In the mid-range subdomain, we explicitly use the composite nature of the trapezoidal rule to write

and then consider the error of the trapezoidal rule in approximating the integral over each separately, where the square/cubic subdomains in the trapezoidal rule are pairwise disjoint and are such that . Since we are away from the origin, on each subdomain the integrand is in which means the standard error bound for the trapezoidal rule on gives an error contribution of no more than for some constant independent of . However, the term does depend on . Since from Lemma 1, the product rule gives . With this we can bound the total error on as

where we have bounded

(up to some geometry-dependent factors that are independent of ) due to concavity of the summand. Therefore, the error contribution from the mid-range subdomain is .

Finally, for the far subdomain, we observe that the integrand is in and its smoothness is independent of in this region, so the standard composite trapezoidal error bound of applies. Therefore, the overall error is . ∎

2.6 The second integral in higher dimensions and final quadrature

We now consider the second integral (IId) in (11). Assuming without loss of generality that and using symmetry and oddness considerations as before, we see that


Defining the constant


and combining this with our quadrature for (Id) gives

or, more generally,

Of course, as written this approximation requires second derivative information in the form of . For smooth , however, we may replace this with a finite-difference stencil involving the neighbors of in the lattice,

just as in the one-dimensional case.

2.7 Summary of quadrature for fractional Laplacian

We briefly summarize our complete approach for discretizing the fractional Laplacian applied to a function . First, we regularize the integrand of (1) by adding to the numerator a windowed Taylor series approximation of about with window function to obtain (Id) in (11). This gives an integral that is nice enough to admit discretization with the trapezoidal rule or related schemes. Then, by exploiting symmetries of the problem, we rewrite the discretization in terms of the constants and in (12), which do not depend on . Finally, we derive an expression for the correction term (IId) in terms of another constant given in (13), which when combined with (Id) and a finite-difference stencil approximation gives a nice expression for as a linear function of evaluated on a regular grid.

A few details of the procedure remain to be discussed. First, there are a number of possibilities for the windowing function . In this paper, we use the piecewise-polynomial window


Of course, this is by no means the only sufficiently smooth choice. Further, we note that the requirement that be compactly supported can be relaxed so long as decays sufficiently quickly as such that the necessary integrals and sums may be computed.

On that note, we also must still compute the constants , , and . For our choice of polynomial window, the integral defining can be computed explicitly; for other choices the integral may be numerically computed to high precision offline using, e.g., adaptive quadrature in MATLAB. For compactly supported , the sum defining has a finite number of nonzero terms and is easily computable. Finally, the infinite lattice sum is given in terms of the Riemann zeta function for and may otherwise be well-approximated using far-field compression techniques related to the fast multipole method (FMM) [18, 47]. We use Chebyshev polynomials for far-field compression in the vein of Fong and Darve [13], though we do not require the full FMM machinery as we are interested only in the lattice sum and not a full approximate operator.

3 Solving the fractional differential equations on a bounded domain

Having developed our trapezoidal rule scheme for evaluating (1) given , we turn now to the fractional differential equations (2) and (3) concerning fractional diffusion on a bounded domain with homogeneous extended Dirichlet conditions. We focus on the case for simplicity.

3.1 The elliptic case: steady-state fractional diffusion

To solve the elliptic problem (3), we discretize using a regular grid of points with linear spacing , where and For notational convenience, we define the index set . Then, replacing the fractional Laplacian with our quadrature-based approximation gives


which is a linear system to be solved for the variables . We remark that the “boundary conditions” affect the system in two ways. First, the center sum has been reduced from an infinite number of terms (in general) to a more manageable finite sum. Second, evaluating the finite-difference stencil for near the boundary of the domain will require the prescribed value of on the boundary, as in the standard (non-fractional) case.

We write (15) in matrix form as


where now and are vectors with corresponding entries and and contains the coefficients implied by (15).

Forward operator and application with FFT

By construction, the approximate fractional Laplacian operator involved in (15) is translation-invariant, which means that the matrix is block Toeplitz with Toeplitz blocks (BTTB) under any natural ordering of the unknowns. As is well known, this in turn implies that may be applied efficiently using the FFT at a cost of FLOPs per application and stored with storage cost .

Further, investigation of the constants , and reveal that is symmetric positive definite. When coupled with the previous observatiton, this leads naturally to the use of the conjugate gradient method (CG) [19] or related iterative methods for solving (16). However, while the FFT ensures low complexity per iteration, the number of iterations required to achieve a specified iteration can be large unless an effective preconditioner is used. This is of particular concern as , whereupon we recover the standard (ill-conditioned) Laplacian.

Preconditioning: Laplacian pattern and fast Poisson solver

To construct an efficient preconditioner for (16), we observe that may be decomposed as the sum of two matrices , where

and we note that . The sparse matrix is (up to a proportionality constant) the typical finite-difference approximation of the negative Laplacian, whereas the matrix has entries that quickly decay away from , particularly for larger . This motivates using itself as a preconditioner when using CG to solve (16). Because is effectively a finite-difference discretization of Poisson’s equation on a regular grid with homogeneous Dirichlet boundary conditions, application of may be accomplished with the FFT at a cost of using typical fast Poisson solver techniques [24, Chapter 12]. For non-rectangular domains, the FFT-based approach is no longer feasible, but the same preconditioner could be used with, e.g., nested dissection [17] or related methods.

We remark that other choices of preconditioner are possible. For example, rather than using as our preconditioner we could instead use , where if and zero otherwise. Preliminary experiments with this approach (not shown) did not seem to show measurable benefit.

3.2 The time-dependent case: time-dependent fractional diffusion

We turn now to the full time-dependent problem (2). For spatial discretization we use the approximate fractional Laplacian just as in Section 3.1, which we combine with a Crank-Nicolson scheme for the discretization of temporal derivatives. This leads to the implicit time-stepping method


to be solved for , where is as in Section 3.1 and now and have entries and for .

Just as in Section 3.1, we exploit BTTB structure to apply such that (17) may be solved efficiently with CG at each time step. Compared to the steady-state problem, the system matrix here is much better conditioned due to the addition of the identity. However, we still find that the number of iterations is reduced substantially via preconditioning, where we use the matrix as preconditioner.

4 Numerical results

To demonstrate and profile our approach to discretizing and solving fractional diffusion problems on bounded Cartesian domains, we implemented a number of examples. All computations were performed in MATLAB R2017a on a 64-bit Ubuntu laptop with a dual-core Intel Core i7-7500U processor at 2.70 GHz and 16GB of RAM.

4.1 Elliptic example: interval in one dimension

We begin with a one-dimensional elliptic example on the interval , discretizing and solving (3) with right-hand side


given in terms of the Gauss hypergeometric function [1]. In this case, the analytic solution on is known and is given up to a known constant of proportionality by [21].

We discretize the interval with regularly-spaced points as in (15) with , recording the total time to construct the discrete operator . We choose in (14) as a function of the number of discretization points , such that is supported inside a ball with a radius of 20 discretization points. Using the known solution for right-hand side (RHS) (18), we record the time to apply our operator using the FFT, , and measure the apply error of our discretization as where is the analytic solution sampled on the discrete grid points and is the discretized RHS. To demonstrate the solution error of our discretization scheme we take the same RHS as before and use CG to solve the resulting linear system (16). This gives a discrete solution that we can compare to by computing the relative solution error These metrics are all shown in Table 1 for three different choices of , with the error metrics plotted in Fig. 1. Note that and do not depend on , so we report these quantities only for .

We show in Table 2 the runtime and iterations required by CG to solve the linear system (16) for two different choices of tolerance in the relative -norm residual. We give results both for the preconditioned system (where the preconditioner is a finite-difference Laplacian as described in Section 3.1) and the unpreconditioned system. Because this is a one-dimensional problem, use of a fast Poisson solver to apply is not stricly necessary for efficiency. Instead, we use a sparse Cholesky factorization, with negligible overhead. The corresponding timing results are plotted in Fig. 2 (left), where we see that our simple preconditioning scheme is effective for reducing the time to solution, especially for larger .

Table 1: Runtimes for the construction of the operator and for application via FFT for the one-dimensional elliptic example, as well as relative apply errors for .