A simplified fast multipole method based on strong recursive skeletonization

Abstract

This work introduces a kernel-independent, multilevel, adaptive algorithm for efficiently evaluating a discrete convolution kernel with a given source distribution. The method is based on linear algebraic tools such as low rank approximation and “skeleton representations” to approximate far-field interactions. While this work is related to previous linear algebraic formulations of the fast multipole method, the proposed algorithm is distinguished by relying on simpler data structures.

The proposed algorithm eliminates the need for explicit interaction lists by restructuring computations to operate exclusively on the near-neighbor list at each level of the tree, thereby simplifying both implementation and data structures. This work also introduces novel translation operators that significantly simplify the handling of adaptive point distributions. As a kernel-independent approach, it only requires evaluation of the kernel function, making it easily adaptable to a variety of kernels. By using operations on the neighbor list (of size at most 27 in 3D) rather than the interaction list (of size up to 189 in 3D), the algorithm is particularly well-suited for parallel implementation on modern hardware.

Numerical experiments on uniform and non-uniform point distributions in 2D and 3D demonstrate the effectiveness of the proposed parallel algorithm for Laplace and (low-frequency) Helmholtz kernels. The algorithm constructs a tailored skeleton representation for the given geometry during a precomputation stage. After precomputation, the fast summation achieves high efficiency on the GPU using batched linear algebra operations.