\section{Overlapped RBF-FD}
\label{sec:overlap}

\subsection{Mathematical description}
\label{sec:orbf_math}

We first present a mathematical description of the overlapped RBF-FD method, recently developed by the first author~\cite{ShankarJCP2017}. Let $X = \{\vx_k\}_{k=1}^N$ be a global set of nodes on the domain. Define the stencil $P_k$ to be the set of nodes containing node $\vx_{\calI^k_1}$ and its $n-1$ nearest neighbors $\{\vx_{\calI^k_2},\hdots,\vx_{\calI^k_n}\}$; here, $\{\calI^k_1,\hdots,\calI^k_n\}$ are indices that map into the global node set $X$. We defer discussion of the number of stencils to the end of this section. For the remainder of this discussion, we will focus without loss of generality on the stencil $P_1$. First, define the \emph{stencil width} $\rho_1$ as

\begin{align}
\rho_1 = \max\limits_{j} \|\vx_{\calI^1_1} - \vx_{\calI^1_j}\|, j=1,\hdots,n.
\label{eq:stencil_rad}
\end{align}

Choosing $\delta \in (0,1]$, we now define the \emph{stencil retention distance} $r_1$ to be

\begin{align}
r_1 = (1-\delta)\rho_1.
\label{eq:ball_rad}
\end{align}

Here, $\delta$ is called the \emph{overlap parameter}. Let $R_1$ be the set of global indices of the $p_1$ nodes in the stencil $P_1$ that are within the distance $r_1$ from $\vx_{\calI^1_1}$. We write $R_1$ as
\begin{align}
R_1 = \{\calR^1_1, \calR^1_2, \hdots, \calR^1_{p_1}\}.
\label{eq:ball_inds}
\end{align}
In general, $R_1$ is some permutation of a subset of the indices of the nodes in $P_1$. Let $\mathbb{B}_1$ contain the nodes whose indices are in $R_1$. Thus, 
\begin{align}
\mathbb{B}_1 = \{\vx_{\calR^1_1}, \hdots, \vx_{\calR^1_{p_1}} \}.
\label{eq:ball_pts}
\end{align}
For convenience, we will refer to $\mathbb{B}_1$ as a \emph{ball}, though it is merely a set of discrete points. The overlapped RBF-FD method involves computing RBF-FD weights for all the nodes in the ball $\mathbb{B}_1$, and repeating this process for each stencil $P_k$. The weights for all the nodes in $\mathbb{B}_1$ with indices in $R_1$ are computed using the following augmented local RBF interpolant on $P_1$:
\begin{align}
s_1(\vx,\vy) = \sum\limits_{j=1}^n w^1_j(\vy) \|\vx - \vx_{\calI^1_j}\|^m + \sum\limits_{i=1}^{M} \lambda^1_i(\vy) \psi^1_i(\vx),
\label{eq:rbf_interp}
\end{align}
where $\|\vx - \vx_{\calI^1_j}\|^m$ is the polyharmonic spline (PHS) RBF of degree $m$ ($m$ is odd), and $\psi^1_i(\vx)$ are the $M$ monomials corresponding to the total degree polynomial of degree $\ell$ in $d$ dimensions. Here, the $n$ overlapped RBF-FD weights associated with the point $\vy$ are $w^1_j(\vy), j=1,\hdots,n$. \revone{Thus, for the stencil $P_1$, we have two sets of points that are relevant to the interpolant: the set $P_1$ itself (with $\vx \in P_1$), and the set $\mathbb{B}_1$ (with $\vy \in \mathbb{B}_1$). In order to compute the weights for the linear operator $\calL$ uniquely at all nodes in $\mathbb{B}_1$ with indices in the set $R_1$, we impose the following two (sets of) conditions:
\begin{align}
\lf.s_1 \rt|_{\vx \in P_1, \vy \in \mathbb{B}_1} &= \lf.\calL \|\vx - \vx_{\calI^1_j}\|^m\rt|_{\vx \in \mathbb{B}_1}, j=1,\hdots,n, \label{eq:interp_constraint}\\
\sum_{j=1}^n \lf.w_j^1(\vy) \psi_i^1(\vx) \rt|_{\vx \in P_1, \vy \in \mathbb{B}_1} &= \lf.\calL \psi^1_i(\vx)\rt|_{\vx \in \mathbb{B}_1}, i=1,\hdots,M. \label{eq:poly_constraint}
\end{align}
}
The first set of conditions enforces that $s_1(\vx,\vy)$ interpolate the derivatives of the PHS RBF at all nodes in $\mathbb{B}_1$. The second set of conditions enforces polynomial reproduction on the overlapped RBF-FD weights. If $\ell$ is the degree of the appended polynomial $\psi^1(\vx)$, $M = {\ell + d \choose d}$; for stability, we also require that $M \leq \lfloor \frac{n}{2} \rfloor$~\cite{ShankarJCP2017,FlyerNS,FlyerPHS}. The interpolant \eqref{eq:rbf_interp} and the two conditions \eqref{eq:interp_constraint}--\eqref{eq:poly_constraint} can be collected into the following block linear system:
\begin{align}
\begin{bmatrix}
A_1 & \Psi_1 \\
\Psi_1^T & O
\end{bmatrix}
\begin{bmatrix}
W_1 \\
\Lambda_1
\end{bmatrix}
=
\begin{bmatrix}
B_{A_1} \\
B_{\Psi_1}
\end{bmatrix},
\label{eq:rbf_linsys}
\end{align}
where
\begin{align}
(A_1)_{ij} &= \|\vx_{\calI^1_i} - \vx_{\calI^1_j} \|^m, i,j=1,\hdots,n, \\
(\Psi_1)_{ij} &= \psi^1_j(\vx_{\calI^1_i}), i=1,\hdots,n, j=1,\hdots,M,\\
(B_{A_1})_{ij} &= \lf.\calL \|\vx - \vx_{\calR^1_j} \|^m \rt|_{\vx = \vx_{\calI^1_i}}, i=1,\hdots,n, j=1,\hdots,p_1, \\
(B_{\Psi_1})_{ij} &= \lf.\calL \psi^1_i(\vx)\rt|_{\vx = \vx_{\calR^1_j}}, i=1,\hdots,M, j=1,\hdots,p_1,\\
O_{ij} &= 0, i,j = 1,\hdots,M.
\end{align}
$W_1$ is the \emph{matrix} of overlapped RBF-FD weights, with each column containing the RBF-FD weights for a point $\vx \in \mathbb{B}_1$. The linear system \eqref{eq:rbf_linsys} has a unique solution if the nodes in $P_1$ are distinct~\cite{Fasshauer:2007,Wendland:2004}. The matrix of polynomial coefficients $\Lambda_1$ is merely a set of Lagrange multipliers that enforces the polynomial reproduction constraint \eqref{eq:poly_constraint}.

\revone{
\subsection{Algorithmic description of differentiation matrix assembly}
\label{sec:assembly}

\begin{algorithm}
\caption{Differentiation matrix assembly using Overlapped RBF-FD}
\label{alg:assembly}
\begin{algorithmic}[1]	
  \Statex{\bf Given}: $X = \{\vx_k\}_{k=1}^N$, the set of nodes in the domain.
	\Statex{\bf Given}: $\delta \in (0,1]$, the overlap parameter.
	\Statex{\bf Given}: $\calL$, the linear differential operator to be approximated.
	\Statex{\bf Given}: $n << N$, the stencil size.
	\Statex{\bf Generate}: $L$, the $N \times N$ differentiation matrix approximating $\calL$ on the set $X$.
	\Statex{\bf Generate}: $N_{\delta}$, the number of stencils.
	\State Build a kd-tree on the set $X$ in $O(N \log N)$ operations.
	\State Initialize $g$, an $N$-long array of flags, to 0.
	\State Initialize the stencil counter, $N_{\delta}=0$.
	\For {$k=1,N$}
		\If {g(k) == 0}		
				\State Use kd-tree to get $\{\vx_{\calI^k_1},\hdots,\vx_{\calI^k_n}\}$. Here, $\calI^k_1 = k$.
				\State Get $R_k = \{\calR^k_1, \hdots, \calR^k_{p_k}\}$ and $\mathbb{B}_k = \{\vx_{\calR^k_1}, \hdots, \vx_{\calR^k_{p_k}} \}$ using \eqref{eq:stencil_rad}--\eqref{eq:ball_pts}.
				\State Use \eqref{eq:rbf_linsys} to compute $W_k$, the $n \times p_k$ matrix of RBF-FD weights.
				\For {$i=1,p_k$}
					\State Set $g\lf(\calR^k_i\rt) = 1$.	
					\For{$j=1,n$}
							\State Set $L\lf(\calR^k_i,\calI^k_j\rt) = W_k(j,i)$.
					\EndFor
				\EndFor
				\State Set $N_{\delta} = N_{\delta}+1$.
		\EndIf
	\EndFor
\end{algorithmic}
\end{algorithm}
We next describe how to assemble an $N \times N$ sparse differentiation matrix from the mathematical description above. The algorithm described here is a \emph{serial} one; parallel algorithms for assembly are a subject of ongoing research. The procedure involves using a kd-tree for nearest neighbor searches, together with a form of breadth-first search, and is described in Algorithm \ref{alg:assembly}. For simplicity, algorithm \ref{alg:assembly} describes the assembly process as directly populating the entries of $L$. In practice, one stores the non-zero elements, and row and column indices/pointers (depending on the sparse matrix format in use), and uses these arrays in place of $L$.

It is also worth briefly examining the cost of this process. Since Algorithm \ref{alg:assembly} uses a kd-tree, it has a preprocessing cost of $O(N \log N)$. Each iteration in the outer loop of the algorithm (line $4$) incurs the following costs: $O\lf((n-1) \log N\rt)$ to find the nearest neighbors, and $O\lf( (n+M)^3 + p_k (n+M)^2 \rt)$ for finding the weight matrix $W_k$. There are $N_{\delta}$ such steps. If the points are quasi-uniform, then $N_{\delta} \approx \frac{N}{p}$, where $p = (1-\delta)^d n$~\cite{ShankarJCP2017}. If $\delta <1$, each iteration of this loop is more expensive than the standard RBF-FD method, but there are far fewer iterations (by a factor $\frac{1}{p}$). The assembly process as described above also requires fewer nearest neighbor searches than the standard RBF-FD method if $\delta <1$.
}

\subsection{Local error estimates}
\label{sec:error_estimate}

The polynomial reproduction property \eqref{eq:poly_constraint} of the augmented local RBF interpolant can be used to develop a local error estimate~\cite{DavydovMinimal2016}. Following~\cite{ShankarJCP2017}, if $\calL$ is a differential operator of order $\theta$ applied to some function $f$, this error estimate can be written as:
{\small
\begin{align}
\lf| \calL f(\vy) - \sum\limits_{j=1}^n w_j(\vy) f(\vx_j) \rt| \leq P(\vy)\max\limits_{|\alpha| = \ell} \|D^{\alpha} f\|_{W^{1,\infty}(\Omega_1)} \lf( h(\vy) \rt)^{\ell+1-\theta}, 
\label{eq:error_estimate}
\end{align}
}
where $\vy$ is some evaluation point, $h(\vy)$ is the largest distance between the point $\vy$ and every point $\vx$ in the stencils, $P(\vy)$ is some \emph{growth function}~\cite{DavydovMinimal2016}, and $\alpha$ is a multiindex. The term $\max\limits_{|\alpha| = \ell} \|D^{\alpha} f\|_{W^{1,\infty}(\Omega_1)}$ is simply the Sobolev $\infty$-norm of $f$ on the stencil $P_1$ with convex hull $\Omega_1$. The error estimate shows \revone{shows that the order of convergence in any RBF-FD method depends primarily on the degree $\ell$, not the location of the evaluation point $\vy$. One can therefore expect the same orders of convergence for the overlapped RBF-FD method as in the standard one, assuming $\vy$ is not too close to the boundary of $\Omega_1$. As was shown in~\cite{ShankarJCP2017}, the difference in accuracy between the standard and overlapped RBF-FD methods is virtually non-existent for higher-order approximations, can be kept small for lower-order approximations using judicious choices of $\delta$.}

\subsection{Parameter selection}

We now describe how the different parameters in the overlapped RBF-FD method are selected.  Based on \eqref{eq:error_estimate}, we know the polynomial degree $\ell$ controls the order of the approximation for a differential operator of order $\theta$. Thus, given a differential operator of order $\theta$, if we require an RBF-FD method with order of accuracy $\xi$, we set
\begin{align}
\ell = \xi + \theta - 1.
\label{eq:select_ell}
\end{align}
For instance, if we are approximating the gradient operator $\nabla$ ($\theta = 1$) to third-order accuracy ($\xi = 3$), this gives $\ell = 3$. On the other hand, if we are approximating the Laplacian operator $\Delta$ ($\theta = 2$) to third-order accuracy ($\xi = 3$), this gives $\ell = 4$. The stencil size $n$ in $d$ dimensions can be deduced from the relationships $M = {\ell + d \choose d}$ and $M \leq \lf \lfloor\frac{n}{2} \rt\rfloor$ to be $n \geq 2M$. In practice, it appears to be beneficial to use a few more nodes per stencil than $2M$, though the optimal choice is likely problem dependent. We use
\begin{align}
n = 2M + \lf \lfloor \ln(2M) \rt\rfloor.
\label{eq:select_n}
\end{align}
We must also choose the overlap parameter $\delta \in (0,1]$. \revone{For small values of $n$ (equivalently of $\ell$), this choice trades off speedup with accuracy as $\delta \to 0$; however, as $n$ increases, smaller values of $\delta>0.2$ barely impact accuracy~\cite{ShankarJCP2017}}. In practice, we have observed that setting $\delta \leq 0.2$ typically completely decouples the stencils, resulting in ill-posed subproblems when solving PDEs. Given these constraints, we use the following heuristic:
\[
 \delta =
  \begin{cases} 
      \hfill 0.7 \hfill & \text{ if $\ell \leq 3$} \\
      \hfill 0.5 \hfill & \text{ if $4 \leq \ell < 6$} \\
			\hfill 0.3 \hfill & \text{ if $\ell \geq 6$} \\
  \end{cases}
\]
The automatic selection of $\delta$ based on stability properties of augmented local RBF interpolants is a subject of ongoing research. It is likely that the optimal choice of $\delta$ depends on the local Lebesgue functions corresponding to the differential operator $\calL$; see Section \ref{sec:stab} and~\cite{ShankarJCP2017,DavydovMinimal2016} for discussions on these local Lebesgue functions.

\revone{
\textbf{Connection to other RBF methods}
As described here, the overlapped RBF-FD method bears some resemblance to both the standard RBF-FD method and the RBF-Partition of Unity (RBF-PU) method~\cite{Wendland:2004,Fasshauer:2007}. In the latter case, one computes a \emph{global} interpolant to the node set $X$ by dividing $X$ into a set of overlapping partitions called \emph{patches} (controlled by an overlap parameter), forming a local RBF approximation on each patch, and then blending these approximations together with weight functions that are compactly-supported on each patch. As such, the standard RBF-FD method can be thought of as a very special RBF-PU method with weight functions that are exactly 1 at each node, and zero elsewhere, resulting in as many patches as there are nodes. The overlapped RBF-FD method can then be thought of as an RBF-PU method where the weight functions are 1 on the balls $\mathbb{B}_k$, zero elsewhere, resulting in fewer patches than nodes. However, in our opinions, the overlapped RBF-FD method is best viewed as an RBF-FD method where one makes a choice of either centered or one-sided approximations for each node in $X$, thereby resulting in stencil-sharing for some nodes. Indeed, the overlapped and standard RBF-FD methods share the same local error estimates from Section \ref{sec:error_estimate}. That said, it may be that viewing the assembly process from the RBF-PU perspective may help with parallelizing our method. This is a topic of ongoing research.
}
{"config": "arxiv", "file": "1806.03798/RBFReview.tex"}