Fast algorithm for computing nonlocal operators with finite interaction distance

TL;DR

This paper introduces a novel fast algorithm for nonlocal operators with finite interaction distance, significantly reducing computational costs and enabling more efficient simulations of nonlocal phenomena.

Contribution

A new class of fast algorithms is developed for nonlocal operators, overcoming previous limitations and matching the efficiency of local numerical methods in some cases.

Findings

Reduces computational complexity of nonlocal operators

Applicable to models like peridynamics and fractional diffusion

Enables efficient large-scale simulations

Abstract

Developments of nonlocal operators for modeling processes that traditionally have been described by local differential operators have been increasingly active during the last few years. One example is peridynamics for brittle materials and another is nonstandard diffusion including the use of fractional derivatives. A major obstacle for application of these methods is the high computational cost from the numerical implementation of the nonlocal operators. It is natural to consider fast methods of fast multipole or hierarchical matrix type to overcome this challenge. Unfortunately the relevant kernels do not satisfy the standard necessary conditions. In this work a new class of fast algorithms is developed and analyzed, which is some cases reduces the computational complexity of applying nonlocal operators to essentially the same order of magnitude as the complexity of standard local numerical methods.

Tables

Table 1. Table 1: The storage and time of matrix-vector multiplication for the 1d nonlocal diffusion in homogenous medium. The domain Ω = ( 0 , 1 ) Ω 0 1 \Omega=(0,1) is discretized using N = 2048 𝑁 2048 N=2048 points. The mesh size h = 1 / 2048 ℎ 1 2048 h=1/2048 . The relative precision of HIF approximation is set to be ϵ = 10 − 8 italic-ϵ superscript 10 8 \epsilon=10^{-8} .

Table 2. Table 2: The storage and time of matrix-vector multiplication for the 1d nonlocal diffusion in homogenous medium. δ = 1 / 4 𝛿 1 4 \delta=1/4 . The relative precision of HIF approximation is set to be ϵ = 10 − 8 italic-ϵ superscript 10 8 \epsilon=10^{-8} .

Table 3. Table 3: The storage and time of matrix-vector multiplication for the 1d nonlocal diffusion in heterogeneous medium. The domain Ω = ( 0 , 1 ) Ω 0 1 \Omega=(0,1) is discretized using N = 2048 𝑁 2048 N=2048 points. The mesh size h = 1 / 2048 ℎ 1 2048 h=1/2048 . The relative precision of HIF approximation is set to be ϵ = 10 − 8 italic-ϵ superscript 10 8 \epsilon=10^{-8} .

Table 4. Table 4: The storage and time of matrix-vector multiplication for 1d nonlocal diffusion in heterogeneous medium. Kernel γ ​ ( s ) 𝛾 𝑠 \gamma(s) is C 3 superscript 𝐶 3 C^{3} away from origin. δ 0 = 1 / 4 subscript 𝛿 0 1 4 \delta_{0}=1/4 . The relative precision of HIF approximation is set to be ϵ = 10 − 8 italic-ϵ superscript 10 8 \epsilon=10^{-8} .

Table 5. Table 5: The storage and time of matrix-vector multiplication for 2d nonlocal diffusion in homogeneous medium. The domain Ω = ( 0 , 1 ) 2 Ω superscript 0 1 2 \Omega=(0,1)^{2} is discretized using N = 64 × 64 𝑁 64 64 N=64\times 64 points. The mesh size h = 1 / 64 ℎ 1 64 h=1/64 . The relative precision of HIF approximation is set to be ϵ = 10 − 3 italic-ϵ superscript 10 3 \epsilon=10^{-3} .

Table 6. Table 6: The storage and time of matrix-vector multiplication for 2d nonlocal diffusion in homogeneous medium. δ = 1 / 2 𝛿 1 2 \delta=1/2 . The relative precision of HIF approximation is set to be ϵ = 10 − 3 italic-ϵ superscript 10 3 \epsilon=10^{-3} .

Table 7. Table 7: Computation time for evaluating the nonlocal operator in 2 ​ d 2 𝑑 2d . k = 0 𝑘 0 k=0 and p 0 = 1 superscript 𝑝 0 1 p^{0}=1 . T multi subscript 𝑇 multi T_{\text{multi}} under the new algorithm is the time for Step 2 and Step 3 combined. Original -1 uses the MATLAB built-in sparse matrix-vector multiplication. Original -2 uses a for-loop to go through 𝒙 i ​ ( i = 1 , 2 , ⋯ , N ) subscript 𝒙 𝑖 𝑖 1 2 ⋯ 𝑁 \bm{x}_{i}(i=1,2,\cdots,N) together with the MATLAB built-in vector-vector multiplication applied in each loop.

Table 8. Table 8: Computation time for evaluating the nonlocal operator in 1d. k = 3 𝑘 3 k=3 and p 3 superscript 𝑝 3 p^{3} is an even polynomial of degree 6 6 6 . δ ​ ( x ) 𝛿 𝑥 \delta(x) is given by Fig 1 with δ 0 = 1 / 4 subscript 𝛿 0 1 4 \delta_{0}=1/4 . T multi subscript 𝑇 multi T_{\text{multi}} under the new algorithm is the time for Step 2 and Step 3 combined. Original -1 uses the MATLAB built-in sparse matrix-vector multiplication. Original -2 uses a for-loop to go through 𝒙 i ​ ( i = 1 , 2 , ⋯ , N ) subscript 𝒙 𝑖 𝑖 1 2 ⋯ 𝑁 \bm{x}_{i}(i=1,2,\cdots,N) together with the MATLAB built-in vector-vector multiplication applied in each loop.

Equations

$$-{\int }_{{\mathcal{H}}_{\mathbf{x}}}\mathbf{f}(\mathbf{u}(\mathbf{y})-\mathbf{u}(\mathbf{x}),\mathbf{y},\mathbf{x})d\mathbf{y}=\mathbf{b}(\mathbf{x}),$$

$$\mathbf{f}(\mathbf{u}(\mathbf{y})-\mathbf{u}(\mathbf{x}),\mathbf{y},\mathbf{x})=\omega (\mathbf{x},\mathbf{y})[\mathbf{e}\cdot (\mathbf{u}(\mathbf{y})-\mathbf{u}(\mathbf{x}))]\mathbf{e},$$

$${\int }_{{\mathcal{H}}_{\mathbf{x}}}|\mathbf{y}-\mathbf{x}{|}^2\omega (\mathbf{x},\mathbf{y})d\mathbf{y}<\infty .$$

$$\mathcal{L}u(\mathbf{x})={\int }_{{\mathcal{H}}_{\mathbf{x}}}\omega (\mathbf{x},\mathbf{y})(u(\mathbf{y})-u(\mathbf{x}))d\mathbf{y}.$$

$$\omega (\mathbf{x},\mathbf{y})=\frac{C(\mathbf{x},\mathbf{y})}{\delta (\mathbf{x},\mathbf{y}{)}^{d+2}}\gamma (\frac{|\mathbf{y}-\mathbf{x}|}{\delta (\mathbf{x},\mathbf{y})}),$$

$$C_0\le C(\mathbf{x},\mathbf{y})\le C_1,\text{ and }0<\delta (\mathbf{x},\mathbf{y})\le {\delta }_1.$$

$$\gamma (|s|)=\frac{c_1}{|s|}\chi (|s|<1),\text{ or }\gamma (|s|)=\frac{c_2}{|s|}(1-|s|)\chi (|s|<1).$$

$$\mathcal{L}u(\mathbf{x})\to C(\mathbf{x})\Delta u(\mathbf{x}).$$

$$\mathcal{L}u(\mathbf{x})\to \nabla (\sigma (\mathbf{x})\nabla u(\mathbf{x})),\text{ where }\sigma (\mathbf{x})=C(\mathbf{x},\mathbf{x}).$$

$$\gamma (|s|)=\frac{1}{|s|}\chi (|s|<1).$$

$$\gamma (|s|)=\kappa (|s|)+p^{2K}(s)\chi (|s|<1),$$

$$p^{2K}(s)=\{\begin{array}{rlr}& 1& K=0\\ & \frac{3}{2}-\frac{1}{2}{s}^2& K=1\\ & \frac{15}{8}-\frac{5}{4}{s}^2+\frac{3}{8}{s}^4& K=2\\ & \frac{17}{8}-2s^2+\frac{9}{8}{s}^4-\frac{1}{4}{s}^6& K=3\end{array}$$

$$\{\begin{array}{rlr}& -{\mathcal{L}}^{P}u(\mathbf{x}):=-{\int }_{{\mathcal{H}}_{\mathbf{x}}}P(\mathbf{x},\mathbf{y})(u(\mathbf{y})-u(\mathbf{x}))d\mathbf{y}=f(\mathbf{x})& x\in \Omega \\ & u(\mathbf{x})=0& x\in {\Omega }^c\end{array}$$

$$P(\mathbf{x},\mathbf{y})=\frac{C(\mathbf{x},\mathbf{y})}{\delta (\mathbf{x},\mathbf{y}{)}^{d+2}}{p}^{2K}(\frac{|\mathbf{y}-\mathbf{x}|}{\delta (\mathbf{x},\mathbf{y})}).$$

$${\mathcal{L}}^{P}u(\mathbf{x})=\frac{C(\mathbf{x})}{\delta (\mathbf{x}{)}^{d+2}}\left({\int }_{{\mathcal{H}}_{\mathbf{x}}\cap \Omega }u(\mathbf{y})d\mathbf{y}-|{\mathcal{H}}_{\mathbf{x}}|u(\mathbf{x})\right).$$

$${\mathcal{N}}_d({\mathbf{x}}_i)=\#\bigcup _{l=0}^L\{Q_l^j\in L_l:(Q_l^j{)}^{\mathrm{o}}\cap \partial {\mathcal{H}}_{{\mathbf{x}}_i}\ne \varnothing \}=\sum _{l=0}^L{\mathcal{N}}_d^l({\mathbf{x}}_i),$$

$${\mathcal{N}}_d^l({\mathbf{x}}_i)=\{\begin{array}{rlr}& O(1)& \text{for }d=1;\\ & O{\left(\frac{N}{(2^d{)}^l}\right)}^{1-1/d}& \text{for }d\ge 2.\end{array}$$

$$\{\begin{array}{r}O(N\log N)\text{for }d=1;\\ C_d\cdot N^{2-1/d}\text{for }d\ge 2,\end{array}$$

$${\int }_{{\mathcal{H}}_{{\mathbf{x}}_i}\cap \Omega }u(\mathbf{y})d\mathbf{y}=\sum _{l=0}^L\sum _{j\in {\mathcal{I}}_l({\mathbf{x}}_i)}{\int }_{Q_l^j}u(\mathbf{y})d\mathbf{y}.$$

$$\{\begin{array}{r}O(KN\log N)\text{for }d=1;\\ O(K{N}^{2-1/d})\text{for }d\ge 2,\end{array}$$

$${\mathcal{L}}_{\delta }^{P}u(\mathbf{x})=\frac{1}{{\delta }^4}{\int }_{B_{\delta }(\mathbf{x})}(u(\mathbf{y})-u(\mathbf{x}))d\mathbf{y}=\frac{1}{{\delta }^4}{\int }_{B_{\delta }(\mathbf{x})\cap \Omega }u(\mathbf{y})d\mathbf{y}-\frac{|B_1|}{{\delta }^2}u(\mathbf{x}).$$

Full text

{AMS}

65R20;82C21;65F05;45K05

1 Introduction

Nonlocal operators as alternatives to the local differential operators have been widely applied to the modeling of a number of physical and biological processes over the last few years. One example is given by the peridynamics models for brittle materials [25, 3, 10, 14, 15, 26], and another is nonstandard diffusion including the use of fractional derivatives [1, 4, 6, 9, 21, 20]. Moreover, even in the case of solving differential equations, nonlocal integral relaxations to local differentiations are sometimes introduced as a numerical technique, see [11, 19, 22, 23].

While the nonlocal models can provide more accurate descriptions of a physical system, the nonlocality also increases the computational cost compared to classical local models based on partial differential equations (PDEs). Thus it is imperative to develop fast algorithms for the computation of nonlocal models. In this work, we focus on nonlocal models with finite nonlocal interaction distances (see [28]). Each point in the domain interacts with points within certain “horizon” of it. That is to say the interaction kernel is compactly supported in certain finite region and may also be singular at zero distance. When the interaction region is spherical, we call the radius of the region the horizon radius or just horizon in short. Such kernels make most fast algorithms that utilizing the low rank property of the far field interactions fail to work, because the kernel is typically not regular around the boundary of its support. In the work of [33, 31], FFT based fast algorithms are developed for solving peridynamic models, in which Toeplitz structures of matrices are greatly exploited. However, such methods fail to work for the more generals cases of inhomogeneuous media and variable horizons resulting in non-translation invariant kernels, which is the concern of this work.

The peridynamics (PD) theory of solid mechanics takes into account of large scale deformations and reformulates the internal force of a deformed body by introducing the nonlocal way of interactions of materials points. Let $\mathbf{u}$ denote the displacement field, the force balance equation in the bond-based peridynamic theory is given by

[TABLE]

where the integration is limited to the finite region $\mathcal{H}_{\bm{x}}$. $\mathbf{b}=\mathbf{b}(\bm{x})$ denotes the body force and $\mathbf{f}$ is the force density function, which contains all the constitutive properties. For a prototype micro-elastic material, the force vector $\mathbf{f}$ points to the direction of the deformed bond with its magnitude a function of the bond relative elongation(stretch). In general, $\mathbf{f}$ depends nonlinearly on the bond relative elongation, which gives the model the flexibility to describe changes of material properties. In this work, we will only consider the linear model obtained from the assumption of small deformation, namely,

[TABLE]

where $\mathbf{e}$ is the unit vector in the bond direction $\bm{y}-\bm{x}$, and $\omega(\bm{x},\bm{y})$ is the micromodulus function satisfying

[TABLE]

For each $\bm{x}$, $\omega(\bm{x},\cdot)$ is a function with support on $\mathcal{H}_{\bm{x}}$. In practice, this is achieved by choosing a fractional type kernel multiplied with the characteristic function or the conical function on $\mathcal{H}_{\bm{x}}$ (see [14]), where both choices create non-smooth transition to zero at the boundary $\partial\mathcal{H}_{\bm{x}}$.

As we can see, the discretization of the nonlocal operator in (1.1) results in a matrix $A$ with high density, for which fast algorithms must be considered in order to lower the cost of matrix multiplication and inversion. There are several classical algorithms for the treatment of non-sparse matrices, including the fast multipole methods (FMMs) [12, 13, 32] and the hierarchical matrix techniques [16, 17]. More recent developments can be found in [18]. We will first consider the existing fast solvers to see whether they are effective in the computation of nonlocal models like (1.1). In Section 2 we will show that the performance of the fast solvers depends largely on the regularity of the kernel across the boundary $\partial\mathcal{H}_{\bm{x}}$. This sends the clear message of the preference of the far field behavior of the kernel function in terms of computation efficiency. Popular choices of kernels in the peridynamics simulation will cause serious problems in the interest of fast computation. We then propose in Section 3 to split the two types of singularities – one at origin and the other at the boundary $\partial\mathcal{H}_{\bm{x}}$– of the kernel function, and treat them separately. Once the kernel is decomposed into two, one that is smooth away from origin and the other that is smooth away from the truncation, we then use different fast algorithms that exploit the different types of smoothness of the two kernels. The classical fast algorithms such as FMMs and hierarchical matrix techniques only resolve the first type of singularity, namely that the singularity set is one point (the origin of the kernel in most cases). Notice that the second type of singularity is on the boundary set $\partial\mathcal{H}_{\bm{x}}$, and it is a set of codimension $1$. In 2d the sets of singularity are boundary curves, and in 3d they are boundary surfaces. Our algorithm that deals with second type of singularity is basically a new FMM type method for kernels that exhibit singularities on codimension 1 sets. There is a reduction of computational complexity in all dimension with optimal rate $O(N\log N)$ for $N$ unknowns in 1d. Other potential applications of this work include computing with the retarded potentials raised in time domain boundary integral equations [24], where the potentials are discontinuous functions defined in space-time.

2 The role of far field behaviors of kernels

We will use the nonlocal diffusion operator to illustrate our methodology in this section. The nonlocal diffusion operator is given by

[TABLE]

See more discussions of details in [6]. Here we assume that the kernel $\omega$ is given by

[TABLE]

where the kernel $\gamma$ is a nonnegative function supported on $[-1,1]$ with $\int_{-1}^{1}s^{2}\gamma(|s|)ds=1$. $C(\bm{x},\bm{y})$ is the diffusion coefficient and $\delta(\bm{x},\bm{y})$ is the horizon function that satisfy

[TABLE]

Here we keep the dependence of $C$ and $\delta$ on $\bm{x}$, $\bm{y}$ so that it could model the inhomogeneous medium, see related works [27, 30, 7].

The common practice in the peridynamics simulation (see [14]) takes ${\gamma}$ to be

[TABLE]

We distinguish between the two cases:

1. (a)

If $C$ and $\delta$ only depends on $x$, namely $C(\bm{x},\bm{y})=C(\bm{x})$ and $\delta(\bm{x},\bm{y})=\delta(\bm{x})$, then (2.3) is the nonlocal diffusion operator of non-divergence type. As $\delta_{1}\to 0$, we have

[TABLE] 2. (b)

If $C$ and $\delta$ are symmetric, namely $C(\bm{x},\bm{y})=C(\bm{y},\bm{x})$ and $\delta(\bm{x},\bm{y})=\delta(\bm{y},\bm{x})$, then (2.3) is the nonlocal diffusion operator of divergence type. As $\delta_{1}\to 0$, we have

[TABLE]

The FFT-based methods mentioned earlier are restrictive in the application to heterogeneous media. Here we seek algorithms that can be applied to the more general cases given by (2.4). We use the recently proposed hierarchical interpolative factorization (HIF) [18] method to factorize the dense matrix $A$ obtained from discretizing the nonlocal operator $\mathcal{L}$. The factorization of the matrix can be used to rapidly apply both $A$ and $A^{-1}$. Therefore it can serve as a direct solver or be stored and reused for iterative methods or time-dependent problems. The HIF takes advantage of the low-rank behavior of the off-diagonal entries of the dense matrix, so the far field behavior of the kernel $\omega$ is the central factor on the performance of the the algorithm. We will illustrate in this section that the performance of the HIF algorithm improves greatly as the smoothness of the kernel away from zero enhances. All computations are performed in MATLAB R2014b on a single core of a 1.1GHz Intel Core M CPU on a 64-bit Mac laptop.

2.1 Nonlocal diffusion in homogeneous media

In this section, we solve the the nonlocal problem on the interval $\Omega=(0,1)$ with $C(x,y)\equiv 1$ and $\delta(x,y)\equiv\delta$ with $\delta$ being the horizon radius of the homogeneous medium. The matrix $A$ is obtained form the discretization of the nonlocal operator $\mathcal{L}$, which is is based on the asymptotically compatible schemes in [29], which ensures uniform discretization error independent of the horizon function $\delta$. In Table 1, the storage and the time of matrix-vector multiplication are listed for the original matrix $A$ as well as its compression after using the HIF method. The original matrix $A$ is stored using the sparse matrix representation with the MATLAB built-in matrix-vector multiplication applied to the test. The kernel $\gamma(s)$ has $O(1/|s|)$ singularity at [math] and is assumed to be smooth expect at origin and $\pm 1$ ($\omega$ has singularity at [math] and $\pm\delta$). The number of discretization nodes $N$ is fixed to be $2048$ in Table 1. The relative precision parameter for the HIF method is set to be $1E-8$.

The left column of Table 1 lists the different types of regularity (at $\pm 1$) of the kernel $\gamma(s)$ and the different choices of the horizon parameter $\delta$. The discontinuous kernel (denoted as $C^{-1}$) is the one that is most often used in peridynamics simulations given by

[TABLE]

The more regular $C^{k}$ ($k=0,1,2,3$) kernels are obtained by subtracting $\gamma$ with polynomials $p^{2k}$ ($k=0,1,2,3$) of degree $2k$. The details of constructing these polynomials will be given in Section 3. When $\delta\geq 1$, the matrix $A$ is fully dense and it does not see the tail of the kernel. As a result, the HIF works well for all the types of kernels. When $\delta$ gets smaller, the regularity of the kernel $\gamma(s)$ at $\pm 1$ plays a decisive role on the performance of the HIF. We see that when $\gamma(s)$ is discontinuous at $\pm 1$, no compression of the matrix $A$ is observed using the HIF. The storage taken and computation time decrease with increasing regularity of $\gamma(s)$ from discontinuous ($C^{-1}$) all the way to three times continuously differentiable ($C^{3}$). In the case of $\delta=1/4$, the critical improvement of the computational efficiency happens between $C^{0}$ kernel functions and $C^{1}$ kernel functions as shown in Table 1. We remark that the location of critical regularity changes with the precision parameter $\epsilon$. When the HIF algorithm runs with higher precision, the higher critical regularity is observed; and when the HIF algorithm runs with lower precision, the lower critical regularity is observed. Similar patterns are also observed in the experiments that will be discussed in Section 2.2 and 2.3. When $\delta$ keeps decreasing, the effect of different regularity of the kernel $\gamma(s)$ will diminish. The extreme case is that $\delta$ equals the mesh size $h$ when the matrix $A$ is in fact tridiagonal. The case of $\delta=h$ is listed in Table 1 to compare the overhead of using the HIF algorithm with the sparse matrix. The overhead of using the HIF algorithm in this case is due to the fact that HIF does not build into itself the ability to detect the horizon and ignore the zeros outside. It will be future work to explore this possibility depending on application.

To account for the complexity of the HIF, we provide the computation time for different discretization nodes $N$ in Table 2. The horizon parameter $\delta$ in this table is fixed to be $1/4$. The scaling results for the choice of the $C^{0}$ kernel is compared with the choice of the $C^{3}$ kernel. The complexity of the matrix-vector multiplication for the original matrix $A$ without compression is approximately $O(N^{2})$ in both cases. On the other hand, the HIF code behaves quite differently for the $C^{0}$ kernel and the $C^{3}$ kernel. The observed complexity of HIF with the $C^{0}$ kernel is still , while the complexity of HIF with the $C^{3}$ kernel is just around $O(N)$.

The simulations show effectiveness of the HIF fast algorithm with the improving far field regularity of nonlocal interaction kernels. Finally, we remark that the cost of constructing the factorization is in general larger than the cost of matrix-vector multiplication. This is acceptable since one only need to construct the factorization once and reuse it for an iterative method or in time-dependent problems.

2.2 Nonlocal diffusion in heterogeneous media

We test the effectiveness of the HIF algorithm for nonlocal diffusion in heterogeneous media in this section. The kernel $\omega$ is defined as (2.4), where $C(x,y)\equiv 1$ and $\delta(x,y)=\delta(x)$. Each point $x$ has its own horizon radius $\delta(x)$. In this case, $\mathcal{L}$ is a nonlocal diffusion operator of non-divergence type. The matrix $A$ obtained form the discretization of the nonlocal operator $\mathcal{L}$ is now non-symmetric, in which case the FFT-based fast algorithm will not work. In Table 3, the storage and the time of matrix-vector multiplication are listed for the original matrix $A$ as well as its compression after using the HIF algorithm.

In the numerical experiments we use $\delta(x)=\delta_{0}(1+e^{-20(x-0.5)^{2}})$ for $x\in\Omega=(0,1)$ so that $\delta(x)$ is a function between $\delta_{0}$ and $2\delta_{0}$, as shown in Fig 1.

The kernel $\gamma(s)$ has $O(1/|s|)$ singularity at [math] and is assumed to be smooth expect at origin and $\pm 1$. The left column of Table 3 lists the types of regularity (at [math] and $\pm 1$) of the kernel $\gamma$ and the horizon parameter $\delta_{0}$. The HIF algorithm for nonlocal diffusion in heterogeneous media works similarly as the homogeneous cases. In particular, we tailored Table 3 to the case where $\delta_{0}=1/4$ because for the extreme cases where $\delta_{0}=1$ or $\delta_{0}=h$, the variation of regularity of the kernel does not result in any change in the computational cost, which is exactly the same behavior observed in Table 1.

To account for the complexity, we provide the computation time for different numbers of discretization nodes $N$ in Table 4. The horizon parameter $\delta_{0}$ is fixed to be $1/4$. Two types of kernels – the $C^{0}$ kernel and the $C^{3}$ kernel – are used in the simulation. The complexity of the matrix-vector multiplication for the original matrix $A$ without compression is $O(N^{2})$ as expected. Similar to the homogeneous case, the observed complexity of HIF in the heterogeneous case is again $O(N^{2})$ for the $C^{0}$ kernel and $O(N)$ for the $C^{3}$ kernel. The numerical examples in Sections 2.1 and 2.2 indicate that the HIF has the same effectiveness for models in heterogeneous media as it has for models in homogeneous media if the kernel functions are sufficiently smooth away from origin.

2.3 Two-dimensional test

In this section, we perform 2d tests on the domain $\Omega=(0,1)\times(0,1)$ to show the effectiveness of the HIF code. We choose $C(x,y)\equiv 1$ and $\delta(x,y)\equiv\delta$. The homogeneous medium has the same the horizon radius $\delta$ for every point. The kernel $\gamma(s)$ is the same as before and has $O(1/|s|)$ singularity at [math] and is assumed to be smooth expect at origin and $\pm 1$.

In Table 5, the storage and the time of matrix-vector multiplication are listed for the original matrix $A$ as well as its compression after using the HIF method. The original matrix $A$ is again stored using the sparse matrix representation and the MATLAB built-in matrix-vector multiplication is used to record the computation time in the right column of Table 5. The relative precision is chosen to be $1E-3$ in the 2d test. We remark that in the case of $\delta=1/2$, the critical improvement of the computational efficiency happens between discontinuous kernel functions and $C^{0}$ kernel functions as a result of the choice of the precision parameter $\epsilon$. Similarly as in 1d, higher critical regularity is observed if $\epsilon$ is getting smaller.

Table 6 shows the storage and time for different value of $N$ with horizon parameter fixed to be $1/2$. The discontinuous ($C^{-1}$) kernel and the $C^{3}$ kernel are used in the simulation. Here we remark that we choose the $C^{-1}$ kernel to be compared with the $C^{3}$ kernel because $C^{-1}$ is right before the critical improvement of efficiency happens with the increase of regularity as seen in Table 5. This is a result of the particular choice of relative precision $\epsilon=1E-3$. As $\epsilon$ becomes smaller, higher regularity is needed in order for the critical improvement of efficiency to happen. We observe that the complexity of the matrix-vector multiplication for the original matrix $A$ without compression is approximately $O(N^{2})$. HIF also gives $O(N^{2})$ computation time for matrix-vector multiplication with the discontinuous kernel being used, while in the case of $C^{3}$ kernel, HIF performs much better and the observed complexity for matrix-vector multiplication is about $O(N\log(N))$.

To sum up, the simulations in Sections 2.1, 2.2 and 2.3 show that the popular choice of the kernel $\gamma$ in (2.6) is problematic in the interest of fast computation since it has non-smooth truncation at $\pm 1$. On the other hand, when the regularity of the kernel gets better the HIF fast algorithm becomes more effective. Those observations motivate us to propose the splitting method that deals kernels with non-smooth truncation in the simulations of the peridynamic model.

3 The splitting of singularities

Popular choices of kernels in the peridynamic simulations contain two types of singularities: the singularity at origin and the non-smooth truncation at the boundary $\partial\mathcal{H}_{\bm{x}}$. The first type singularity is handled by the HIF algorithm through the idea of far-field compression. In 1d, the compression is performed through a block matrix with small block sizes near the diagonal. The algorithm is thus effective when the kernel is sufficiently smooth away from origin, as we have seen in the previous numerical examples. Next, we propose to split the two types of singularities and treat them respectively. In Fig 2, the left plot is the kernel $\gamma(|s|)=\frac{1}{|s|}\chi{(|s|<1)}$. The kernel splits into two parts, the middle plot which has singularity at zero but $C^{3}$ away from zero, and the right plot which is discontinuous at $s=\pm 1$ but smooth in the interior. In practice, the kernel $\gamma(|s|)$ is written into

[TABLE]

where $\kappa(|s|)$ is $C^{K}$ regular away from zero, and $p^{2K}$ is an even polynomial of degree $2K$ that matches $\gamma$ and its derivatives at $|s|=1$ up to the $K$-th order. For the choice of $\gamma(|s|)=\frac{1}{|s|}\chi{(|s|<1)}$, $p^{2K}$ is explicitly given by the following formulas for $K=0,1,2,3$.

[TABLE]

As suggested by Section 2, the nonlocal operator corresponding to the kernel $\kappa(s)$ can be treated by the HIF fast algorithm effectively. So we only need to deal with the second part corresponding to the kernel $p^{2K}(s)\chi(|s|<1)$. Taking $\Omega=[0,1]^{d}$ in the $d$-dimensional Euclidean space, we consider the following nonlocal diffusion problem with Dirichlet type boundary

[TABLE]

where $\mathcal{H}_{\bm{x}}$ is some neighborhood of $\bm{x}$, and $P$ is given by

[TABLE]

For convenience of illustration, we take $C(\bm{x},\bm{y})\equiv C(\bm{x})$, $\delta(\bm{x},\bm{y})\equiv\delta(\bm{x})$, $K=0$ and $p^{0}\equiv 1$ in the next. The generalization to $p^{2K}$ for $K>0$ is straightforward. Taking into account of the boundary condition in (3.7), the nonlocal operator is to be evaluated at each $x$:

[TABLE]

We therefore need only to find a fast algorithm compute $\int_{\mathcal{H}_{\bm{x}}\cap\Omega}u(\bm{y})d\bm{y}$. Our algorithm starts by the hierarchical decomposition of the space $\Omega$, which is similar to the one in standard FMM [12]. Assume that the domain $\Omega$ is uniformly $N=n^{d}$ be the total number of discretization nodes. Let $L=\log_{2}(n)$, then from level $1$ to $L$, the computation domain $\Omega=[0,1]^{d}$ is hierarchically subdivided into panels. Each panel in the $l$-th level can be represented by one of the cubes $\prod_{j=1}^{d}[\frac{k_{j}}{2^{l}},\frac{k_{j}+1}{2^{l}}]$ , $k_{j}=0,1,\cdots,2^{l}-1$.

The hierarchical subdivision of $\Omega$ forms a tree structure. Fig 3 is the quadtree formed by the subdivision of the two-dimensional domain $[0,1]^{2}$. The blue dot in Fig 3 represents the original domain $[0,1]^{2}$, and it has four children $\{[\frac{i}{2},\frac{i+1}{2}]\times[\frac{j}{2},\frac{j+1}{2}]\}_{i,j=0}^{1}$ represented by the red dots in level $1$. The subdivision of the $4$ domains in level $1$ results in $16$ subdomains in level $2$ represented by the black dots. The hierarchical subdivision is performed until reaching the $L$-th level. Similarly, in 1d a binary tree will be formed by the subdivision of $[0,1]$ and in 3d an octree will be formed by the subdivision of $[0,1]^{3}$.

We now present the algorithm. It contains the following $3$ steps that will be explained in detail.

Initialization. Traverse the tree from top to bottom to obtain the decomposition of $\mathcal{H}_{\bm{x}}\cap\Omega$ for each grid point $\bm{x}$. 2. 2.

Traverse the tree from bottom to top to obtain the partial sums $\sum_{\bm{x}_{k}\in Q_{l}^{j}}u(\bm{x}_{k})$ for each panel $Q_{l}^{j}$ in the $l$-th level. 3. 3.

*Approximate the integral $\int_{\mathcal{H}_{\bm{x}}\cap\Omega}u(\bm{y})d\bm{y}$ by using (1) and (2). *

Step 1. Initialization step. For each $\bm{x}$ on the grid point, decompose the domain of integration $\mathcal{H}_{\bm{x}}\cap\Omega$ using the panels in $L_{l}$ ($l=0,1,\cdots,L$). The decomposition of $\mathcal{H}_{\bm{x}}\cap\Omega$ for each $\bm{x}$ is stored for reuse. More specifically, we traverse the tree from top to bottom to determine whether a panel $Q_{l}^{j}$ from level $l$ is included or not, and then express $\mathcal{H}_{\bm{x}}\cap\Omega=\bigcup_{l=0}^{L}\bigcup_{j\in\mathcal{I}_{l}(\bm{x})}Q_{l}^{j}$, where $\mathcal{I}_{l}(\bm{x})$ is the set of indices at the $l$-th to be included in the decomposition. We may do this recursively by calling the function recur($\Omega$), where recur is defined in Algorithm 1 below.

Let us now discuss on the complexity of performing the recur function when $\mathcal{H}_{\bm{x}}$ is a $d$-dimensional ball centered at $\bm{x}$. The main cost of running the recur function comes from the intersection test of balls and cubes. To test whether a cube $Q$ is fully contained in a ball, we only need to test whether the total of $2^{d}$ corners of $Q$ are all contained in the ball. This requires a total of $O(2^{d})$ operations, where the constant in $O(2^{d})$ is independent of the dimension $d$. Now to test whether a cube $Q$ has nontrivial intersection with a ball, we use the algorithm given in [2], which is an $O(d)$ algorithm that determines whether an axis-aligned bounding box (AABB) intersects a ball. See the function defined in Algorithm 2, where the inputs are three vectors $\textbf{\text{Bmin}},\textbf{\text{Bmax}},\textbf{\text{C}}$ and a positive number $r$. The vector $\textbf{\text{Bmin}}\in\mathbb{R}^{d}$ stores the minima of the AABB for each axis, $\textbf{\text{Bmax}}\in\mathbb{R}^{d}$ stores the maxima of the AABB for each axis, and $\textbf{\text{C}}\in\mathbb{R}^{d}$ and $r$ are the center and radius of the ball resepctively. Combining the above discussions, we know that each time calling the recur function requires $C\cdot 2^{d}$ operations, where $C$ is independent of dimension.

Fig 4 shows an example of decomposing a circular region in 2d. Now the total complexity of Step 1 is obviously equal to $C\cdot 2^{d}$ times the total number of times the recur function is called. Now for each discretization node $\bm{x}_{i}$ ($i=1,2,\cdots,N$), let $\mathscr{N}_{d}(\bm{x}_{i})$ denote the number of times the recur function is called in order to decompose the domain $\mathcal{H}_{\bm{x}_{i}}\cap\Omega$. The total number of times the recur function is called is then given by $\sum_{i=1}^{N}\mathscr{N}_{d}(\bm{x}_{i})$. Now observe that $\mathscr{N}_{d}(\bm{x}_{i})$ is also equal to the number of panels that intersect with the boundary $\partial\mathcal{H}_{\bm{x}_{i}}$ non-trivially, namely

[TABLE]

where $\mathscr{N}^{l}_{d}(\bm{x}_{i})$ denotes the number of panels in the $l$-th level that intersect with $\partial\mathcal{H}_{\bm{x}_{i}}$ non-trivially. Assume that for each point $\bm{x}_{i}$, the set $\partial\mathcal{H}_{\bm{x}_{i}}$ has uniform bounded mean curvature, then $\mathscr{N}^{l}_{d}(\bm{x}_{i})$ can be estimated by the following formula

[TABLE]

To sum up, the total complexity of Step 1, which equals $C\cdot 2^{d}\sum_{i=1}^{N}\sum_{l=0}^{L}\mathscr{N}^{l}_{d}(\bm{x}_{i})$ is given by

[TABLE]

where the dimensional dependent constant $C_{d}$ is of order $O(2^{d})$.

Step 2. Compute the the partial sums for each tree node . For a given vector $\{u(\bm{x}_{i})\}_{i=1}^{N}$, we assign the value $u(\bm{x}_{i})$ to a leaf node. Then traverse the tree bottom to top, we assign each parent node the value of the sum of all its children. The completion of this step gives us the partial sums $\sum_{\bm{x}_{k}\in Q_{l}^{j}}u(\bm{x}_{k})$ for every panel $Q_{l}^{j}$. The partial sums obtained can then be used to approximate the integral $\int_{Q_{l}^{j}}u(\bm{y})d\bm{y}$. We remark that if the kernel function given by (3.8) is generated from the polynomial $p^{2K}$ with $K>0$, we then need to compute the partial sums in the form of $\sum_{x_{k}\in Q_{l}^{j}}x_{k}^{m}u(x_{k})$ for all $m=0,1,\cdots 2K$, so that they can be used to approximate the integral $\int_{Q_{l}^{j}}P(\bm{x},\bm{y})u(\bm{y})d\bm{y}$.

Fig 5 exemplifies the way partial sums are computed for a two-dimensional problem. Each of the black dots on the left figure represents a leaf node that contains the partial sum of the function in the little black box that includes it. Then the partial sums are passed to the coarser level where each parent node presented by the red dots stores the sum of the values from its four children. The total sum of the function on the whole domain is stored in the top tree node presented by the blue dot in the right figure in Fig 5. Therefore, the computational cost of Step 2 depends only on the number of discretization nodes $N$ and the degree $2K$ of the polynomial $p^{2K}$. To sum up, the complexity of Step 2 is given by $O(KN)$, where the constant in $O(KN)$ is independent of dimension.

Step 3. For each discretization node $\bm{x}_{i}$ ($i=1,2\cdots N$), use the decomposition of $\mathcal{H}_{\bm{x}_{i}}\cap\Omega$ obtained in Step 1 and the partial sums obtained in step 2 to approximate the integral

[TABLE]

Finally, the value of $\mathcal{L}_{\delta}^{P}u(\bm{x}_{i})$ is obtained through the relation (3.9). We again remark that in the case of $K>0$ for the polynomial $p^{2K}$, all the partial sums of the form $\sum_{x_{k}\in Q_{l}^{j}}x_{k}^{m}u(x_{k})(m=0,1,\cdots,2K)$ will be needed to approximate the desired integral. With the precomputed sums for the approximate integrals of the form $\int_{Q_{l}^{j}}u(\bm{y})d\bm{y}$ in Step 2, the number of summations we need to take in (3.11) is essentially equal to $\mathscr{N}_{d}(\bm{x}_{i})$ defined in Step 1. By the calculations in Step 1 and taking into the account of the general case that $p^{2K}$ being a polynomial of degree $2K>0$, complexity of Step 3 is then given by the following formulas,

[TABLE]

and the constants in the above estimates are independent of dimension.

3.1 Numerical tests

We perform numerical tests based on the proposed algorithm.

Example 1. In this example, we take $\Omega=[0,1]^{2}$ in the 2-dimensional space. Assume $k=0$, $p^{0}(x)\equiv 1$, $\delta(\bm{x})=\frac{1}{2}$ and $C(\bm{x})=1$. The nonlocal operator to be evaluated is given by

[TABLE]

Table 7 contains the experimental data of the new algorithm in comparison with two versions of matrix-vector multiplication for the original matrix. The scaling results are shown in Fig. 6. Since the initialization step (Step 1) is only needed to be performed once in an iterative method or time-dependent problem, the computation time recorded for the new algorithm is given by the cost for Step 2 and Step 3 combined together. The computation time of the new algorithm grows at the rate $O(N^{1.5})$, which verifies that the the formula given in (3.10) for $d=2$. In comparison, the two versions of matrix-vector multiplication for the original matrix shows $O(N^{2})$ complexity. The first version (Original - 1) is performed with the MATLAB built-in matrix-vector multiplication. Because of the nature of MATLAB, the second version (Original -2), which goes through the evaluation of the nonlocal operator at each node $\bm{x}_{i}$ $(i=1,2,\cdots N)$ by using a for-loop, is also listed for a fair comparison with the new algorithm.

Example 2. To show that the algorithm can be applied to more generalized cases, we perform the numerical test on the 1d domain $\Omega=[0,1]$ and the kernel function $P(x,y)$ is given by (3.8) with $C(x,y)=1$, $\delta(x,y)=\delta(x)$ as shown in Fig 1. We take $K=3$ and $p^{6}$ a polynomial of degree $6$. Now in Step 2 of the algorithm, we not only need to compute the partial sums $\sum_{x_{k}\in Q_{l}^{j}}u(x_{k})$, but also $\sum_{x_{k}\in Q_{l}^{j}}x_{k}^{m}u(x_{k})$ for all $m=1,2,\cdots 6$ so that they can be used to approximate $\int_{\mathcal{H}_{x}\cap\Omega}P(x,y)u(y)dy$.

The computation time for the new algorithm is given in Table 8 in comparison with two versions of matrix-vector multiplication for the original matrix, with scaling results shown in Fig. 7. The computation time of the new algorithm grows at the rate which is slightly higher than $O(N)$, which relates well to our theoretical estimate of $O(N\log(N))$ for $d=1$.

4 Conclusion

In this work, we develop and investigate fast algorithms for a class of nonlocal models featured by finite domain of nonlocal interactions. Important applications include the computation of the peridynamics model of solid mechanics and models of nonstandard diffusion. The main message is that the far field behavior of the nonlocal interaction kernel is decisive in order for the fast computation of nonlocal models. We use numerical tests for the nonlocal diffusion model as a proof of concept, and demonstrate that better regularity of the kernel away from zero is better in the interest of fast computation. This is reasonable since most fast algorithms that can be applied to non-translation invariant operators and general meshes are based on low rank approximation of the far field interactions. Based on this observation, the kernels that are most often used in the peridynamics simulations fail to work under existing fast algorithms since they all have non-smooth truncations at a finite distance. We propose a remedy for such kernels, and that is the splitting of the singularities at origin and at the finite distance and treat them separately. The splitting is done by constructing polynomials that match the derivatives of the original kernels at finite distances. The second type of singularity after splitting is on larger sets than just point discontinuity. We propose a new FMM like fast algorithm to deal with such singularities which exploits the smoothness of the kernels away from the truncation. Such idea that uses low rank approximations of kernels for different regions is also presented in [8] for a completely different problem. Note that we only focus on the fast matrix-vector multiplication in this work, and it can be used in a Krylov subspace method if computing the inverse of a matrix is the interest. Notice that the kernel function that uses to regularize the original kernel typically is not small, as indicate in the bottom right picture in Fig. 2. If the original kernel is a small perturbation of a kernel that have smooth decay in the far field, then the inverse of the corresponding matrix could be simply computed using the series expansion under the smallest assumption. The goal of this paper is to show that fast algorithms for nonlocal operators with finite interaction distances are possible. Even if the techniques developed here substantially reduces the computational cost, there is room for improvements from future research. It would be natural to pursue further research as nonlocal modeling becomes more common. The essentially optimal computational complexity of $O(N\log N)$ is for 1d and in higher dimensions there is algebraic complexity higher than linear. For further complexity reduction in higher dimensions a better representation of precomputed quantities near the horizon would be required and it would be natural to use different geometries than boxes, for example, curvelets [5] in 2d deals efficiently with line discontinuities.

Based on our findings, several more points have the value of further study. First, our work offers important messages to nonlocal modeling from the perspective of scientific computing. It is worthwhile to explore in the future whether the kernels with smooth decays at far field can be used in replace of non-smooth kernels in the nonlocal models without compromising the physics. Second, in order to understand the properties of the algorithms, we separated the global HIF technique from the handling of discontinuous at the horizon by using a simple splitting. The more interesting question on the optimal way to split the singularities could be asked in the future. Moreover, it will be more desirable to build a tighter coupling of the singularities into the original HIF technique, and the question on detecting the zero values of the kernel could be asked as mentioned earlier. Finally, we note that our case study is restricted to nonlocal diffusion type problem. It is reasonable that for the peridynamics system of equations similar studies can be done.

Acknowledgements

The first author want to thank Lexing Ying and Yuwei Fan for providing with helpful discussions and HIF codes during her visit to Stanford University. The authors also thank Qiang Du for helpful discussions on the subject. The Oden Institute is acknowledged for its support. The research of Bjorn Engquist is supported in part by the U.S. NSF grant DMS-1620396. The research of Xiaochuan Tian is supported in part by the U.S. NSF grant DMS-1819233. We thank the referees for the helpful comments to improve this work.