Fast multipole method

In its simplest form, the fast multipole method seeks to evaluate the following function:

${\displaystyle f(y)=\sum _{\alpha =1}^{N}{\frac {\phi _{\alpha }}{y-x_{\alpha }}}}$,

where ${\displaystyle x_{\alpha }\in [-1,1]}$ are a set of poles and ${\displaystyle \phi _{\alpha }\in \mathbb {C} }$ are the corresponding pole weights on a set of points ${\displaystyle \{y_{1},\ldots ,y_{M}\}}$ with ${\displaystyle y_{\beta }\in [-1,1]}$. This is the one-dimensional form of the problem, but the algorithm can be easily generalized to multiple dimensions and kernels other than ${\displaystyle (y-x)^{-1}}$.

Naively, evaluating ${\displaystyle f(y)}$ on ${\displaystyle M}$ points requires ${\displaystyle {\mathcal {O}}(MN)}$ operations. The crucial observation behind the fast multipole method is that if the distance between ${\displaystyle y}$ and ${\displaystyle x}$ is large enough, then ${\displaystyle (y-x)^{-1}}$ is well-approximated by a polynomial. Specifically, let ${\displaystyle -1<t_{1}<\ldots <t_{p}<1}$ be the Chebyshev nodes of order ${\displaystyle p\geq 2}$ and let ${\displaystyle u_{1}(y),\ldots ,u_{p}(y)}$ be the corresponding Lagrange basis polynomials. One can show that the interpolating polynomial:

${\displaystyle {\frac {1}{y-x}}=\sum _{i=1}^{p}{\frac {1}{t_{i}-x}}u_{i}(y)+\epsilon _{p}(y)}$

converges quickly with polynomial order, ${\displaystyle |\epsilon _{p(y)}|<5^{-p}}$, provided that the pole is far enough away from the region of interpolation, ${\displaystyle |x|\geq 3}$ and ${\displaystyle |y|<1}$. This is known as the "local expansion".

The speed-up of the fast multipole method derives from this interpolation: provided that all the poles are "far away", we evaluate the sum only on the Chebyshev nodes at a cost of ${\displaystyle {\mathcal {O}}(Np)}$, and then interpolate it onto all the desired points at a cost of ${\displaystyle {\mathcal {O}}(Mp)}$:

${\displaystyle \sum _{\alpha =1}^{N}{\frac {\phi _{\alpha }}{y_{\beta }-x_{\alpha }}}=\sum _{i=1}^{p}u_{i}(y_{\beta })\sum _{\alpha =1}^{N}{\frac {1}{t_{i}-x_{\alpha }}}\phi _{\alpha }}$

Since ${\displaystyle p=-\log _{5}\epsilon }$, where ${\displaystyle \epsilon }$ is the numerical tolerance, the total cost is ${\displaystyle {\mathcal {O}}((M+N)\log(1/\epsilon ))}$.

To ensure that the poles are indeed well-separated, one recursively subdivides the unit interval such that only ${\displaystyle {\mathcal {O}}(p)}$ poles end up in each interval. One then uses the explicit formula within each interval and interpolation for all intervals that are well-separated. This does not spoil the scaling, since one needs at most ${\displaystyle \log(1/\epsilon )}$ levels within the given tolerance.

• Barnes–Hut simulation
• Multipole expansion
• n-body simulation

• Gibson, Walton C. The Method of Moments in Electromagnetics (3rd ed.), Chapman & Hall/CRC, 2021. ISBN 9780367365066.
• Abstract of Greengard and Rokhlin's original paper
• A short course on fast multipole methods by Rick Beatson and Leslie Greengard.
• JAVA Animation of the Fast Multipole Method Nice animation of the Fast Multipole Method with different adaptations.