Fast multipole method

<h2 id="sketch-of-the-algorithm">Sketch of the Algorithm</h2>

<p>In its simplest form, the fast multipole method seeks to evaluate the following function:
</p>

f
        (
        y
        )
        =
        
          ∑
          
            α
            =
            1
          
          
            N
          
        
        
          
            
              ϕ
              
                α
              
            
            
              y
              −
              
                x
                
                  α
                
              
            
          
        
      
    
    {\displaystyle f(y)=\sum _{\alpha =1}^{N}{\frac {\phi _{\alpha }}{y-x_{\alpha }}}}
  
,
<p>where 
  
    
      
        
          x
          
            α
          
        
        ∈
        [
        −
        1
        ,
        1
        ]
      
    
    {\displaystyle x_{\alpha }\in [-1,1]}
  
 are a set of poles and 
  
    
      
        
          ϕ
          
            α
          
        
        ∈
        
          C
        
      
    
    {\displaystyle \phi _{\alpha }\in \mathbb {C} }
  
 are the corresponding pole weights on a set of points 
  
    
      
        {
        
          y
          
            1
          
        
        ,
        …
        ,
        
          y
          
            M
          
        
        }
      
    
    {\displaystyle \{y_{1},\ldots ,y_{M}\}}
  
 with 
  
    
      
        
          y
          
            β
          
        
        ∈
        [
        −
        1
        ,
        1
        ]
      
    
    {\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 
  
    
      
        (
        y
        −
        x
        
          )
          
            −
            1
          
        
      
    
    {\displaystyle (y-x)^{-1}}
  
.
</p><p>Naively, evaluating 
  
    
      
        f
        (
        y
        )
      
    
    {\displaystyle f(y)}
  
 on 
  
    
      
        M
      
    
    {\displaystyle M}
  
 points requires 
  
    
      
        
          
            O
          
        
        (
        M
        N
        )
      
    
    {\displaystyle {\mathcal {O}}(MN)}
  
 operations.  The crucial observation behind the fast multipole method is that if the distance between 
  
    
      
        y
      
    
    {\displaystyle y}
  
 and 
  
    
      
        x
      
    
    {\displaystyle x}
  
 is large enough, then 
  
    
      
        (
        y
        −
        x
        
          )
          
            −
            1
          
        
      
    
    {\displaystyle (y-x)^{-1}}
  
 is well-approximated by a <a href="/facts/Polynomial/Lzak8VVx">polynomial</a>. Specifically, let 
  
    
      
        −
        1
        <
        
          t
          
            1
          
        
        <
        …
        <
        
          t
          
            p
          
        
        <
        1
      
    
    {\displaystyle -1<t_{1}<\ldots <t_{p}<1}
  
 be the <a href="/facts/Chebyshev_nodes/4onexv5e">Chebyshev nodes</a> of order 
  
    
      
        p
        ≥
        2
      
    
    {\displaystyle p\geq 2}
  
 and let 
  
    
      
        
          u
          
            1
          
        
        (
        y
        )
        ,
        …
        ,
        
          u
          
            p
          
        
        (
        y
        )
      
    
    {\displaystyle u_{1}(y),\ldots ,u_{p}(y)}
  
 be the corresponding <a href="/facts/Lagrange_polynomial/MnpMmhtC">Lagrange basis polynomials</a>. One can show that the interpolating polynomial:
</p>

1
            
              y
              −
              x
            
          
        
        =
        
          ∑
          
            i
            =
            1
          
          
            p
          
        
        
          
            1
            
              
                t
                
                  i
                
              
              −
              x
            
          
        
        
          u
          
            i
          
        
        (
        y
        )
        +
        
          ϵ
          
            p
          
        
        (
        y
        )
      
    
    {\displaystyle {\frac {1}{y-x}}=\sum _{i=1}^{p}{\frac {1}{t_{i}-x}}u_{i}(y)+\epsilon _{p}(y)}

<p>converges quickly with polynomial order, 
  
    
      
        
          |
        
        
          ϵ
          
            p
            (
            y
            )
          
        
        
          |
        
        <
        
          5
          
            −
            p
          
        
      
    
    {\displaystyle |\epsilon _{p(y)}|<5^{-p}}
  
, provided that the pole is far enough away from the region of interpolation, 
  
    
      
        
          |
        
        x
        
          |
        
        ≥
        3
      
    
    {\displaystyle |x|\geq 3}
  
 and 
  
    
      
        
          |
        
        y
        
          |
        
        <
        1
      
    
    {\displaystyle |y|<1}
  
. This is known as the "local expansion".
</p><p>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 
  
    
      
        
          
            O
          
        
        (
        N
        p
        )
      
    
    {\displaystyle {\mathcal {O}}(Np)}
  
, and then interpolate it onto all the desired points at a cost of 
  
    
      
        
          
            O
          
        
        (
        M
        p
        )
      
    
    {\displaystyle {\mathcal {O}}(Mp)}
  
:
</p>

∑
          
            α
            =
            1
          
          
            N
          
        
        
          
            
              ϕ
              
                α
              
            
            
              
                y
                
                  β
                
              
              −
              
                x
                
                  α
                
              
            
          
        
        =
        
          ∑
          
            i
            =
            1
          
          
            p
          
        
        
          u
          
            i
          
        
        (
        
          y
          
            β
          
        
        )
        
          ∑
          
            α
            =
            1
          
          
            N
          
        
        
          
            1
            
              
                t
                
                  i
                
              
              −
              
                x
                
                  α
                
              
            
          
        
        
          ϕ
          
            α
          
        
      
    
    {\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 }}

<p>Since 
  
    
      
        p
        =
        −
        
          log
          
            5
          
        
        ⁡
        ϵ
      
    
    {\displaystyle p=-\log _{5}\epsilon }
  
, where 
  
    
      
        ϵ
      
    
    {\displaystyle \epsilon }
  
 is the numerical tolerance, the total cost is 
  
    
      
        
          
            O
          
        
        (
        (
        M
        +
        N
        )
        log
        ⁡
        (
        1
        
          /
        
        ϵ
        )
        )
      
    
    {\displaystyle {\mathcal {O}}((M+N)\log(1/\epsilon ))}
  
.
</p><p>To ensure that the poles are indeed well-separated, one recursively subdivides the unit interval such that only 
  
    
      
        
          
            O
          
        
        (
        p
        )
      
    
    {\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 
  
    
      
        log
        ⁡
        (
        1
        
          /
        
        ϵ
        )
      
    
    {\displaystyle \log(1/\epsilon )}
  
 levels within the given tolerance.
</p>
<h2 id="see-also">See also</h2>
<ul><li><a href="/facts/Barnes%E2%80%93Hut_simulation/bsxHpaYb">Barnes–Hut simulation</a></li>
<li><a href="/facts/Multipole_expansion/Kao2WRBA">Multipole expansion</a></li>
<li><a href="/facts/N-body_simulation/xrHwX7j3"><i>n</i>-body simulation</a></li></ul>

<h2 id="further-readings">Further readings</h2>
<ul><li>Yijun Liu: <i>Fast Multipole Boundary Element Method: Theory and Applications in Engineering</i>, Cambridge Univ. Press, ISBN 978-0-521-11659-6 (2009).</li></ul>
<h2 id="external-links">External links</h2>
<ul><li>Gibson, Walton C. <i>The Method of Moments in Electromagnetics</i> (3rd ed.), Chapman & Hall/CRC, 2021. <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 9780367365066.</li>
<li><a href="http://portal.acm.org/citation.cfm?id=36901">Abstract of Greengard and Rokhlin's original paper</a></li>
<li><a href="http://math.nyu.edu/faculty/greengar/shortcourse_fmm.pdf">A short course on fast multipole methods</a> by Rick Beatson and Leslie Greengard.</li>
<li><a href="http://www.umiacs.umd.edu/~ramani/cmsc878R/fmmdemo2.jar">JAVA Animation of the Fast Multipole Method</a> Nice animation of the Fast Multipole Method with different adaptations.</li></ul>
<h3>Free software</h3>
<ul><li><a href="http://sourceforge.net/projects/puma-em/">Puma-EM</a> A high performance, parallelized, open source Method of Moments / Multilevel Fast Multipole Method electromagnetics code.</li>
<li><a href="http://www.harperlangston.com/kifmm3d/documentation/index.html">KIFMM3d</a> The Kernel-Independent Fast Multipole 3d Method (kifmm3d) is a new FMM implementation which does not require the explicit multipole expansions of the underlying kernel, and it is based on kernel evaluations.</li>
<li><a href="https://github.com/dmalhotra/pvfmm">PVFMM</a> A optimized parallel implementation of KIFMM for computing potentials from particle and volume sources.</li>
<li><a href="http://www.yijunliu.com/Software">FastBEM</a> Free fast multipole boundary element programs for solving 2D/3D potential, elasticity, stokes flow and acoustic problems.</li>
<li><a href="http://www.fastfieldsolvers.com">FastFieldSolvers</a> maintains the distribution of the tools, called FastHenry and FastCap, developed at M.I.T. for the solution of Maxwell equations and extraction of circuit parasites (inductance and capacitance) using the FMM.</li>
<li><a href="https://github.com/exafmm/exafmm">ExaFMM</a> ExaFMM is a CPU/GPU capable 3D FMM code for Laplace/Helmholtz kernels that focuses on parallel scalability.</li>
<li><a href="http://scalfmm-public.gforge.inria.fr/doc/">ScalFMM</a> <a href="https://web.archive.org/web/20170502145137/http://scalfmm-public.gforge.inria.fr/doc/">Archived</a> 2017-05-02 at the <a href="/facts/Wayback_Machine/nmQ3a6JC">Wayback Machine</a> ScalFMM is a C++ software library developed at <a href="/facts/Inria/ruIkcapf">Inria</a> Bordeaux with high emphasis on genericity and parallelization (using <a href="/facts/OpenMP/69IpI80o">OpenMP</a>/<a href="/facts/Message_Passing_Interface/BEuBgDpm">MPI</a>).</li>
<li><a href="https://jacksondebuhr.github.io/dashmm/">DASHMM</a> DASHMM is a C++ Software library developed at Indiana University using Asynchronous Multi-Tasking HPX-5 runtime system. It provides a unified execution on shared and distributed memory computers and provides 3D Laplace, Yukawa, and Helmholtz kernels.</li>
<li><a href="https://zhang416.github.io/recfmm/">RECFMM</a> Adaptive FMM with dynamic parallelism on multicores.</li></ul>

<h2 id="references">References</h2>

<ol>
<li id="fn:1"><p>Rokhlin, Vladimir (1985). "Rapid Solution of Integral Equations of Classic Potential Theory." J. Computational Physics Vol. 60, pp. 187–207. <a href="https://www.sciencedirect.com/science/article/pii/0021999185900026" target="_blank">https://www.sciencedirect.com/science/article/pii/0021999185900026</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></p></li>
<li id="fn:2"><p>Nader Engheta, William D. Murphy, Vladimir Rokhlin, and Marius Vassiliou (1992), “The Fast Multipole Method for Electromagnetic Scattering Computation,” IEEE Transactions on Antennas and Propagation 40, 634–641. <a href="/wiki/Nader_Engheta" target="_blank">/wiki/Nader_Engheta</a> <a href="#fnref:2" class="footnote-back-ref">↩</a></p></li>
<li id="fn:3"><p>"The Fast Multipole Method". Archived from the original on 2011-06-03. Retrieved 2010-12-10.
 <a href="https://web.archive.org/web/20110603231158/http://www-theor.ch.cam.ac.uk/people/ross/thesis/node97.html" target="_blank">https://web.archive.org/web/20110603231158/http://www-theor.ch.cam.ac.uk/people/ross/thesis/node97.html</a> <a href="#fnref:3" class="footnote-back-ref">↩</a></p></li>
<li id="fn:4"><p>Cipra, Barry Arthur (May 16, 2000). "The Best of the 20th Century: Editors Name Top 10 Algorithms". SIAM News. 33 (4). Society for Industrial and Applied Mathematics: 2. Retrieved February 27, 2019. <a href="/wiki/Barry_Arthur_Cipra" target="_blank">/wiki/Barry_Arthur_Cipra</a> <a href="#fnref:4" class="footnote-back-ref">↩</a></p></li>
</ol>

Fast multipole method open-in-new

Fast multipole method