Summation of smooth and singular kernels

We are interested in the fast evaluation of linear combinations of radial functions, i.e. the computation of

$$g\left(\ensuremath{\boldsymbol{y}}_j\right) := \sum_{k=1}^N \alph... ...ensuremath{\boldsymbol{y}}_j-\ensuremath{\boldsymbol{x}}_k\right\Vert _2\right)$$

for $j=1,\hdots,M$ and nodes $\ensuremath{\boldsymbol{x}}_k, \ensuremath{\boldsymbol{y}}_j \in \ensuremath{\mathbb{R}}^d$ . For smooth kernels $K$ with an additional parameter $c>0$ , e.g. the Gaussian $K(x)={\,\rm {e}}^{-x^2/c^2}$ , the multiquadric $K(x)=\sqrt{x^2+c^2}$ or the inverse multiquadric $K(x)=1/\sqrt{x^2+c^2}$ our algorithm requires ${\cal O}(N+M)$ arithmetic operations. In the case of singular kernels $K$ , e.g.,

$$\frac{1}{x^2},\; \frac{1}{\vert x\vert}, \; \log \vert x\vert, \;... ...ert, \; \frac{1}{x}, \; \frac{\sin(c x)}{x},\; \frac{\cos(c x)}{x}, \; \cot(cx)$$

an additional regularisation procedure must be incorporated and the algorithm has the arithmetic complexity ${\cal O}(N\log N+M)$ or ${\cal O}(M\log M + N)$ if either the target nodes $\ensuremath{\boldsymbol{y}}_j$ or the source nodes $\ensuremath{\boldsymbol{x}}_k$ are ``reasonably uniformly distributed''.

Note that the proposed fast algorithm [49,50,23] generalises the diagonalisation of convolution matrices by Fourier matrices to the setting of arbitrary nodes. In particular, this yields nearly the same arithmetic complexity as the FMM [5] while allowing for an easy change of various kernels. A recent application in particle simulation is given in [43]. The directory applications/fastsum contains C and MATLAB programs that show how to use the fast summation method.