---
author:
- 'Christian Vollmann[^1]'
- Volker Schulz
bibliography:
- 'literatur.bib'
title: '**Exploiting multilevel Toeplitz structures in high dimensional nonlocal diffusion**'
---

**Abstract.** We present a finite element implementation for the steady-state nonlocal Dirichlet problem with homogeneous volume constraints. Here, the nonlocal diffusion operator is defined as integral operator characterized by a certain kernel function. We assume that the domain is an arbitrary $d$-dimensional hyperrectangle and the kernel is translation invariant. Under these assumptions, we carefully analyze the structure of the stiffness matrix resulting from a continuous Galerkin method with multilinear elements and exploit this structure in order to cope with the curse of dimensionality associated to nonlocal problems. For the purpose of illustration we choose a particular kernel, which is related to space-fractional diffusion and present numerical results in 1d, 2d and for the first time also in 3d.\
 \
**Keywords.** Nonlocal diffusion, finite element method, translation invariant kernel, multilevel Toeplitz, fractional diffusion.

Introduction
============

The field of nonlocal operators attracts increasing attention from the mathematical society. This is due to the steadily growing pool of applications where nonlocal models are in use; including e.g., image processing [@imageproc1; @imageproc2], machine learning [@machinelearning], peridynamics [@peridym1; @peridym2], fractional diffusion [@fractlapl] or nonlocal Dirichlet Forms [@nonlocaldirichletform] and jump processes [@jumpprocess]. In contrast to local diffusion problems, interactions can occur at distance in the nonlocal case. This relies on the definition of the nonlocal diffusion operator $-{\mathcal{L}}$, which acts on a function ${u\colon {\mathbb{R}}^d \to {\mathbb{R}}}$ by $$-{\mathcal{L}}u(x) {\coloneqq}2\int_{{\mathbb{R}}^d} (u(x)-u(y)){\gamma}(x,y)dy,$$ where ${\gamma}\colon {\mathbb{R}}^d \times {\mathbb{R}}^d \to {\mathbb{R}}$ is a nonnegative and symmetric function characterizing the precise nonlocal diffusion. In this paper we are interested in the *steady-state nonlocal Dirichlet problem with volume constraints* given by $$\begin{aligned}
-{\mathcal{L}}u(x) = f(x) ~~~&(x\in \Omega) ,\nonumber\\
u(x) = g(x)     ~~~&(x\in \Omega_I), \label{prob}\end{aligned}$$ where $\Omega \subset {\mathbb{R}}^d$ is a bounded domain. Here, the constraints are defined on a volume $\Omega_I$, the so called interaction domain, which is disjoint from $\Omega$.

The recently developed vector calculus by Gunzburger et al. [@nonlocalveccal] builds a theoretic foundation for the description of these nonlocal diffusion phenomena. In particular, this framework allows us to consider finite-dimensional approximations using Galerkin methods similar to the analysis of (local) partial differential equations. However, in contrast to local finite element problems, here we are faced with two basic difficulties. On the one hand, the assembling procedure may require sophisticated numerical integration tools in order to cope with possible singularities of the kernel function. On the other hand, discretizing nonlocal problems leads to densely populated systems. The latter tremendously affects the solving procedure, especially in higher dimensions. Thus, numerical implementations are challenging and in order to lift the concept of nonlocal diffusion from a theoretical standpoint to an applicable approach in practice, the development of efficient algorithms, which go beyond preliminary cases in 1d and 2d, is essential in this context.

In recent works, several approaches for discretizing problem (\[prob\]) have been presented. We want to mention for instance a 1d finite element code by D’Elia and Gunzburger [@fractlapl], where the fractional kernel ${\gamma}(x,y) =\frac{c_{d,s}}{2||y-x||^{d+2s}}$ has been used, but which can easily be extended to general (singular) kernels. Also for the fractional kernel, Acosta, Bersetche and Borthagaray [@acosta] developed a finite element implementation for the 2d case. In general, a lot of work has been done for the discretization of fractional diffusion problems and fractional derivatives in various definitions (mainly via finite difference schemes) not only for 1d and 2d [@wangToeplitz; @wang2d], but also for the 3d case [@wang3d]. However, to the best of the authors’ knowledge, for general kernel functions 3d finite element implementations for problem (\[prob\]) are not yet available.

In this paper we study a finite element approximation for problem (\[prob\]) on an arbitrary $d$-dimensional hyperrectangle (parallel to the axis) for translation invariant kernel functions. More precisely, we analyze from a computational point of view a continuous Galerkin discretization with multilinear elements of the following setting:

-   We set $\Omega {\coloneqq}\prod_{i=0}^{d-1}[a_i,b_i]$, where $[a_i, b_i]$ are compact intervals on ${\mathbb{R}}$.

-   We assume that the kernel ${\gamma}$ is *translation invariant*, such that $${\gamma}(b+Qx, b+Qy) = {\gamma}(x,y)$$ for all $b\in{\mathbb{R}}^d$ and $Q \in {\mathbb{R}}^{d \times d}$ orthonormal.

As a consequence, these structural assumptions on the underlying problem are reflected in the stiffness matrix; we obtain a $d$-level Toeplitz matrix, which has two crucial advantages. On the one hand, we only need to assemble (and store) the first row (or column) of the stiffness matrix. On the other hand, we can benefit from an efficient implementation of the matrix-vector product for solving the linear system. This result is presented in Theorem \[multilevel\_toeplitz\] and is crucial for this work, since it finally enables us to solve the discretized system in an affordable way. For illustrative purposes we choose the fractional kernel and exploit a third assumption on the interaction horizon for simplifying the implementation:

-   We assume that interactions only occur at a certain distance $R$, which we assume to be larger or equal the diameter of the domain $\Omega$, such that $\Omega \subset B_R(x)$ for all $x \in \Omega$.

The paper is organized as follows. In Section \[sec:theory\] we cite the basic results about existence and uniqueness of weak solutions and finite-dimensional approximations. In Section \[sec:femsetting\] we give details about the precise finite element setting and proof our main result, that the stiffness matrix is multilevel Toeplitz. In Sections \[sec:assembling\] and \[sec:solving\] we explain in detail the implementation of the assembling and solving procedure, respectively. In Section \[sec:numres\] we round off these considerations by presenting numerical results with application to space-fractional diffusion.

Nonlocal diffusion problems {#sec:theory}
===========================

We review the relevant aspects of nonlocal diffusion problems as they are introduced in [@wellposedness] which constitute the theoretic fundamentals of this work.

Let $\Omega \subset {\mathbb{R}}^d$ be a bounded domain with piecewise smooth boundary. Further let ${\gamma}\colon {\mathbb{R}}^d \times {\mathbb{R}}^d \to {\mathbb{R}}$ be a nonnegative and symmetric function (i.e., ${\gamma}(x,y) = {\gamma}(y,x) \geq 0$), which we refer to as *kernel*. Then we define the action of the *nonlocal diffusion operator* $-{\mathcal{L}}{\coloneqq}-{\mathcal{L}}_{\gamma}$ on a function $u \colon {\mathbb{R}}^d \to {\mathbb{R}}$ by $$-{\mathcal{L}}u(x) {\coloneqq}2\int_{{\mathbb{R}}^d} (u(x)-u(y)){\gamma}(x,y)dy$$ for $x \in \Omega$. In addition to that, we assume that there exists a constant ${\gamma}_0 > 0$ and a finite *interaction horizon* or *radius* $R > 0$ such that for all $x \in \Omega$ we have $$\begin{aligned}
{\gamma}(x,y) \geq 0 ~~~&\forall~ y \in B_R(x), \nonumber\\
{\gamma}(x,y) \geq {\gamma}_0 > 0 ~~~&\forall ~y \in B_{R/2}(x),\nonumber\\
{\gamma}(x,y) = 0   ~~~&\forall ~y \in (B_{R}(x))^c, \label{gen_assumption}\end{aligned}$$ where $B_R(x) {\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}y \in {\mathbb{R}}^d \colon ||x-y||_2< R\rk$. Then we define the *interaction domain* by $$\Omega_I {\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}y \in \Omega^c\colon \exists x\in \Omega\colon{\gamma}(x,y) \neq 0 \rk$$ and finally introduce the *steady-state nonlocal Dirichlet problem with volume constraints* as $$\begin{aligned}
-{\mathcal{L}}u(x) = f(x) ~~~&(x\in \Omega) ,\nonumber\\
u(x) = g(x)     ~~~&(x\in \Omega_I), $$ where $f\colon \Omega \to {\mathbb{R}}$ is called the *source* and $g: \Omega_I \to {\mathbb{R}}$ specifies the Dirichlet volume constraints. In the remainder of this paper we assume $g\equiv 0$.

Weak formulation
----------------

For the purpose of constructing a finite element framework for nonlocal diffusion problems, we introduce the concept of weak solutions as it is presented in [@wellposedness].

We define the bilinear form $$\begin{aligned}
a(u,v) {\coloneqq}& \int_{{\Omega \cup \Omega_I}} v(-{\mathcal{L}}u )  dx,\end{aligned}$$ and the associated linear functional $$\begin{aligned}
\ell(v) {\coloneqq}\int_{\Omega} fvdx.\end{aligned}$$ By establishing a nonlocal vector calculus it is shown in [@wellposedness], that the following equality holds: $$\begin{aligned}
a(u,v)= &\int_{{\Omega \cup \Omega_I}} \int_{{\Omega \cup \Omega_I}\cap B_R(x)} (u(y)-u(x))(v(y)-v(x)){\gamma}(x,y)  dydx.\end{aligned}$$ This implies that $a$ is symmetric and nonnegative, or equivalently, the linear nonlocal diffusion operator $-{\mathcal{L}}$ is self-adjoint with respect to the $L^2$-product and nonnegative. Furthermore, we define the *nonlocal energy space* $$V({\Omega \cup \Omega_I}) {\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}u \in L^2({\Omega \cup \Omega_I}) \colon |||u||| <\infty \rk,$$ where $ |||u|||{\coloneqq}\sqrt{\tfrac{1}{2}a(u,u)}$ and the *nonlocal constrained energy space* $$V_c({\Omega \cup \Omega_I}) {\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}u \in L^2({\Omega \cup \Omega_I}) \colon |||u||| <\infty  ~\text{and}~ u_{|\Omega_I} \equiv 0 ~\text{a.e.}\rk.$$ We note that $||| \cdot |||$ constitutes a semi-norm on $V({\Omega \cup \Omega_I}) $ and due to the volume constraints a norm on $V_c({\Omega \cup \Omega_I}) $. With these preparations at hand, a weak formulation of (\[prob\]) can be formulated as $$\begin{aligned}
\text{\textit{Find $u \in V_c({\Omega \cup \Omega_I}) $ such that }} a(u,\cdot) \equiv \ell(\cdot) ~ \text{\textit{on}}~ V_c({\Omega \cup \Omega_I}) \label{weakform}.\end{aligned}$$

In order to make statements about the existence and uniqueness of weak solutions we have to further specify the kernel. In [@wellposedness] the authors consider, among others, a certain class of kernel functions, on which we will focus in the remainder of this section and in our numerical experiments. More precisely, we require that there exists a fraction $s\in (0,1)$ and constants ${\gamma}_1,{\gamma}_2 > 0$, such that for all $x \in \Omega$ it holds that $$\begin{aligned}
{\gamma}(x,y)||y-x||^{d+2s} \in [{\gamma}_1,{\gamma}_2] ~~~&\forall ~y \in B_R(x) \label{case_1}.\end{aligned}$$ Then it is shown in [@wellposedness] that the nonlocal energy space $V({\Omega \cup \Omega_I})$ is equivalent to the fractional-order Sobolev space $$\begin{aligned}
H^s({\Omega \cup \Omega_I}) {\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}u\colon ||u||_{H^s({\Omega \cup \Omega_I})} {\coloneqq}||u||_{L^2({\Omega \cup \Omega_I})} +  |u|_{H^s({\Omega \cup \Omega_I})} < \infty   \rk,\end{aligned}$$ where $$|u|_{H^s({\Omega \cup \Omega_I})}^2 {\coloneqq}\int_{{\Omega \cup \Omega_I}} \int_{{\Omega \cup \Omega_I}}  \frac{(u(x)-u(y))^2}{||x-y||^{d+2s}}dydx.$$ Hence, there exist two positive constants $C_1$ and $C_2$ such that $$\begin{aligned}
C_1 ||u||_{H^s({\Omega \cup \Omega_I})} \leq |||u||| \leq C_2 ||u||_{H^s({\Omega \cup \Omega_I})} ~~~\forall~ u \in V_c({\Omega \cup \Omega_I}).\end{aligned}$$

This equivalence implies that $(V_c({\Omega \cup \Omega_I}),|||\cdot|||)$ is a Banach space for kernel functions of this class. Applying Lax-Milgram Theorem finally brings in the well posedness of problem (\[weakform\]).

Let $\Omega \subset {\mathbb{R}}^d$ be a bounded domain with piecewise smooth boundary and let the kernel function ${\gamma}$ satisfy the general requirement (\[gen\_assumption\]) and the specific assumption in (\[case\_1\]). Then for any linear functional $\ell$ there exists a unique $u \in V_c({\Omega \cup \Omega_I})$ such that $a(u,v) = \ell(v) $ for all $v$ in $V_c({\Omega \cup \Omega_I})$.

See [@wellposedness].

Finite-dimensional approximation
--------------------------------

With the concept of weak solutions we can proceed as in the local case to develop finite element approximations of (\[weakform\]).

Therefore, let ${\left\{}
\newcommand{\rk}{\right\}}V_c^N\rk_N$ be a sequence of finite-dimensional subspaces of $V_c({\Omega \cup \Omega_I})$, where $N = \dim (V_c^N)$, and let $u_N$ denote the solution of $$\begin{aligned}
\emph{\text{Find $u_N \in V_c^N$ such that } $a(u_N, v_N) = \ell(v_N)$ \text{ for all $v_N$ in $V_c^N$} }. \label{finiteelementproblem}\end{aligned}$$ Then from [@wellposedness] we recall the following result.

\[fem\_theorem\] Let $m$ be a nonnegative integer and $s \in (0,1)$. Further, let $ \Omega$ and ${\Omega \cup \Omega_I}$ be polyhedral domains and let $V_c^N$ consist of piecewise polynomials of degree no more than $m$. Assume that the triangulation is shape-regular and quasi-uniform as $h \to 0$ and suppose that ${u^* \in V_c({\Omega \cup \Omega_I}) \cap H^{m+t}({\Omega \cup \Omega_I})}$, where $t \in [s,1]$. Then there exists a constant $C>0$ such that, for sufficiently small $h$, $$\begin{aligned}
||u^*-u_N|| _{H^s({\Omega \cup \Omega_I})} \leq C h^{m+t-s}||u^*||_{H^{m+t}({\Omega \cup \Omega_I})}. \label{fem_estimate_fract}\end{aligned}$$

See [@wellposedness Theorem 6.2].

The derivation of the discretized problem then relies on the construction of the stiffness matrix. Therefore let ${\left\{}
\newcommand{\rk}{\right\}}\varphi_0, \ldots, \varphi_{N-1} \rk$ be a basis of $V_c^N$, such that the finite element solution $u^N \in V_c^N$ can be expressed as a linear combination $u^N = \sum_{k=0}^{N-1} u_k^N {{\varphi_k}}$. If the basis functions are chosen in a way such that ${{\varphi_k}}(x_k) = 1$ on appropriate grid points $x_k$, then the coefficients satisfy $u_k^N = u^N(x_k)$. We test for all basis functions, such that the finite element problem (\[finiteelementproblem\]) reads as $$\begin{aligned}
\textit{Find $u^N\in {\mathbb{R}}^{N}$ such that }\sum_{k=0}^{N-1} u_k^N a({{\varphi_k}}, {{\varphi_j}}) = \ell({{\varphi_j}}) =: b_j  \textit{ for } 0 \leq j < N.\end{aligned}$$ The stiffness matrix $A^N=(a_{kj})_{kj} \in {\mathbb{R}}^{N \times N}$ is given by $$\begin{aligned}
a_{kj} {\coloneqq}\int_{{\Omega \cup \Omega_I}} \int_{{\Omega \cup \Omega_I}}  ({{\varphi_k}}(y)-{{\varphi_k}}(x))({{\varphi_j}}(y)-{{\varphi_j}}(x)){\gamma}(x,y)dydx\end{aligned}$$ and we finally want to solve the discretized *Galerkin system* $$\begin{aligned}
A^Nu^N = b^N, \label{discretizedproblem}\end{aligned}$$ where $u^N, b^N \in {\mathbb{R}}^{N}$. The properties of the bilinear form $a$ imply that $A^N$ is symmetric and positive definite, such that there exists a unique solution $u^N$ of the finite-dimensional problem (\[discretizedproblem\]).

Finite element setting {#sec:femsetting}
======================

In this section we study a continuous Galerkin discretization of the homogeneous nonlocal Dirichlet problem, given by $$\begin{aligned}
2\int_{{\Omega \cup \Omega_I}} (u(x)-u(y)){\gamma}(x,y)dy = f(x) ~~~&(x\in \Omega) ,\nonumber\\
u(x) = 0    ~~~~~~~&(x\in \Omega_I),\end{aligned}$$ under the following assumptions:

-   We set $\Omega {\coloneqq}\prod_{i=0}^{d-1}[a_i,b_i]$, where $[a_i, b_i]$ are compact intervals on ${\mathbb{R}}$.

-   We assume that the kernel ${\gamma}$ is *translation invariant*, such that $${\gamma}(b+Qx, b+Qy) = {\gamma}(x,y)$$ for all $b\in{\mathbb{R}}^d$ and $Q \in {\mathbb{R}}^{d \times d}$ orthonormal.

Assumption (A1) allows for a simple triangulation of $\Omega$, which we use to define a finite-dimensional energy space $V^{N}_c$. Together with (A2) we can show that this discretization yields the multilevel Toeplitz structure of the stiffness matrix, where the order of the matrix is determined by the number of grid points in each respective space dimension.

Definition of the finite-dimensional energy space
-------------------------------------------------

We decompose the domain $\Omega = \prod_{i=0}^{d-1}[a_i,b_i]$ into $d$-dimensional hypercubes with sides of length $ h>0 $ in each respective dimension. Note that we can omit a discretization of $\Omega_I$ since we assume homogeneous Dirichlet volume constraints. Let $\operatorname{\mathbf{N}}= (N_i)_{0 \leq i < d} {\coloneqq}(\frac{b_i-a_i}{ h })_{0 \leq i < d}$ and ${\operatorname{\mathbf{L}}{\coloneqq}(N_i-1)_{0 \leq i < d}}$, then for the interior of $\Omega$ this procedure results in $\operatorname{\mathbf{L}}^d {\coloneqq}\prod_{i=0}^{d-1} L_i$ degrees of freedom. Due to the simple structure of the domain we can choose a canonical numeration for the resulting grid $\prod_{i=0}^{d-1} \left( a_i +  h {\left\{}
\newcommand{\rk}{\right\}}0, \ldots, L_i -1\rk \right) $ of inner points. More precisely, we will employ the map $$\begin{aligned}
E^{\operatorname{\mathbf{n}}}(z) {\coloneqq}\sum_{i=0}^{d-1} z_i p_i(\operatorname{\mathbf{n}}),\end{aligned}$$ where $p_i(\operatorname{\mathbf{n}}) {\coloneqq}\prod_{j>i}n_j$, for establishing an order on a structured grid $\prod_{i=0}^{d-1} {\left\{}
\newcommand{\rk}{\right\}}0, \ldots n_i-1 \rk$, where $\operatorname{\mathbf{n}}=(n_0,\ldots,n_{d-1}) \in \mathbb{N}^d$. Its inverse is given by $$\begin{aligned}
(E^{-\operatorname{\mathbf{n}}}(k))_{0 \leq i < d} =  (\lfloor \tfrac{k}{p_i(n)}\rfloor - \lfloor \tfrac{k}{p_{i-1(n)}}\rfloor n_i )_{0 \leq i < d}.\end{aligned}$$ Let $e {\coloneqq}(1,\ldots, 1)\in {\mathbb{R}}^{d}$ and $a {\coloneqq}(a_0,\ldots,a_{d-1})$, then we define the ordered array of inner grid points $
(x_k)_{0 \leq k < \operatorname{\mathbf{L}}^d} \in {\mathbb{R}}^{\operatorname{\mathbf{L}}^d \times d} 
$ by $$x_k   {\coloneqq}a  + h ( E^{-\operatorname{\mathbf{L}}}(k)+e)$$ for $0 \leq k < \operatorname{\mathbf{L}}^d$. We further define finite elements $S_k {\coloneqq}b_k +  h \square$, where $ (b_k^i)_{0 \leq i < d}   {\coloneqq}a +  h  E^{-\operatorname{\mathbf{N}}}(k) $ and $\square {\coloneqq}[0,1]^d$, such that $\Omega = \bigcup_{k=0}^{\operatorname{\mathbf{N}}^d - 1} S_k.$ Next we aim to define appropriate element basis functions on the reference element $\square$. Therefore we denote by $$(v_k)_{0 \leq k < 2^d} \in {\mathbb{R}}^{2^d \times d}$$ the vertices of the unit cube $\square$ ordered according to $
 v_k  {\coloneqq}E^{-(2,\ldots, 2)}(k).
$ Then for each vertex $v_k$, $0 \leq k < 2^d$, we define an *element basis function* $\psi_k \colon \square \to [0,1]$ by $$\begin{aligned}
\psi_k(x) &= \left(\prod_{i=0, v_k^i = 0}^{d-1} (1-x_i) \right)\left(\prod_{i=0, v_k^i = 1}^{d-1} x_i \right)   .\end{aligned}$$ For dimensions $d \in {\left\{}
\newcommand{\rk}{\right\}}1,2,3 \rk$ respectively, these are the usual linear, bilinear and trilinear element basis functions (see e.g. [@fem_book Chapter 1]). They are defined in a way such that $0 \leq \psi_k \leq 1$ and $\psi_k(v_k) = 1$. Moreover, we define the *reference basis function* $\varphi \colon {\mathbb{R}}^d \to [0,1]$ by $$\begin{aligned}
\varphi(x) &{\coloneqq}\begin{cases}
\psi_i(v_i + x)&: x \in (\square - v_i)\\
0&:else.
\end{cases} \label{reference_basisfunction}\end{aligned}$$ We note that $J{\coloneqq}[-1,1]^d = \dot\bigcup_{i=0}^{2^d} (\square -v_i ) $ (disjoint union), such that $\varphi$ is well defined and $supp(\varphi) =  J$. Now let the physical support be defined as $$I_k {\coloneqq}\bigcup {\left\{}
\newcommand{\rk}{\right\}}S_i \colon  x_k \in S_i, 0 \leq i < \operatorname{\mathbf{N}}^d \rk ,$$ which is a patch of the elements touching the node $x_k$. We associate to each element $S_k$ the transformation $T_k \colon  J \to I_k$, $T_k(v) {\coloneqq}x_k +  h  v $. We note that $\det dT_k(x) \equiv h^d  $. Then for each node $x_k$ we define a basis function $\varphi_k \colon {\Omega \cup \Omega_I}\to [0,1]$ by $$\begin{aligned}
\varphi_k(x) &{\coloneqq}\begin{cases}
\varphi(T_k^{-1}(x))&: x \in I_k\\
0&:else
\end{cases}\\
&=\begin{cases}
\psi_i(v_i + T_k^{-1}(x))&: T_k^{-1}(x)\in (\square - v_i)\\
0&:else,
\end{cases}\end{aligned}$$ which satisfies $0\leq \varphi_k \leq 1$ and $\varphi_k(x_k) = 1$. Figure \[fig:basisfunction\] illustrates the latter considerations for $d=2$, $a = (0,0)$ and $b=(1,1)$.

Finally, we can define a constrained finite element space by $$\begin{aligned}
V_c^{ h } {\coloneqq}\operatorname{span}{\left\{}
\newcommand{\rk}{\right\}}\varphi_k \colon 0 \leq k < \operatorname{\mathbf{L}}^d \rk \label{fem_space},\end{aligned}$$ such that each linear combination consisting of a set of these basis functions fulfills the homogeneous Dirichlet volume constraints. Notice that we parametrize these spaces by the grid size $ h $ indicating the dimension $\operatorname{\mathbf{L}}^d$, which is by definition a function of $ h $. Finally, we close this subsection with the following observations, which we exploit above.

Let $x \in {\mathbb{R}}^d$, then:

-   $\varphi(x) = \varphi(|x|)$, where $|x| {\coloneqq}( |x_i|)_i$.

-   Let $R_i \colon {\mathbb{R}}^d \to {\mathbb{R}}^d$, for $0 \leq i < d$, denote the reflection $
            R_i(x) = (x_0, \ldots, -x_i, \ldots, x_{d-1}),
            $ then $$\begin{aligned}
            \varphi(x)  = \varphi (|x| ) \Leftrightarrow \varphi(x)  = \varphi (R_i(x)) ~~\forall~ 0 \leq i < d. \label{reflection_invariance}
            \end{aligned}$$

-   $\varphi(x)    = \varphi((x_{\sigma(i)})_i)$ for all permutations $\sigma \colon {\left\{}
    \newcommand{\rk}{\right\}}0,\ldots, d-1 \rk \to {\left\{}
    \newcommand{\rk}{\right\}}0,\ldots, d-1 \rk$.

We first show i). Since $\varphi(x) = 0 = \varphi(|x|)$ for $x$ in $int({J})^c$, let $x \in int({J})$. Thus, there exists an index $0 \leq k < d $ such that $x \in int(\square) - v_k,$ which implies that $x_i < 0$ if and only if $v_k^i = 1. $ Hence, we can conclude that $$\begin{aligned}
    \varphi(x) &= \psi_k(x + v_k) \\
    &= \left(\prod_{i=0, v_k^i = 0}^{d-1} (1-x_i) \right)\left(\prod_{j=0, v_k^i = 1}^{d-1} (1+x_i) \right) \\
    &= \left(\prod_{i=0, v_k^i = 0}^{d-1} (1-|x_i|) \right)\left(\prod_{j=0, v_k^i = 1}^{d-1} (1-|x_i|) \right) \\
    &= \prod_{i=0}^{d-1} (1-|x_i|) \\
    &=\psi_0 (|x| + v_0) \\
    &= \varphi(|x|).
    \end{aligned}$$ Then, on the one hand, we have that $|x| = |R_i(x)|$ and therefore $\varphi(x)  = \varphi (R_i(x))$ for all $0 \leq i < d$. On the other hand, we note that the operation $|\cdot|$ is a composition of reflections $R_i$, more precisely $$|x| = \left(\prod_{x_i < 0}R_i\right)(x) .$$ Thus, we obtain the equivalence stated in ii). Statement iii) follows from the representation $
\varphi(x)  = \varphi(|x|) 
    =\prod_{i=0}^{d-1} (1-|x_i|)
$ due to the commutativity of the product.

Multilevel Toeplitz structure of the stiffness matrix
-----------------------------------------------------

Now we aim to show that the stiffness matrix $A$ owns the structure of a $d$-level Toeplitz matrix. This is decisive for this work, since it finally enables us to solve the discretized system (\[discretizedproblem\]) in an affordable way.

From now on the assumption (A2) on the kernel function becomes relevant. At this point we note that the mapping $(x,y) \mapsto \mathcal{X}_{B_R(x)}(y) $ is translation invariant, since translations of the upper form $x \mapsto b + Qx$ are isometric in Euclidean space. This also holds if we define the ball $B_R$ with respect to another than the $||\cdot||_2$-norm. Hence, we can regard the kernel as a function $${\gamma}(x,y) = \mathcal{X}_{B_R(x)}(y)\sigma(x,y),~~x,y \in {\mathbb{R}}^d ,$$ for some translation invariant function $\sigma$. In order to analyze the multilevel structure of ${A \in {\mathbb{R}}^{\operatorname{\mathbf{L}}^d \times \operatorname{\mathbf{L}}^d}}$ it is convenient to introduce an appropriate multi-index notation. To this end, we choose $E^{\operatorname{\mathbf{L}}}$ from above as index bijection and we identify $a_{\operatorname{\mathbf{i}}\operatorname{\mathbf{j}}} = a_{E^{\operatorname{\mathbf{L}}}(\operatorname{\mathbf{i}}),E^{\operatorname{\mathbf{L}}}(\operatorname{\mathbf{j}})}. $ We call the matrix $A$ *$d$-level Toeplitz* if $$a_{\operatorname{\mathbf{i}}\operatorname{\mathbf{j}}} = a(\operatorname{\mathbf{i}}-\operatorname{\mathbf{j}}).$$ If even $a_{\operatorname{\mathbf{i}}\operatorname{\mathbf{j}}} = a(|\operatorname{\mathbf{i}}-\operatorname{\mathbf{j}}|)$, then each level is symmetric and we can reconstruct the whole matrix from the first row (or column). For a more general and detailed consideration of multilevel Toeplitz matrices see for example [@multileveltoeplitz]. However, with this notation at hand we can now formulate

\[multilevel\_toeplitz\] Let $\Omega = \prod_{i=0}^{d-1}[a_i,b_i]$ and assume that the kernel ${\gamma}$ is translation invariant and satisfies the general assumptions (\[gen\_assumption\]) as well as the specific assumption (\[case\_1\]). Further let the finite element space $V_c^{ h }$ be defined as in (\[fem\_space\]) for a grid size $ h >0$. Then the stiffness matrix $A$ associated to problem (\[finiteelementproblem\]) is $d$-level Toeplitz, where each level is symmetric.

The key point in the proof is the relation $$\begin{aligned}
a_{kj} = a(|{ h }^{-1}(x_k-x_j)|),\end{aligned}$$ which we show in two steps. First we show that $a_{kj} = a({ h }^{-1}(x_k-x_j)) $ and then we proof $a(z) = a(|z|).$ Therefore let us recall that the entry $a_{kj}$ of the stiffness matrix $A$ is given by $$a_{kj} = \int_{{\Omega \cup \Omega_I}} \int_{{\Omega \cup \Omega_I}}  ({{\varphi_k}}(y)-{{\varphi_k}}(x))({{\varphi_j}}(y)-{{\varphi_j}}(x)){\gamma}(x,y)dydx.$$ Having a closer look at the support of the integrand, we find that $$\begin{aligned}
({{\varphi_k}}(y)-{{\varphi_k}}(x))({{\varphi_j}}(y)-{{\varphi_j}}(x)) = 0  \Leftrightarrow (x,y) \in (I_k^c \times I_k^c) \cup (I_j^c \times I_j^c) \cup {\left\{}
\newcommand{\rk}{\right\}}(x,x)\colon x \in {\Omega \cup \Omega_I}\rk.\end{aligned}$$ Since ${\left\{}
\newcommand{\rk}{\right\}}(x,x)\colon x \in {\Omega \cup \Omega_I}\rk$ has null $\lambda_{2d}$-Lebesgue measure we can neglect it in the integral and obtain $$\begin{aligned}
a_{kj} = \int_{\left( I_k^c \times I_k^c\right)^c \cap \left(I_j^c \times I_j^c \right)^c} ({{\varphi_k}}(y)-{{\varphi_k}}(x))({{\varphi_j}}(y)-{{\varphi_j}}(x)){\gamma}(x,y)dydx.\end{aligned}$$\
Aiming to show $a_{kj} = a({ h }^{-1}(x_k-x_j))$ we need to carry out some basic transformations of this integral. Since by definition $\varphi_j = \varphi \circ T_j^{-1}$ and also $\det T_j(x) \equiv h^d$ we find $$\begin{aligned}
a_{kj} &= \int_{\left( I_k^c \times I_k^c\right)^c \cap \left(I_j^c \times I_j^c \right)^c} 
({{\varphi_k}}(y)-{{\varphi_k}}(x))({{\varphi_j}}(y)-{{\varphi_j}}(x)){\gamma}(x,y)dydx\\
&=  h ^{2d}\int_{T_j^{-1}({\mathbb{R}}^d) \times T_j^{-1}({\mathbb{R}}^d) }
(1-\mathcal{X}_{I_j^c \times I_j^c }(T_j(v),T_j(w)))
(1-\mathcal{X}_{I_k^c \times I_k^c }(T_j(v),T_j(w)))\\
&~~~~~~~~~~~~ ((\varphi \circ  T_k^{-1})(T_j(w))-(\varphi \circ T_k^{-1})(T_j(v)))(\varphi(w)-\varphi(v)){\gamma}(T_j(v),T_j(w))dwdv.\end{aligned}$$ Now we make a collection of observations. Due to assumption (A2) we have $${\gamma}(T_j(v),T_j(w)) = {\gamma}(x_j +  h  v,x_j +  h  w) = {\gamma}( h  v,  h  w).$$ Furthermore, by definition of the transformations $T_j, T_k$ we find that $T_j^{-1}({\mathbb{R}}^d)={\mathbb{R}}^d$ as well as $(T_k^{-1}\circ T_j)(v) =   h ^{-1}(x_j +  h  v - x_k) =  h ^{-1}(x_j - x_k)+ v $. Since these transformations are bijective we also have that $\mathcal{X}_{M^c\times M^c }(T_j(x),T_j(y)) = \mathcal{X}_{(T_j^{-1}(M))^c\times (T_j^{-1}(M))^c }(x,y)$ for a set $M\subset {\mathbb{R}}^d$. Hence, defining $x_{jk}  {\coloneqq}h ^{-1}(x_j - x_k) = - x_{kj} $ and recognizing $T_j^{-1}(I_k) = x_{kj} + J$ we finally obtain $$\begin{aligned}
a_{kj}
&= { h }^{2d}\int_{\left( J^c \times J^c\right)^c \cap \left((x_{kj}+J )^c \times (x_{kj}+J )^c\right)^c}
\\
&~~~~~~~~~~~~ (\varphi(w-x_{kj})-\varphi(v-x_{kj}))(\varphi(w)-\varphi(v)){\gamma}( h  w, h  v)dwdv\\
&=a(x_{kj}).\end{aligned}$$

Next, we proof that this functional relation fulfills $a(z) = a(|z|)$. Let us for this purpose define $F(x,y;z) {\coloneqq}(\varphi(y-z)-\varphi(x-z))(\varphi(y)-\varphi(x)){\gamma}( h  y, h  x)$ such that $$a(z) = { h }^{2d} \int_{\left( J^c \times J^c\right)^c \cap \left((z+J )^c \times (z+J )^c \right)^c}
F(x,y;z)dydx.$$ Let $z \in {\left\{}
\newcommand{\rk}{\right\}}x_{kj}: 0 \leq k,j<\operatorname{\mathbf{L}}^d \rk$. Then there exists an orthonormal matrix $R= R(z)  \in {\mathbb{R}}^{d \times d}$, which is a composition of reflections $R_i$ from (\[reflection\_invariance\]), such that $Rz = |z|$. Then from (\[reflection\_invariance\]) and the assumption (A2) on the kernel, we obtain for $x,y \in {\mathbb{R}}^d$ that $$\begin{aligned}
F(Rx,Ry;|z|) & = (\varphi(Ry-Rz)-\varphi(Rx-Rz))(\varphi(Ry)-\varphi(Rx)){\gamma}( h  Ry, h  Rx) \\
&= (\varphi(y-z)-\varphi(x-z))(\varphi(y)-\varphi(x)){\gamma}( h  y, h  x) \\
&=F(x,y;z).\end{aligned}$$ Since $R(J ) = J  $ and therefore $$\begin{aligned}
&R\left(\left( J^c \times J^c\right)^c \cap \left((z+J )^c \times (z+J )^c\right)^c\right) \\
=&\left( J^c \times J^c\right)^c \cap \left((|z|+J )^c \times (|z|+J )^c\right)^c,\end{aligned}$$ we eventually obtain $$\begin{aligned}
a(|z|)&={ h }^{2d} \int_{\left( J^c \times J^c\right)^c \cap \left((|z|+J )^c \times (|z|+J )^c\right)^c}
F(x,y;|z|)dydx \\
&={ h }^{2d} \int_{\left( J^c \times J^c\right)^c \cap \left((z+J )^c \times (z+J )^c\right)^c}
F(Rx,Ry;|z|)dydx \\
&={ h }^{2d} \int_{\left( J^c \times J^c\right)^c \cap \left((z+J )^c \times (z+J )^c\right)^c}
F(x,y;z)dydx \\
&=a(z).\end{aligned}$$ Finally, we can show that $A$ carries the structure of a $d$-level Toeplitz matrix. By having a closer look at the definitions of $E^{\operatorname{\mathbf{L}}}$ and the grid points $x_k$ we can conclude that $$\begin{aligned}
a_{\operatorname{\mathbf{i}}\operatorname{\mathbf{j}}} &= a_{E^{\operatorname{\mathbf{L}}}(\operatorname{\mathbf{i}})E^L(\operatorname{\mathbf{j}})}
  = a(|{ h }^{-1}(x_{E^{\operatorname{\mathbf{L}}}(\operatorname{\mathbf{i}})}-x_{E^{\operatorname{\mathbf{L}}}(\operatorname{\mathbf{j}})})|)
= a(|\operatorname{\mathbf{i}}- \operatorname{\mathbf{j}}|).\end{aligned}$$ Thus, the entry $a_{\operatorname{\mathbf{i}}\operatorname{\mathbf{j}}}$ only depends on the difference $\operatorname{\mathbf{i}}- \operatorname{\mathbf{j}}$.

As a consequence, in order to implement the matrix-vector product, it is sufficient to assemble solely the first row or column $$M {\coloneqq}(a_{\ell0})_\ell = (a({ h }^{-1}(x_\ell-x_0)))_\ell = (a(E^{\operatorname{\mathbf{L}}}(\ell)))_\ell$$ of the stiffness matrix $A$, since for $\ell(k,j) {\coloneqq}E^L( h ^{-1}(|x_k-x_j|))$ we get $$a_{kj} = a( h ^{-1}(|x_k-x_j|)) = a(E^{\operatorname{\mathbf{L}}}(\ell(k,j))) = M_{\ell(k,j)}.$$ Note that $\ell(k,j) {\coloneqq}E^{\operatorname{\mathbf{L}}}( h ^{-1}(|x_k-x_j|))$ is well defined, since $ h ^{-1}(|x_k-x_j|) $ lies in the domain of definition of $E^{\operatorname{\mathbf{L}}}$.

\[L1\_radial\] Exploiting that $\varphi$ is invariant under permutations, the same proof (by composing the reflection $R$ with a permutation matrix) shows that $a(z)   = a((z_{\sigma(i)})_i)$ for all permutations $\sigma \colon {\left\{}
\newcommand{\rk}{\right\}}0,\ldots, d-1 \rk \to {\left\{}
\newcommand{\rk}{\right\}}0,\ldots, d-1 \rk$. We will use this observation to accelerate the assembling process.

Assembling procedure {#sec:assembling}
====================

In this section we aim to analyze the entries $a_{kj} = a(x_{kj})$ of the stiffness matrix $A$ more closely and derive a representation which can be efficiently implemented.

We first characterize the domain of integration occurring in the integral in $a(x_{kj})$. Let us define ${J_{kj} {\coloneqq}(x_{kj}+J)} $ then $$\begin{aligned}
&\left( J^c \times J^c \cup J_{kj}^c \times J_{kj}^c \right)^c =\left( J^c \times J^c\right)^c \cap \left(J_{kj}^c \times J_{kj}^c \right)^c\\
=&\left( (J \times J) \cup (J^c \times J) \cup (J \times J^c) \right) \cap \left((J_{kj} \times J_{kj}) \cup (J_{kj}^c \times J_{kj}) \cup (J_{kj} \times J_{kj}^c)\right) \\
=&(C \times C) \cup (D_k \times C) \cup (C \times D_k)
\cup (D_j \times C) \cup (J^c \cap J_{kj}^c \times C) \cup (D_j \times D_k)\\
&\cup (C \times D_j) \cup (D_k \times J_{kj}) \cup (C \times J^c \cap J_{kj}^c),
$$ where we set $C {\coloneqq}J \cap J_{kj}$, $D_k {\coloneqq}J  \cap J_{kj}^c  $ and $D_j {\coloneqq}J^c \cap J_{kj} $. By exploiting the symmetry of the integrand we thus get $$\begin{aligned}
a_{kj}/ h ^{2d} = & \int_{J \cap J_{kj} } \int_{(J \cap J_{kj})}  F(x,y;x_{kj})dydx\\
&+2\int_{J \cap J_{kj}} \int_{J_{kj}^c \cap J }  F(x,y;x_{kj})dydx \\
&+2\int_{J \cap J_{kj}} \int_{J_{kj} \cap J^c }  F(x,y;x_{kj})dydx \\
&+2\int_{J \cap J_{kj}} \int_{J^c \cap J_{kj}^c}  F(x,y;x_{kj})dydx \\
&+2\int_{J \cap J_{kj}^c} \int_{J^c \cap J_{kj} }  F(x,y;x_{kj})dydx ,\end{aligned}$$ where $F(x,y;z) =  (\varphi(y-z)-\varphi(x-z))(\varphi(y)-\varphi(x)){\gamma}( h  y, h  x)$. Note that this representation holds for a general setting without assuming (A1) and (A2). For implementation purpose we additionally require from now on:

-   We assume that $R \geq diam(\Omega) = ||\operatorname{b}- a ||_2$, such that $\Omega \subset B_R(x)$ for all $x \in \Omega$.

This third assumption simplifies the domain of integration in the occurring integrals in the sense that we can omit the intersection with the ball $B_{R}( h  x)$ which is hidden behind the kernel function. This coincides with the application to space-fractional diffusion problems where we aim to model $R\to \infty$ (see Section \[sec:numres\]). Furthermore, since we can construct the whole stiffness matrix $A$ from the first row $M$, it is convenient to introduce the following $\operatorname{\mathbf{L}}^d$-dimensional vectors: $$\begin{aligned}
{\texttt{sing}}_k &{\coloneqq}{ h }^{2d} \int_{J \cap J_k  } \int_{J \cap J_k}  F(x,y;E^{\operatorname{\mathbf{L}}}(k))dydx \\
&+2{ h }^{2d}\int_{J \cap J_k } \int_{J_k^c \cap J }  F(x,y;E^{\operatorname{\mathbf{L}}}(k))dydx \\
&+2{ h }^{2d}\int_{J \cap J_k } \int_{J_k \cap J^c }  F(x,y;E^{\operatorname{\mathbf{L}}}(k))dydx, \\
{\texttt{rad}}_k &{\coloneqq}2{ h }^{2d}\int_{J \cap J_k } \int_{  (J \cup J_k )^c}  F(x,y;E^{\operatorname{\mathbf{L}}}(k))\mathcal{X}_{B_R( h  x)}( h  y) dydx ,\\
{\texttt{dis}}_k &{\coloneqq}2{ h }^{2d}\int_{J \cap J_{k}^c} \int_{J^c \cap J_k }  F(x,y;E^{\operatorname{\mathbf{L}}}(k))dydx ,\end{aligned}$$ for $0 \leq k < \operatorname{\mathbf{L}}^d$, where $J_k {\coloneqq}J_{0k} = E^{\operatorname{\mathbf{L}}}(k)+J$. Since $x_{k0} = \tfrac{x_k-x_0}{ h } = E^{\operatorname{\mathbf{L}}}(k)$, we have ${M = {\texttt{sing}}+{\texttt{rad}}+ {\texttt{dis}}}$. As we will see in the subsequent program, each of these vectors requires a different numerical handling which justifies this separation. In these premises we point out, that on the one hand we may touch possible singularities of the kernel function along the integration in ${\texttt{sing}}_k$. On the other hand, the computation of ${\texttt{rad}}_k$ may require the integration over a “large” domain if $R \to \infty$ (e.g. fractional kernel). Both are numerically demanding tasks and complicate the assembling process. In contrast to that, the computation of ${\texttt{dis}}_k$ turns out to be numerically viable without requiring a special treatment. However, we fortunately find for $k$ with $J \cap J_k =  \emptyset$ that ${\texttt{sing}}_k=0 ={\texttt{rad}}_k$. Hence, it is worth identifying those indices and treat them differently in the assembling loop. Therefore, from $J = \dot\bigcup_{i=0}^{2d} (\square -v_i )$ we can deduce that $J \cap J_k \neq \emptyset$ if and only if $k  = E^{\operatorname{\mathbf{L}}}(v_i)$ for an index $0 \leq i < 2^d$. As a consequence, we only have to compute ${\texttt{sing}}_k$ and ${\texttt{rad}}_k$ for $k \in {\texttt{idx}}_0 {\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}E^{\operatorname{\mathbf{L}}}(v_i) \colon 0 \leq i < 2^d  \rk  $. We can cluster these indices even more. For that reason let us define on ${\texttt{idx}}_0$ the equivalence relation $$k \sim j :\Leftrightarrow \exists~ \text{permutation matrix}~ P \in {\mathbb{R}}^{\operatorname{\mathbf{L}}^d \times \operatorname{\mathbf{L}}^d}\colon ~E^{-\operatorname{\mathbf{L}}}(k)  =P E^{-\operatorname{\mathbf{L}}}(j).$$ Then due to Remark \[L1\_radial\] after Theorem \[multilevel\_toeplitz\] we have to compute the values ${\texttt{sing}}_k$ and ${\texttt{rad}}_k$ only for $k \in  {\texttt{idx}}_s^k{\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}[j]_\sim \colon j \in {\texttt{idx}}_0 \rk$. In order to make this more precise, we figure out that the quotient set can further be specified as $${\texttt{idx}}_s^k = {\left\{}
\newcommand{\rk}{\right\}}[E^{\operatorname{\mathbf{L}}}(z)]_\sim  \colon z \in S \rk,$$ where $$S {\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}(0, \ldots, 0), (1,0, \ldots, 0),(1,1,0, \ldots, 0), \ldots,  (1,1, \ldots, 1)\rk \subset {\left\{}
\newcommand{\rk}{\right\}}v_i \colon 0 \leq i < 2^d\rk$$ with $|S| = d+1 \leq 2^d$. With other words, we group those $v_i$ which are permutations of one another. The associated indices $0 \leq i < 2^d$ are thus given by ${\left\{}
\newcommand{\rk}{\right\}}E^{(2,\ldots,2)}(z) \colon z \in S \rk =:{\texttt{idx}}_s^i$.

In addition to the preceding considerations, we also want to partition the integration domains $J \cap J_k$, $J \cap J_k^c$ and $J^c \cap J_k$, where $k\in{\texttt{idx}}_s^k$, into cubes $(\square - v_\nu)$, such that we can express the reference basis function $\varphi$ with the help of the element basis functions $\psi_i$. This is necessary in order to compute the integrals in a vectorized fashion and obtain an efficient implementation. Let us start with $J \cap J_k$, where $k = E^{\operatorname{\mathbf{L}}}(v_i)$ for $0 \leq i < 2^d $ such that $J_k = v_i + J$. Then we define the set $$\begin{aligned}
D_i &{\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}0 \leq \mu < 2^d \colon  \square - v_\mu \in J \cap J_k \rk 
= {\left\{}
\newcommand{\rk}{\right\}}0 \leq \mu < 2^d \colon  \exists \kappa \colon v_\mu + v_i = v_\kappa \rk\\
&= {\left\{}
\newcommand{\rk}{\right\}}0 \leq \mu < 2^d \colon (v_\mu + v_i)_j < 2 ~\forall~ 0 \leq j < d \rk .\end{aligned}$$ Since $J = \dot\bigcup_{i=0}^{2^d} (\square -v_i )$, we find that $J \cap J_k = \bigcup_{\nu\in D_i}( \square - v_\nu)$. With this, we readily recognize that $J \cap J_k^c = \bigcup_{\nu\in{\left\{}
\newcommand{\rk}{\right\}}0, \ldots, 2^d - 1 \rk  \backslash D_i} (\square - v_\nu)$. Similarly, one can derive a set $D_i^c$, such that ${J^c \cap J_k = v_i + \bigcup_{\nu\in D_i^c}( \square - v_\nu)}$. Figure \[fig:indexsets\] illustrates the latter considerations and in Table \[indexsets\] these sets are listed for dimensions $d\in {\left\{}
\newcommand{\rk}{\right\}}1,2,3\rk$, respectively.

 \

  -------------------------------------------------------------------------------------------------------------------------------------------------------------------
   $d$   $i \in {\texttt{idx}}_s^i$   $k = E^L(v_i)\in {\texttt{idx}}_s^k$          $D_i$                      $D_i^c = {\left\{}                 $\kappa(D_i, i)$
                                                                                                   \newcommand{\rk}{\right\}}0\leq j < i \rk$  
  ----- ---------------------------- -------------------------------------- --------------------- -------------------------------------------- ----------------------
    1               $0$                                0                           $(0,1)$                        $\emptyset$                         $ (0,1)$

                    $1$                                1                             $0$                              $0$                               $ 1$

    2               $0$                               $0$                        $(0,1,2,3)$                      $\emptyset$                       $ (0,1,2,3)$

                    $2$                              $p_0$                         $(0,1)$                          $(0,1)$                           $ (2,3)$

                    $3$                           $p_0 + p_1$                        $0$                           $(0,1,2)$                            $ 3$

    3               $0$                               $0 $                   $(0,1,2,3,4,5,6,7)$                  $\emptyset$                   $ (0,1,2,3,4,5,6,7)$

                    $4$                              $p_0$                       $(0,1,2,3)$                      $(0,1,2,3)$                       $ (4,5,6,7)$

                    $6$                           $p_0 + p_1 $                     $(0,1)$                      $(0,1,2,3,4,5)$                       $ (6,7)$

                    $7$                         $p_0+p_1+p_2  $                      $0$                       $(0,1,2,3,4,5,6)$                        $ 7$
  -------------------------------------------------------------------------------------------------------------------------------------------------------------------

  : Index sets for the implementation.[]{data-label="indexsets"}

With this at hand, we can now have a closer look at the vectors ${\texttt{sing}}$, ${\texttt{rad}}$ and ${\texttt{dis}}$ and put them into a form which is suitable for the implementation.

### sing {#sing .unnumbered}

Let $i \in {\texttt{idx}}_s^i$, such that $E^{\operatorname{\mathbf{L}}}(k) = v_i$. Since $J \cap J_k = \bigcup_{\nu\in D_i} \square - v_\nu$, we can transform the first integral in ${\texttt{sing}}$ as follows $$\begin{aligned}
  \int_{J \cap J_k } \int_{J \cap J_k}  F(x,y;E^{\operatorname{\mathbf{L}}}(k))dydx
&= \sum_{\nu \in D_i} \sum_{\mu \in D_i} \int_{\square - v_\nu} \int_{\square - v_\mu}  F(x,y; v_i)dydx \\
&= \sum_{\nu \in D_i} \sum_{\mu \in D_i} \int_{\square } \int_{\square }  F(x- v_\nu,y- v_\mu; v_i)dydx .\end{aligned}$$ By definition of $D_i$ we have for $\mu \in D_i$ that $v_\mu + v_i = v_{\kappa(\mu, i)}$ with $\kappa(\mu,i) = E^{(2,\ldots,2)}(v_i+v_\mu)$ and since $$F(x- v_\nu,y- v_\mu; v_i) =  (\varphi(y-v_{\kappa(\mu, i)})-\varphi(x-v_{\kappa(\nu, i)}))(\varphi(y-v_\mu)-\varphi(x-v_\nu)){\gamma}( h (y-v_\mu), h (x-v_\nu))$$ we find due to (\[reference\_basisfunction\]) that $$\begin{aligned}
        & \sum_{\nu \in D_i} \sum_{\mu \in D_i} \int_{\square } \int_{\square }  F(x- v_\nu,y- v_\mu; v_i)dydx \\
        =& \sum_{\nu \in D_i} \sum_{\mu \in D_i} \int_{\square } \int_{\square }  ((\psi_{\kappa(\mu, i)}(y)-\psi_{\kappa(\nu, i)}(x))(\psi_{\mu}(y)-\psi_{\nu}(x))){\gamma}( h (y- v_\mu), h (x- v_\nu))dydx.     \end{aligned}$$ We separate the case $\mu = \nu$ since the kernel may have singularities at $(x,y)$ with $x=y$ and therefore these integrals need a different numerical treatment. At this point we find that $$\begin{aligned}
&\sum_{\nu \in D_i}\int_{\square } \int_{\square }  ((\psi_{\kappa(\nu, i)}(y)-\psi_{\kappa(\nu, i)}(x))(\psi_{\nu}(y)-\psi_{\nu}(x))){\gamma}( h  y, h  x)dydx\\
=&|D_i|\int_{\square } \int_{\square }  ((\psi_{i}(y)-\psi_{i}(x))(\psi_{0}(y)-\psi_{0}(x))){\gamma}( h  y, h  x)dydx .\end{aligned}$$ This simplification follows from some straightforward transformations exploiting assumption (A2) and the fact that each element basis function $\psi_\nu$ can be expressed by $\psi_0$ through the relation $\psi_\nu \circ g_\nu = \psi_0 ,$ where $g_\nu(x) {\coloneqq}v_\nu + R_\nu x$ for an appropriate rotation matrix $R_\nu$. With this observation we also find that the other two integrals in ${\texttt{sing}}$ are equal, i.e., $$\begin{aligned}
\int_{J \cap J_k} \int_{J_k^c \cap J }  F(x,y;E^{\operatorname{\mathbf{L}}}(k))dydx 
=\int_{J \cap J_k} \int_{J_k \cap J^c }  F(x,y;E^{\operatorname{\mathbf{L}}}(k))dydx.\end{aligned}$$ By exploiting again that $J \cap J_k = \bigcup_{\nu\in D_i} \square - v_\nu$ and thus $J \cap J_k^c = \bigcup_{\nu\in{\left\{}
\newcommand{\rk}{\right\}}0, \ldots, 2^d - 1 \rk  \backslash D_i} \square - v_\nu$ we obtain $$\begin{aligned}
&\int_{J \cap J_k} \int_{J_k^c \cap J }  F(x,y;E^{\operatorname{\mathbf{L}}}(k))dydx \\
=&-\sum_{ \nu \in D_i} \sum_{\mu \in   {\left\{}
\newcommand{\rk}{\right\}}0, \ldots,  2^d \rk  \backslash D_i}\int_{\square } \int_{\square } \psi_{\kappa(\nu, i)}(x) (\psi_{\mu}(y)-\psi_{\nu}(x)){\gamma}( h (v_\nu- v_\mu), h (x-y ))dydx.\end{aligned}$$ Note that by definition of $D_i$ we have that for $\mu \in   {\left\{}
\newcommand{\rk}{\right\}}0, \ldots, 2^d \rk  \backslash D_i$ there is no $\kappa \in {\left\{}
\newcommand{\rk}{\right\}}0, \ldots,  2^d \rk$ such that $v_\mu + v_i = v_\kappa$ and therefore $\varphi(y - (v_\mu + v_i)) = 0$ for all $y \in \square$. All in all we have

$$\begin{aligned}
&{\texttt{sing}}_k/ h ^{2d} =\\
& |D_i|\int_{\square } \int_{\square }  ((\psi_{i}(y)-\psi_{i}(x))(\psi_{0}(y)-\psi_{0}(x))){\gamma}( h  y, h  x)dydx  \\
&+ \sum_{\nu \in D_i}[ \sum_{\mu \in D_i \mu \neq \nu} \int_{\square } \int_{\square }  ((\psi_{\kappa(\mu, i)}(y)-\psi_{\kappa(\nu, i)}(x))(\psi_{\mu}(y)-\psi_{\nu}(x))){\gamma}( h (v_\mu- v_\nu), h (y- x))dydx \\     
&-4 \sum_{\mu \in   {\left\{}
\newcommand{\rk}{\right\}}0, \ldots,  2^d\rk  \backslash D_i}\int_{\square } \int_{\square } \psi_{\kappa(\nu, i)}(x) (\psi_{\mu}(y)-\psi_{\nu}(x)){\gamma}( h (v_\mu- v_\nu), h (y-x ))dydx ].\end{aligned}$$

### rad {#rad .unnumbered}

Let $i \in {\texttt{idx}}_s^i$, such that $E^{\operatorname{\mathbf{L}}}(k) = v_i$. Proceeding as above we obtain $$\begin{aligned}
{\texttt{rad}}_k &= 2 h ^{2d}\int_{J \cap J_k} \int_{ (J \cup J_k)^c}  F(x,y;v_i)dydx \\
&= 2 h ^{2d}\sum_{\nu \in D_i}\int_{\square} \psi_\nu(x)\psi_{\kappa(\nu,i)}(x)  \int_{   (J \cup J_k)^c} {\gamma}( h (x-v_\nu),  h  y) \mathcal{X}_{B_{R}( h (x-v_\nu))}( h  y)dydx.\end{aligned}$$ Let us define $$\begin{aligned}
P_\nu(x) {\coloneqq}&  \int_{   (J \cup J_k)^c} {\gamma}( h (x-v_\nu),  h  y) \mathcal{X}_{B_{R}( h (x-v_\nu))}( h  y)dy.\end{aligned}$$ Again, we can use $\psi_\nu \circ g_\nu = \psi_0$ in order to show by some straightforward transformations that $$\int_{\square} \psi_0(x)\psi_{i}(x) P_0(x)dx = \int_{\square}\psi_\nu(x)\psi_{\kappa(\nu,i)}(x)  P_\nu(x)dx$$ for all $\nu \in D_i$. Hence, we get

$${\texttt{rad}}_k = 2  h ^{2d}|D_i|\int_{\square} \psi_0(x)\psi_{i}(x) P_0(x)dx .$$

### dis {#dis .unnumbered}

Now we distinguish between the case where ${\texttt{rad}}_k$ and ${\texttt{sing}}_k$ are zero and the complement case. First let $i \in {\texttt{idx}}_s^i$, such that $E^{-\operatorname{\mathbf{L}}}(k) = v_i$, then we find $$\begin{aligned}
{\texttt{dis}}_k &{\coloneqq}2{ h }^{2d}\int_{J_k^c \cap J} \int_{J_k \cap J^c }  F(x,y;E^{-\operatorname{\mathbf{L}}}(k))dydx \\
& = -2{ h }^{2d} \sum_{\nu \in {\left\{}
\newcommand{\rk}{\right\}}0,\ldots,  2^d\rk  \backslash D_i} \sum_{\mu \in  D_i^c} \int_\square \int_\square   \psi_\mu(y)\psi_\nu(x) {\gamma}( h (v_\mu - v_\nu - v_i), h (y-x)) dydx.\end{aligned}$$ Now let $k \neq E^{-\operatorname{\mathbf{L}}}(v_i)$ for any $0 \leq i < 2^d$, then $J_k^c \cap J = J$ and $J_k \cap J^c =J_k$ such that $$\begin{aligned}
\frac{{\texttt{dis}}_k}{2{ h }^{2d}}= \int_{J} \int_{J_k}  F(x,y;E^{-\operatorname{\mathbf{L}}}(k))dydx 
= -\int_{J} \int_{J_k}  \varphi(y+E^{-\operatorname{\mathbf{L}}}(k))\varphi(x){\gamma}( h (y+E^{-\operatorname{\mathbf{L}}}(k)), h  x)dydx.\end{aligned}$$ Since $J_k = E^{-\operatorname{\mathbf{L}}}(k) + J $ by definition and $J = \dot\bigcup_{i=0}^{2^d} (\square -v_i )$ we obtain $$\begin{aligned}
{\texttt{dis}}_k
&=-2{ h }^{2d} \sum_{0 \leq \nu < 2^d} \sum_{0 \leq \mu < 2^d} \int_{\square} \int_{\square} \psi_\nu(x)\psi_\mu(y){\gamma}( h (v_\mu-E^{-\operatorname{\mathbf{L}}}(k)-v_\nu), h (y-x))dydx.\end{aligned}$$

All in all we conclude

$$\begin{aligned}
&{\texttt{dis}}_k/-2{ h }^{2d} \\
= &\begin{cases}
 \sum_{\nu \in {\left\{}
\newcommand{\rk}{\right\}}0,\ldots,  2^d\rk  \backslash D_i} \sum_{\mu \in  D_i^c} \int_\square \int_\square  \psi_\nu(x) \psi_\mu(y) {\gamma}( h (v_\mu - v_\nu - E^{-\operatorname{\mathbf{L}}}(k)), h (y-x)) dydx: &k = E^{\operatorname{\mathbf{L}}}(v_i)\\
  \sum_{0 \leq \nu < 2^d} \sum_{0 \leq \mu < 2^d} \int_{\square} \int_{\square} \psi_\nu(x)\psi_\mu(y){\gamma}( h (v_\mu-v_\nu-E^{-\operatorname{\mathbf{L}}}(k)), h (y-x))dydx: &else.
\end{cases}\end{aligned}$$

### Source term {#source-term .unnumbered}

We compute $b^{ h } \in {\mathbb{R}}^{\operatorname{\mathbf{L}}^d} $ by

$$\begin{aligned}
b^{ h }_k = \int_{\Omega} f {{\varphi_k}}dx
=h^d \sum_{\nu = 0}^{2^d - 1} \int_{\square} f(x_k +  h (v- v_\nu)) \psi_{\nu}(v)dv
=h^d 2^d \int_{\square} f(x_k  +  h  v) \psi_{0}(v)dv,\end{aligned}$$

 \
where the last equality follows again from considering $\psi_\nu \circ g_\nu = \psi_0$.

Solving procedure {#sec:solving}
=================

Now we discuss how to solve the discretized, fully populated multilevel Toeplitz system. The fundamental procedure uses an efficient implementation for the matrix-vector product of multilevel Toeplitz matrices, which is then delivered to the conjugate gradient (CG) method.

Let us first illuminate the implementation of the matrix-vector product $Tx , $ where $T \in {\mathbb{R}}^{\operatorname{\mathbf{L}}^d \times \operatorname{\mathbf{L}}^d}$ is a symmetric $d$-level Toeplitz matrix of order $\operatorname{\mathbf{L}}= (L_0, \ldots, L_{d-1})$ and $x$ a vector in ${\mathbb{R}}^{\operatorname{\mathbf{L}}^d}$. The crucial idea is to embed the Toeplitz matrix into a circulant matrix for which matrix-vector products can be efficiently computed with the help of the discrete fourier transform (DFT) [@structuredmatrices]. Here, by a $d$-level circulant matrix, we mean a matrix $C\in {\mathbb{R}}^{\operatorname{\mathbf{L}}^d \times \operatorname{\mathbf{L}}^d}$, which satisfies $$C_{\operatorname{\mathbf{i}}\operatorname{\mathbf{j}}} = C((\operatorname{\mathbf{i}}-\operatorname{\mathbf{j}})\mod \operatorname{\mathbf{L}}),$$ where $(\operatorname{\mathbf{i}}\mod \operatorname{\mathbf{L}}){\coloneqq}(i_k \mod L_k)_k$. In the real symmetric case, such that $a_{\operatorname{\mathbf{i}}\operatorname{\mathbf{j}}} = a(|\operatorname{\mathbf{i}}- \operatorname{\mathbf{j}}|)$ as it is present in our setting, $T$ can be reconstructed from its first row $R {\coloneqq}(T_{0i})_i \in {\mathbb{R}}^{\operatorname{\mathbf{L}}^d}$. Since circulant matrices are special Toeplitz matrices, the same holds for these matrices as well. Due to the multilevel structure it is convenient to represent $T$ by a tensor $t$ in ${\mathbb{R}}^{L_0 \times \cdots \times L_{d-1}}$, which is composed of the values contained in $R$. More precisely we define $$t(\operatorname{\mathbf{i}}) {\coloneqq}R_{E^{\operatorname{\mathbf{L}}}(\operatorname{\mathbf{i}})}$$ for $ \operatorname{\mathbf{i}}\in \prod_{i=0}^{d-1} {\left\{}
\newcommand{\rk}{\right\}}0,\ldots, L_i-1\rk$ with $E^{\operatorname{\mathbf{L}}}$ from above. Now $t$ can be embedded into the tensor representation $c\in {\mathbb{R}}^{2L_0 \times \cdots \times 2L_{d-1}}$ of the associated $d$-level circulant matrix by $$c(\operatorname{\mathbf{i}}) {\coloneqq}t(\hat{i}_0, \ldots, \hat{i}_{d-1}),$$ where $$\hat{i}_k {\coloneqq}\begin{cases}
i_k: &i_k<L_k, \\
0: &i_k = L_k, \\
2L_k-i_k: &else.
\end{cases}$$ We note that $t = c([0:L_0-1], \ldots, [0:L_{d-1}-1]) $. Thus, we can use Algorithm \[product\] to compute the product $T \cdot x$, where the DFT is carried out by the fast fourier transform (FFT).

 \
INPUT: $t \in {\mathbb{R}}^{L_0 \times \cdots \times L_{d-1}}$ representing $T \in {\mathbb{R}}^{\operatorname{\mathbf{L}}^d \times \operatorname{\mathbf{L}}^d}$, $x \in {\mathbb{R}}^{\operatorname{\mathbf{L}}^d} $\
OUTPUT: $y = Tx$\
\
1. Construct $c\in {\mathbb{R}}^{2L_0 \times \cdots \times 2L_{d-1}}$ by $c(\operatorname{\mathbf{i}}) {\coloneqq}t((\hat{i}_0, \ldots, \hat{i}_{d-1}) $\
2. Construct $x' \in {\mathbb{R}}^{2L_0 \times \cdots \times 2L_{d-1}}$ by $$x'(\operatorname{\mathbf{i}}){\coloneqq}\begin{cases}
    x_{E^L(\operatorname{\mathbf{i}})}: & \operatorname{\mathbf{i}}\in \prod_{i=0}^{d-1} {\left\{}
\newcommand{\rk}{\right\}}0,\ldots, L_i-1\rk \\
    0: else
    \end{cases}$$\
3. Compute $\Lambda = FFT_L(c) $\
4. Compute $z = FFT_L(x')$\
5. Compute $w = \Lambda  z $ (pointwise)\
6. Compute $y' = FFT^{-1}_L(w)$\
7. Construct $y \in {\mathbb{R}}^{L_0 \times \cdots \times L_{d-1}}$ by $y(\operatorname{\mathbf{i}}) = y'(\operatorname{\mathbf{i}})$ for $ \operatorname{\mathbf{i}}\in \prod_{i=0}^{d-1} {\left\{}
\newcommand{\rk}{\right\}}0,\ldots, L_i-1\rk $\
8. Return $y.reshape(\operatorname{\mathbf{L}}^d)$

In the Python code we use the library pyFFTW (<https://hgomersall.github.io/pyFFTW/>) to perform a parallelized multidimensional DFT, which is a pythonic wrapper around the C subroutine library FFTW (<http://www.fftw.org/>). Furthermore, we want to note, that an MPI implementation for solving multilevel Toeplitz systems in this fashion is presented in [@MPIproduct], which inspired us to apply the upper procedure.\
Finally, with this algorithm at hand, we employ a CG method, as it can be found for example in [@nocedale:numopt], to obtain the solution of the discretized system (\[discretizedproblem\]).

Numerical Experiments {#sec:numres}
=====================

In this last section we want to complete the previous considerations by presenting numerical results in 1d, 2d and for the first time also in 3d. We now specify the nonlocal diffusion operator $-{\mathcal{L}}$ by choosing the fractional kernel ${\gamma}(x,y) = \frac{c_{d,s}}{2||y-x||^{d+2s}}$ and shortly recall how the fractional Laplace operator $(-\Delta)^s$ crystallizes out as special case of $-{\mathcal{L}}$. Furthermore we describe in detail the employed numerical integration and finally discuss the results of the implementation.

Relation to space-fractional diffusion problems
-----------------------------------------------

We follow [@fractlapl] and define the action of the *fractional Laplace operator $(-\Delta)^s$* on a function ${u:{\mathbb{R}}^d \to {\mathbb{R}}}$ by $$\begin{aligned}
(-\Delta)^s u(x) {\coloneqq}c_{d,s} \int_{{\mathbb{R}}^d} \frac{u(x)- u(y)}{||y-x||^{d+2s}}dy,\end{aligned}$$ where $c_{d,s} {\coloneqq}s 2^{2s} \frac{\Gamma(\tfrac{d+2}{2})}{\Gamma(\tfrac{1}{2}) \Gamma(1-s)}  .$ The homogeneous steady-state *space-fractional diffusion problem* then reads as $$\begin{aligned}
(-\Delta)^s u(x) &= f(x), ~~~x \in \Omega, \nonumber\\
u(x) &= 0, ~~~~~~~ x \in \Omega^c \label{fractprob}.\end{aligned}$$ Thus, choosing the fractional kernel ${\gamma}(x,y) {\coloneqq}\frac{c_{d,s}}{2||y-x||^{d+2s}}\mathcal{X}_{B_R(x)}(y)$, which satisfies the conditions (\[case\_1\]) and (\[gen\_assumption\]), the nonlocal Dirichlet problem (\[prob\]) can be considered as a truncated version of problem (\[fractprob\]). In [@fractlapl] the authors show that the weak solution $u_R$ of the truncated problem (\[prob\]), for some interaction radius $R>0$, converges to the weak solution $u_\infty$ of (\[fractprob\]) as $R \to \infty$. We repeat the corresponding result stated in [@fractlapl].

Let $u_R \in V_c({\Omega \cup \Omega_I})$ and $u_\infty \in H_\Omega^s({\mathbb{R}}^d) {\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}u \in H^s({\mathbb{R}}^d) \colon u_{|\Omega^c} \equiv 0 \rk$ denote the weak solutions of (\[prob\]) and (\[fractprob\]) respectively. Then $$||u_\infty - u_R ||_{H^s({\Omega \cup \Omega_I})} \leq \frac{K_d}{C_1^2 s(R-I)^{2s}}||u_\infty||_{L^2(\Omega)},$$ where $I {\coloneqq}\min{\left\{}
\newcommand{\rk}{\right\}}L\in {\mathbb{R}}\colon \Omega \subset B_L(x)~~ \forall ~x \in \Omega \rk$, $C_1$ is the equivalence constant from above and $K_d$ is a constant depending only on the space dimension $d$.

See [@fractlapl Theorem 3.1].

Combining this result with the finite element estimate (\[fem\_estimate\_fract\]) we arrive at the estimate $$\begin{aligned}
||u_\infty - u_{R}^N  ||_{H^s({\Omega \cup \Omega_I})}  \leq \frac{K_d}{C_1^2 s(R-I)^{2s}}||u_\infty||_{L^2(\Omega)} + C h^{m+t-s}||u_R||_{H^{m+t}({\Omega \cup \Omega_I})}.\end{aligned}$$ Thus, the accuracy of finite element solutions for the weak formulation of problem (\[fractprob\]) depend on both, the discretization quality and the interaction horizon [@fractlapl].

Numerical computation of the integrals
--------------------------------------

In this subsection we want to point out how we numerically handle the occurring integrals.

### Nonsingular integrals

We mainly have to compute integrals of the form $
\int_{\square} \int_{\square} g(x,y) dydx \label{referenceintegral},
$ where $\square = [0,1]^d$ and $g \colon {\mathbb{R}}^d \times {\mathbb{R}}^d \to {\mathbb{R}}$ is a (typically smooth) function, which we assume to have no singularities in the domain $\square \times \square$. We approximate the value of this integral by employing a $n$-point Gauss-Legendre quadrature rule in each dimension. More precisely, we built a $d$-dimensional tensor grid ${\texttt{X}}\in {\mathbb{R}}^{d \times n^{d}}$ with associated weights $\texttt{W}^{single} \in {\mathbb{R}}^{n^d}$, such that $\int_{\square} g(x,y) dy \approx \sum_{i=0}^{n^{d}-1} g(x, {\texttt{X}}_i)  \texttt{W}^{single}_i$ for $x \in {\texttt{X}}$. Finally, we define the arrays $$\begin{aligned}
\texttt{V} &{\coloneqq}({\texttt{X}}, {\texttt{X}}, \ldots, {\texttt{X}}) \in {\mathbb{R}}^{d \times n^{2d}}, \\
\texttt{Q} &{\coloneqq}({\texttt{X}}_0, \ldots, {\texttt{X}}_0,{\texttt{X}}_1, \ldots, {\texttt{X}}_1, \ldots , {\texttt{X}}_{n^d - 1}, \ldots, {\texttt{X}}_{n^d - 1} ) \in {\mathbb{R}}^{d \times n^{2d}} \\end{aligned}$$ with associated weights $\texttt{W}^{double}\in {\mathbb{R}}^{n^{2d}}$, such that we finally arrive at the following quadrature rule: $$\begin{aligned}
\int_{\square} \int_{\square} g(x,y) dydx \approx \sum_{i=0}^{n^{2d}-1} g(\texttt{V}_i, \texttt{Q}_i)  \texttt{W}^{double}_i = (g(\texttt{V},\texttt{Q}) \cdot \texttt{W}^{double}).sum().\end{aligned}$$

### Singular integrals

As we have seen for the space-fractional diffusion problem, the kernel function may come along with singularities at $(x,y)$ with $x=y$. Therefore we start with a general observation, which paves the way for numerically handling these singularities. Let $f \colon {\mathbb{R}}^d \times {\mathbb{R}}^d \to {\mathbb{R}}$ be a symmetric function, i.e., $f(x,y) = f(y,x)$, and let us further define the sets $M {\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}(x,y) \in \square\times \square \colon y_d \in [0,x_d] \rk$ and $M' {\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}(x,y) \in \square\times \square \colon (y,x) \in M\rk.$ Then it is straightforward to show that $M \cup M' = \square \times \square$ and $$\begin{aligned}
M \cap M' = &{\left\{}
\newcommand{\rk}{\right\}}(x,y)\in \square \times \square \colon x_d = y_d \rk \\
= &{\left\{}
\newcommand{\rk}{\right\}}(v^{d-1},z,w^{d-1},z) \in {\mathbb{R}}^{2d} \colon (v^{d-1},w^{d-1},z) \in [0,1]^{2d-1} \rk,\end{aligned}$$ such that $\lambda_{2d}(M \cap M')= 0$. Hence, we find that $$\begin{aligned}
\int_{\square \times \square} f d\lambda_{2d}   &= \int_{M} f d\lambda_{2d} + \int_{M'} f d\lambda_{2d} - \int_{M \cap M'} f d\lambda_{2d}\\
&=\int_{M} f d\lambda_{2d} + \int_{M'} f d\lambda_{2d} .\end{aligned}$$ From the symmetry of $f$ we additionally deduce that $\int_{M} f d\lambda_{2d} = \int_{M'} f d\lambda_{2d} $ and therefore we finally obtain $$\begin{aligned}
\int_{\square \times \square} f d\lambda_{2d}   = 2 \int_{M} f d\lambda_{2d}.\end{aligned}$$ This observation can now be applied to the singular integrals occurring in the vector ${\texttt{sing}}$ such that $$\begin{aligned}
&\int_{\square}\int_{\square}((\psi_{i}(y)-\psi_{i}(x))(\psi_{0}(y)-\psi_{0}(x))){\gamma}( h  y, h  x)dydx   \\
= &2\int_{[0,1]^d}\int_{[0,1]^{d-1}\times [0,x_{d-1}]} ((\psi_{i}(y)-\psi_{i}(x))(\psi_{0}(y)-\psi_{0}(x))){\gamma}( h  y, h  x)dydx .\end{aligned}$$ The essential advantage of this representation relies on the fact that the singularities are now located on the boundary of the integration domain. Thus, we do not evaluate the integrand on its singularities while using quadrature points which lie in the interior. For the outer integral we employ the $d$-dimensional Gauss quadrature as pointed out above. Then for the inner integral with a fixed $x \in {\texttt{X}}$ we extend the one-dimensional adaptive (G7,K15)-Gauss-Kronrod quadrature rule to $d$-dimensional integrals by again tensorising the one-dimensional quadrature points. Moreover, in order to take full advantage of the Gauss-Kronrod quadrature, we divide the set $[0,1]^{d-1}\times [0,x_{d-1}]$ into $2^{d-1}$ disjoint rectangular subsets such that the singularity $x$ is located at a vertex. The latter partitioning reinforces the adaptivity property of the Gauss-Kronrod quadrature rule.

### Integrals with large interaction horizon

Now we discuss the quadrature of $$\begin{aligned}
P_0(x) =~ &\int_{ (J \cup J_k)^c} {\gamma}( h  x,  h  y)\mathcal{X}_{B_{R}( h  x)}( h  y) dy \\
=~&(1/h^d)\int_{ (I_0 \cup I_k )^c} {\gamma}( h  x, y-( a + h  ))\mathcal{X}_{B_{R}( h  x)}(y-( a + h  )) dy.\end{aligned}$$ We carry out two simplifications for the implementation, which are mainly motivated by the fact that ${\gamma}(x,y) \to 0$ as $y  \to \infty$ for kernels such as the fractional one. First, since $h \to 0$, we set $B_{R}( h  x) \equiv B_{R}(0)$. Especially when $R$ is large and $h$ small, this simplification does not significantly affect the value of the integral. Second, we employ the $||\cdot||_\infty$-norm for the ball $B_R(0)$ instead of the $||\cdot||_2$-norm. Hereby we also loose accuracy in the numerical integration but it simplifies the domain of integration in the sense that we can use our quadrature rules for rectangular elements. The latter simplification can additionally be justified by the fact that we want to model $R \to \infty$ for the fractional kernel and therefore have to truncate the $||\cdot||_2$-ball in any case. Consequently, we are concerned with the quadrature of the integral $$\begin{aligned}
\int_{  B_{R}^{||\cdot||_\infty}( a  +  h )   \backslash (I_0 \cup I_k )} {\gamma}( h  x, y-(a+ h )) dy\end{aligned}$$ for $x \in {\texttt{X}}$. For this purpose, we define the box $\mathcal{B}{\coloneqq}\prod_{i=0}^{d-1} [a_i-\lambda,a_i+\lambda]$ for a constant $R\geq \lambda > 2h$ such that $(I_0 \cup I_k ) \subset \mathcal{B}$ and we partition $$B_{R}^{||\cdot||_\infty}( a  +  h )  = \mathcal{B} \cup B_{R}( a  +  h ) \backslash \mathcal{B} .$$ Thus, we obtain $$\begin{aligned}
P_0(x)h^d&\approx\int_{  \mathcal{B}  \backslash (I_0 \cup I_k )} {\gamma}( h  x, y-( a  +  h )) dy + \int_{  B_{R}^{||\cdot||_\infty}( a  +  h )  \backslash  \mathcal{B}} {\gamma}( h  x, y-( a  +  h )) dy.\end{aligned}$$ We discretize $ \mathcal{B}$ with the same elements which we used for $\Omega$. This is convenient for two reasons. On the one hand we capture the critical values of the kernel, which in case of singular kernels typically decrease as $y\to \infty$. On the other hand, we have to leave out the integration over $(I_0 \cup I_k )$, which then can easily be implemented since we use the same discretization. However, this results in $\operatorname{\mathbf{N}}_2^d$, where $N_2^i {\coloneqq}2\lambda /  h $, hypercubes with base points $ y_j {\coloneqq}( a -\lambda e) +   h  E^{-\operatorname{\mathbf{N}}_2}(j)$ such that $$\begin{aligned}
\int_{   \mathcal{B}  \backslash (I_0 \cup I_k )} {\gamma}( h  x, y) dy&= \sum_{j=0,j \notin R_i}^{\operatorname{\mathbf{N}}_2^d -1} \int_{ y_j +  h \square} {\gamma}( h  x,  y-( a  +  h )) dy\\
&= h^d \sum_{j=0,j \notin R_i}^{\operatorname{\mathbf{N}}_2^d - 1} \int_{ \square} {\gamma}(\lambda e + h (e-E^{-\operatorname{\mathbf{N}}_2}(j)),  h (y-x)) dy,\end{aligned}$$ where $R_i$ contains the indices for those elements, which are contained in $(I_0 \cup I_k )$. This set can be characterized as $R_i {\coloneqq}{\left\{}
\newcommand{\rk}{\right\}}0\leq j < \operatorname{\mathbf{N}}_2^d - 1 \colon  y_j +  h \square \subset  (I_0 \cup I_k )\rk$. Note that by definition we have $I_0 =  a + h (e + \bigcup_{0 \leq \nu < 2^d} \square - v_\nu)$ and since $E^{-\operatorname{\mathbf{L}}}(k) = v_i$ for $k \in {\texttt{idx}}_s^k$ such that $x_k =  a + h (e +v_i)$ we know that $I_0^c \cap I_k =  a + h (e +v_i + \bigcup_{\nu \in   D_i^c} \square - v_\nu)$. Hence, $ y_j +  h \square \subset  (I_0 \cup I_k )$ if and only if $ y_j =  a + h (e - v_\nu)$ for $0 \leq \nu < 2^d$ or $ y_j =  a +   h (e +v_i- v_\nu)$ for $\nu \in   D_i^c$. Since $ y_j =  a -\lambda e + h  E^{-\operatorname{\mathbf{N}}_2}(j) $ we find by equating $ y_j$ with these requirements that $$\begin{aligned}
R_i ={\left\{}
\newcommand{\rk}{\right\}}E^{\operatorname{\mathbf{N}}_2}(e -  v_\nu +\lambda  h ^{-1} e) : 0 \leq \nu < 2^d \rk 
\cup {\left\{}
\newcommand{\rk}{\right\}}E^{\operatorname{\mathbf{N}}_2}(e +v_i -  v_\nu +\lambda  h ^{-1} e): \nu \in   D_i^c \rk.\end{aligned}$$ Now we discuss the quadrature of the second integral $\int_{  B_{R}^{||\cdot||_\infty}( a  +  h )  \backslash  \mathcal{B}} {\gamma}( h  x, y-( a  +  h )) dy$ for $x \in {\texttt{X}}$. Since we want to apply the algorithm to fractional diffusion, we have to get around the computational costs that occur when $R$ is large. In order to alleviate those costs we follow the idea in [@fractlapl] and apply a coarsening rule to discretize the domain $ B_{R}^{||\cdot||_\infty}( a  +  h ) \backslash \mathcal{B}$. Assume we have a procedure which outputs a triangulation $(z, \hat{ h })$ of $  B_{R}^{||\cdot||_\infty}( a  +  h )  \backslash \mathcal{B}$ consisting of $ N_3 $ hyperrectangles with base points $ z_j \in {\mathbb{R}}^d$ and sides of length $ \hat{ h }_j   \in {\mathbb{R}}^d$. Then we obtain $$\begin{aligned}
\int_{  B_{R}^{||\cdot||_\infty}( a  +  h )  \backslash \mathcal{B}} {\gamma}( h  x, y-( a  +  h )) dy &= \sum_{j=0}^{ N_3-1 } \int_{ z_j +  \hat{ h }_j \square} {\gamma}( h  x,  y-( a  +  h ) ) dy\\
&= \sum_{j=0}^{ N_3 -1} \hat{ h }_j ^d\int_{ \square} {\gamma}(- z_j+ a  +  h ,  \hat{ h }_j y - h  x ) dy.\end{aligned}$$ All in all we thus have

$$\begin{aligned}
    {\texttt{rad}}_k  &\approx 2|D_i|  { h }^{2d}\sum_{j=0,j \notin R_i}^{\operatorname{\mathbf{N}}_2^d-1} \int_{\square}\int_{ \square}\psi_0(x)\psi_{i}(x){\gamma}(\lambda e + h (e -E^{-\operatorname{\mathbf{N}}_2}(j)),  h (y-x)) dydx \\
    &+2|D_i| h ^{d} \sum_{j=0}^{ N_3-1 } \hat{ h }_j ^d\int_{\square}\int_{ \square}\psi_0(x)\psi_{i}(x){\gamma}(- z_j+ a  +  h ,  \hat{ h }_j  y- h  x) dydx.
    \end{aligned}$$

This leaves space for discussion concerning the choice of an optimal coarsening rule. We use the following preliminary approach in our code: We decompose $  B_{R}^{||\cdot||_\infty}( a  +  h )  \backslash \mathcal{B}$ into $(3^d-1)$ $d$-dimensional hyperrectangles surrounding the box $\mathcal{B}$. Then we build a tensor grid by employing in each dimension the coarsening strategy $ v + i^qh_{min} $ where $v$ is a vertex of $\mathcal{B}$, $q \geq 1$ the coarsening parameter and $h_{min}$ a minimum grid size. By concatenating all arrays we obtain a triangulation $(z, \hat{ h })$ of $  B_{R}^{||\cdot||_\infty}( a  +  h )  \backslash \mathcal{B}$.

Numerical results
-----------------

The implementation has been carried out in Python and the examples were run on a HP Workstation Z240 MT J9C17ET with Intel Core i7-6700 - 4 x 3.40GHz. Since we started from an arbitrary dimension throughout the whole analyzes, the codes for each dimension $d \in {\left\{}
\newcommand{\rk}{\right\}}1,2,3\rk$ own the same structure. Depending on the dimension, one only has to adapt the index sets ${\texttt{idx}}_s^i$, ${\texttt{idx}}_s^k$, $D_i$ and $D_i^c$, the implementation of the quadrature rules discussed above, the $2^d$ element basis functions and the plotting procedure. The rest can be implemented in a generic way. In addition to that, we framed the relevant representations for implementing the assembly of the first row. The implementation of the solving procedure only consists of delivering Algorithm \[product\] to the CG method. Moreover, the codes are parallelized over 8 threads on the four Intel cores. Within the assembling process of the first row, we first compute the more challenging $(d+1)$ entries $M_k = {\texttt{sing}}_k + {\texttt{rad}}_k + {\texttt{dis}}_k$ for $k \in {\texttt{idx}}_s^k$. Here we parallelize the computations of the integrals in ${\texttt{rad}}_k$ over the base points $y_j$ for the box and $z_j$ for the coarsening strategy. For the remaining $\operatorname{\mathbf{L}}^d - (d+1)$ indices we have that $M_k = {\texttt{dis}}_k$ and we can simply parallelize the loop over ${\left\{}
\newcommand{\rk}{\right\}}0\leq k < \operatorname{\mathbf{L}}^d \rk \backslash {\texttt{idx}}_s^k$. As mentioned above, the solving process is parallelized via the parallel fourier transform pyFFTW.

In all examples we use the fractional kernel $${\gamma}(x,y) = \frac{c_{d,s}}{2||y-x||_2^{d+2s}}\mathcal{X}_{B_R(x)}(y)$$ for $R=T+\lambda$ where $T=2^{10}$ with coarsening parameter $q=1.5$, minimum grid size $h_{min}=10^{-2}$ and a parameter $\lambda > 2h$ for the box $\mathcal{B}$. The CG method stops if a sufficient decrease of the residual $||Ax_k-b||/||b|| < 10^{-12}$ is reached. We present numerical examples for $d \in {\left\{}
\newcommand{\rk}{\right\}}1,2,3\rk$. For each grid size $ h $ we report on the number of grid points (“dofs”) and the number of CG iterations (“cg its”) as well as the CPU time (“CPU solving”) needed for solving the discretized system. Furthermore we compute the energy error $||u_R^h - u_\infty||_{H^s({\Omega \cup \Omega_I})}$, where $u_\infty$ is a numerical surrogate taken to be the finite element solution on the finest grid, and the rate of convergence.

1d {#d .unnumbered}
--

For the 1d example we choose the following set of parameters:

-   Domain: $a=-1$, $b=1$, $\Omega = [-1,1]$

-   Fraction: $s=0.7$

-   Source term: $f(x) \equiv 1$

-   Parameter for box: $\lambda = 5$

-   Gauss points: $n=7$ for unit interval $[0,1]$

The results are presented in Figure \[fig:1d\] and Table \[tab:1d\].

[ |c| r | r | c |c |c|c|]{} $h$ & dofs & cg its & energy error &rate &CPU solving \[s\]&CPU solving direct \[s\]\
$2^{-9}$ & 1,023 &201 &3.15e-04 & 0.71 &0.04 &0.00\
$2^{-10}$ & 2,047 &328 &1.86e-04 &0.71 &0.11 &0.01\
$2^{-11}$ & 4,095 &535 &1.08e-04 & 0.72 &0.20 &0.03\
$2^{-12}$ & 8,191 &932 &6.13e-05 & 0.72 &1.22 &0.12\
$2^{-13}$ & 16,383 &1,519 &3.34e-05 & 0.75 & 2.96 &0.50\
$2^{-14}$ & 32,767 & 2,475 &1.70e-05 & 0.79 & 13.99 &2.03\
$2^{-17}$ & 65,535 & 13,755 & - & - & 540.71 &143.22\

2d {#d-1 .unnumbered}
--

For the 2d example we choose the following set of parameters:

-   Domain: $a=(0,0)$, $b=(1,0.2)$, $\Omega = [0,1] \times [0,0.2]$

-   Fraction: $s=0.4$

-   Source term: $f(x) = x_0$

-   Parameter for box: $\lambda = 1$

-   Gauss points in each dimension: $n=6$, i.e., $36$ quadrature points for the unit square $[0,1]^2$

The results are presented in Figure \[fig:2d\] and Table \[tab:2d\].

[ |c | r | c| c |c |c|]{} $h$ & dofs & cg its& energy error &rate &CPU solving \[s\]\
$2^{-8}$ & 12,750 &40 & 4.10e-05 & 0.61 &0.07\
$2^{-9}$ & 51,611 &54 & 2.49e-05 & 0.60 &0.55\
$2^{-10}$ & 207,669 &72 & 1.42e-05 & 0.62 &2.33\
$2^{-11}$ & 835,176 &96 & 7.17e-06 & 0.70 &15.13\
$2^{-12}$ & 3,349,710 &128 & 2.17e-06 & 0.75 & 81.57\
$2^{-13}$ & 13,408,667 & 169 & - & - & 485.13\

3d {#d-2 .unnumbered}
--

For the 3d example we choose the following set of parameters:

-   Domain: $a=(0,0,0)$, $b=(1,1,1)$, $\Omega = [0,1]^3$

-   Fraction: $s= 0.4$

-   Source term: $f(x) \equiv 1$

-   Parameter for box: $\lambda = 0.5$

-   Gauss points in each dimension: $n=4$, i.e., $64$ quadrature points for the unit cube $[0,1]^3$

The results are presented in Figure \[fig:3d\] and Table \[tab:3d\].

[ |c | r |c| c |c |c|]{} $h$ & dofs & cg its & energy error &rate &CPU solving \[min\]\
$2^{-5}$ & 29,791 &21 & 1.81e-04 &0.77 &0.01\
$2^{-6}$ & 250,047 & 23 & 3.08e-05 & 0.71 &0.03\
$2^{-7}$ & 2,048,383 & 30 & 1.55e-05 & 0.73 &0.53\
$2^{-8}$ & 16,581,375 & 41 & 5.90e-06 & 0.75 &4.73\
$2^{-9}$ & 133,432,831 & 55 & - &- &36.28\

Discussion
----------

Concerning the number of CG iteration we want to point out several observations. First, we generally observed for various parameters that the number of CG iterations increases as we emphasize the singularity, i.e., as $s \to 1$. This coincides with Theorem 6.3 in [@wellposedness], stating that for shape-regular and quasi-uniform meshes the following estimate for the condition number $cond(A)$ holds: $$cond(A)\leq ch^{-2s},$$ where $c>0$ is a generic constant. A second general observation throughout various experiments was, that the ratio $\tfrac{k{(h/2)}}{k(h)}$, where $k(h)$ denotes the number of CG iterations needed for grid size $h$, depends only on $s$. More precisely, we observe that $\tfrac{k{(h/2)}}{k(h)} = c_s $, where $c_s$ is a constant independent of the dimension and the domain. For $s=0.4$ for instance we observe that $k(h/2) \approx 1.33 k(h)$, whereas for $s=0.7$ we find that $k(h/2) \approx 1.63 k(h)$. This is also reflected in the results above. However, since we consider a very “structured” case, where the domains are hyperrectangles triangulated with regular grids, we want to leave it an observation. Finally a third and last observation in this regard is the surprisingly low number of CG iterations needed for the 3d case. Running the code for different dimensions, but for a fixed comparable parameter setting, where we set the domain to be the unit hypercube $[0,1]^d$ and the source term to be $f \equiv 1$, we find that for problems of the same size, i.e., with the same number of degrees of freedom, the number of CG iterations decreases as the dimension increases. Specifically, let $s=0.4$, $\lambda = 0.5$ and $T=2^{10}$. Then in Table \[tab:compare\_dofs\] we find the results where the size of the discretized system is fixed to $dofs = \operatorname{\mathbf{L}}^d \approx 250000$.

[ |c | c |c|]{} $d$ & h & cg its\
$1$ & $1./250,048 \approx 2^{-18} $ &615\
$2$ & $2^{-9}$ & 62\
$3$ & $2^{-6} $ &23\

This might be due to the additional Toeplitz levels affecting the condition number of the stiffness matrix. In contrast to that, while fixing the grid size, the number of CG iterations seems to vary less comparing the different dimensions. In Table \[tab:compare\] the reader finds the results for a fixed grid size $h= 2^{-6}$.

[ |c | c |c|]{} $d$ & dofs & cg its\
$1$ & 63 &17\
$2$ & 3,969 & 24\
$3$ & 250,047 &23\

However, this has to be analyzed more concretely and a thorough investigation is not the intention at this point.

Besides this, we also want to mention how the assembling procedure is affected by the choice of the fraction $s$. For $s \to 1$ the treatment of the singularity becomes a more delicate issue and one has to increase the number of quadrature points $n$ to prevent instabilities. However, since the kernel decays more rapidly as $y\to \infty$ these costs can be balanced by choosing a smaller $\lambda$. In contrast to that, as $s \to 0$, a moderate number of quadrature points suffices, but $\lambda$ has to be chosen larger in order to prevent a degenerated solution.

Finally, we also note, that the library pyFFTW needs a lot of memory for building the FFT object, such that we had to move the 3d computations for the finest grid to a machine with a larger RAM. However, one can circumvent this problem by using the sequential FFT implementation available in the NumPy library.

Concluding remarks
==================

We presented a finite element implementation for the steady-state nonlocal Dirichlet problem with homogeneous volume constraints on an arbitrary $d$-dimensional hyperrectangle and for translation invariant kernel functions. We use a continuous Galerkin method with multilinear element basis functions and theoretically back up our numerics with the framework for nonlocal diffusion developed by Gunzburger et al. [@wellposedness]. The key result showing the multilevel Toeplitz structure of the stiffness matrix is proven for arbitrary dimension and paves the way for the first 3d implementations in this area. Furthermore, we comprehensively analyze the entries of the stiffness matrix and derive representations which can be efficiently implemented. Since throughout the whole analysis we start from an arbitrary dimension, we can almost generically implement the code; we mainly have to adapt the implementation of the quadrature rules, the element basis functions as well as the index sets.

An important extension of this work is to incorporate the case where the interaction horizon is smaller than the diameter of the domain. This complicates the integration and with that the assembling procedure, but the stiffness matrix is no more fully populated and its structure still remains multilevel Toeplitz. Having that, one can model the transition to local diffusion and access a greater range of kernels. The resulting code for 2d would then present a fast and efficient implementation, which could be used for example in image processing. Also, since certain kernel functions allow for solutions with jump discontinuities, also discontinuous Galerkin methods are conforming [@wellposedness]. In this case, one has to carefully analyze the structure of the resulting stiffness matrix, which might differ from a multilevel Toeplitz one. Moreover, an aspect concerning the solving procedure, which is not examined above, is that of an efficient preconditioner for the discretized Galerkin system. A multigrid method might be a reasonable candidate due to the simple structure of the grid (see also [@chen2016convergence; @PANG2012693]). In general, a lot of effort has been put in the research of preconditioning structured matrices (see e.g., [@olshans; @serra2007toeplitz; @chan1988optimal; @Capizzano2000AnyCP]). Since we observe a moderate number of CG iterations in our numerical examples, a preconditioner has not been implemented yet.

The main drawback of our approach relies on the fact that the code is strictly limited to regular grids and is thus not applicable to more complicated domains. It is crucial that each element has the same geometry in order to achieve the multilevel Toeplitz structure of the stiffness matrix; meaning that only rectangular domains are reasonable.

However, one could think about a coupling strategy, which allows to decompose a general domain into rectangular parts and the remaining parts. This is beyond the scope of this paper and left to future work. In contrast to that, the restriction to translation invariant kernels appears to be rather weak, since a lot of kernels treated in literature are even radial.

Acknowledgement {#acknowledgement .unnumbered}
===============

The first author has been supported by the German Research Foundation (DFG) within the Research Training Group 2126: ‘Algorithmic Optimization’.

[^1]: Universitaet Trier, D-54286 Trier, Germany, Email: vollmann@uni-trier.de, volker.schulz@uni-trier.de
