椰蓉面包，Barnes-Hut 算法和快速多极子方法（Fast Multipole Method）

八六最近刷到一个视频讲FMM的，链接：https://www.youtube.com/watch?v=FhMftauQZqU

想到最近八六发文都太水了，所以我们今天来介绍一下FMM~

（没错这种语气是看筱泽广讲解物理学合同志学的~小筱泽太棒了）

（关键就在于说话先不说主语吧，对吧？）

（看完这篇文给我去看筱泽物理学讲义啊给我去看，第一部已经汉化了喜：https://e-hentai.org/g/3983370/2f2b6c0b0c/）

那为了模仿合同志里的小短剧讲解方式，我们也来做吧~

正文开始

---

（下午1点04分，学院走廊里）

八六：诶，这不是evilKT同学吗，好巧啊~

Evil：八六同学！叫我evil就好啦~看这个方向，又是要去买面包来糊弄过去午餐了吧

八六：啊哈哈，真是什么都瞒不过evil啊……

Evil：那我也陪你去吧~

八六：嗯！听说有新口味的面包，一起尝尝吧~

（一段路上温馨的对话后，学院面包店门口）

Evil：啊！刚出炉的椰蓉面包！

八六：哇！难道evil也是椰蓉爱好者？

Evil：难道八六也？！

八六：嗯！老板，请来两个椰蓉面包！

（两人拿着面包，坐在了餐桌两侧）

八六：说起来，evil同学，你不觉得这些椰蓉很像一个个小小的小惑星吗？

Evil：？嘛，八六同学这么一说的话好像是有点……

八六：呵呵…要是这些椰蓉之间的引力需要用计算机模拟，evil同学该怎么办呢？

Evil：诶？突然问这种问题……那只能利用牛顿先生的万有引力定律对每一对椰蓉之间的引力进行计算来获得每个椰蓉受到的合外力，再利用这个力来确定各个椰蓉的加速度……

八六：嗯！很直接的想法呢！倒不如说感觉不会有其他的方法了吧！

Evil：什么嘛！每次八六同学你这样说就肯定是对我藏着掖着！快说吧，肯定有更好的方法吧！

八六：要是要求完全符合引力定律的话确实没有哦？

Evil：诶，怎么这样……

八六：（邪魅的笑容）不过，要是可以允许一点点误差的话，更快的方法也是有的！不如说既然决定使用计算机进行模拟的话数值计算误差比如浮点数误差就是不可避免的吧~

Evil：我就知道！八六同学真是坏心眼！

八六：好啦好啦！正好最近看到了有趣的方法，名字叫做快速多极子方法（Fast Multipole Method）

Evil：快速？多级子？何意味？

八六：啊！面包都要凉了！我们等会去自习室再说吧~（双手合十）我开动了~

Evil：啊啊！果然八六还是坏心眼！

Barnes-Hut 算法

---

（享用了椰蓉面包的两人，一路上聊着天，到了无人的自习室内）

八六：啊，没有人呢。

Evil：当然的吧，现在是午休时间呢！

八六：啊，那正好给evil同学助眠呢，讲解算法的话~

Evil：不要太欺负我啦……

八六：呵呵~evil同学还记得我们要讨论什么算法吗？

Evil：不是快速多极子算法吗？

八六：嗯！蒸蚌！不过在那之前，我想先讲讲另一个算法呢，叫做Barnes-Hut 算法。

Evil：又来了！奇怪的算法名字！

八六：嗯？有这么奇怪吗？虽然确实是人名命名的，没有FMM的名字这么易懂呢~

Evil：原来在八六的心里，快速多极子这种名词很易懂吗……

八六：呵呵~evil同学这么一说，不是显得我很喜欢炫耀吗？

Evil：事实如此吧？！

八六：嘛，回到我们的话题，这么多的椰蓉，要是按照引力定律平方反比的式子来说，面包上两个角的椰蓉之间的引力会比相对更近的椰蓉之间的引力小吧？

Evil：（完全没有反省自己啊八六同学！）

八六：所以说，对于我们一个角上的这个椰蓉，要是我们把对角这一块的椰蓉都当成一个大椰蓉团，再计算角上的椰蓉被这个大椰蓉团所吸引的引力，就可以省去很多对椰蓉之间的计算吧？

Evil：啊……有点听不懂了……

八六：嗯……（环顾四周）啊！自习室有白板真是帮大忙了。

（八六在白板上画着什么，Evil则在旁边静静地看着）

八六：怎么样，看的还清楚吗？

Evil：嗯！有图就好理解多了~

Evil：不过这椰蓉面包也太方方正正，椰蓉也太少了吧~

八六：看来是黑心老板呢~

八六：就像这样，将椰蓉面包分成四份，再把小份的面包再分成四份，这样的操作哦。

Evil：嗯？那A和H的椰蓉为什么还是一个大块呢？

八六：是因为当一份面包里只剩一个椰蓉的话，再分下去也只会多出三个没有椰蓉的贫瘠面包块吧？

Evil：没有椰蓉的椰蓉面包，确实有点可怜呢……

八六：所以，当一个块里只剩一个椰蓉的时候，就可以停下四分的过程啦。

Evil：一直要分到不存在有一块里面有多个椰蓉的时候吗？

八六：没错！evil同学，一下就理解了呢~

Evil：虽然是这样……但是这样分完之后，又该怎么像八六之前说的那样一个椰蓉只用对远处的一大团椰蓉计算引力就好了呢？现在这样分完，还不是要每一个带有一个椰蓉的小块之间成对计算引力吗？

八六：这时候就是将图形结构化的好时候啦~等我一下。

（八六又在图上画了一些东西）

Evil：啊，仔细一看，确实可以将之前的平面图变成这种树状图呢。

八六：没错！对于我们的椰蓉面包表面宇宙是四分，也就是四叉树。同理，对于我们所在的三维空间，也可以将一个大立方体分成八个小立方体，从而获得对应的八叉树哦。

Evil：很巧妙的思路转换呢~

八六：看着这个树形结构，evil的话，应该已经看出来哪些椰蓉之间需要怎么计算引力了吧~

Evil：嗯……同级的，比如B和C之间和EFG之间就需要老方法两两成对计算吧……那么对于B对D，C对D这种不在同一级的计算……

八六：evil同学加油！

Evil：啊！就可以把B和C的质量加在一起，然后放到它们上级区域的中心，当作上级的那个母节点的等效椰蓉团吧！这样B-C母节点就是一个和D同级的椰蓉了，也就可以之间用引力定律两两成对计算！

八六：evil果然很聪明啊~不过补充一点，与其用BC质量之和放到上级区域中心，也可以用BC质量与位置的加权求和得到更准确的等效引力中心哦，这在不是四等分而是方块有大有小的时候也很好用~

Evil：所以计算一个节点的引力，就可以递归的从子叶一直合并到同级的节点再计算引力呢，比如BCD对A的引力。然后再根据牛顿第三定律也可以得到A对BCD的引力。这样的话……确实可以省去不少计算量呢！

八六：啊，其实并不能用牛顿第三定律得到A对BCD的引力哦？因为只能得到A对同层的东北节点，或者说是A对BC合并和D合并的质心的引力嘛。

Evil：啊，对哦……

八六：所以说，要是计算的对象在叶子节点里面，比如对于B来说就有A和D，这种还是要直接用万有引力定律计算呢。不过还是会比最开始的方法快不少吧？特别是有很多稀疏的椰蓉可以聚成很多小团的时候。

Evil：是呢……不过要是有一团比较紧密的椰蓉和一个离得很远的椰蓉，那这一团的椰蓉对这个离得很远的椰蓉的引力也要一对一对计算呢，明明感觉可以当成差不多的力……

八六：是吧？不过，一般来说，相比起原来算法的$O(N^2)$的时间复杂度，Barnes-Hut 算法一般能有$O(N\log N)$的表现呢，虽然最坏的情况就像刚刚的例子一样，还是会退化到$O(N^2)$。

Evil：已经很棒了呢！大概。

八六：当然，还有更棒的方法哦？甚至能做到$O(N)$的时间复杂度呢。

Evil：$O(N)$？！也就是快速多极子方法吧？快说啦，八六同学！

八六：啊，快要到下午上课时间了，还是要好好上课才行哦evil同学~

Evil：诶，怎么这样……那下午放学后，记得等我哦！

八六：（目视着跑出自习室的evil）呵呵，就是认真好学这点才让人感觉到十分可爱啊，evil同学~

快速多极子算法

---

（下课铃响，八六正收拾着桌上的数学，计算机和物理课本，抬头看到了Evil正站在教室门口）

八六：啊，这不是evil同学吗？有何贵干？

Evil：骗人的吧……这么快就忘了中午和我的约定吗八六酱……

八六：呵呵，evil同学，想稍微捉弄一下呢。

Evil：啊啊，果然还是回去问AI吧，八六同学太坏了。

八六：诶……对不起啦，evil同学，看来做过头了呢……

Evil：好啦，要道歉的话，就快点讲清楚吧，快速多极子什么的。

八六：呵呵，那教室的黑板，就借用一下啦~

Evil：什么吗，这种坏坏的笑容，完全不像是有好好反省自己的样子……

（八六整理了下裙摆，走上了讲台，而evil很自然的坐在了八六的座位上）

八六：那么，我们开始吧。

（八六拿起粉笔在黑板上画着，evil在下面认真的看着）

八六：（在八六还在画图的时候）首先，我们先解释一下这个方法的名字吧~

八六：快速的部分很好理解，毕竟实际上确实在N比较大的时候，FMM是比B-H算法快的~而多极子……evil同学知道偶极子和四极子之类的吗，比如电偶极子的概念。

Evil：知道哦，电偶极子不就是两个等电量的异性点电荷构成的结构吗。好像在距离远超于两个点电荷相隔距离之处它们所产生的电场可以近似为它们的一个叫电偶极矩的量所产生的电场。要是它们的相隔距离趋向于 0 ，好像就是点电偶极子了，也是最纯粹的电偶极子，好像性质只由电偶极矩决定呢。

八六：诶？！（震惊地看向evil）这些可是我们课本上还没学到的知识呢！没想到evil酱居然已经知道了，我还想好好炫耀一下呢，什么的。

Evil：什么嘛，看不起我的样子……

八六：对不起！太得意忘形了……

Evil：嘛，好歹也是稍微看过一点费恩曼物理学讲义的，虽然看第三卷的时候看到一半就看不下去了……

八六：evil同学，意外的很强呢！

Evil：把意外去掉啊！

八六：嘛，不过快速多极子的多极子，其实是指多极子展开（Multipole Expansion，ME）哦。

Evil：啊，我猜猜，多极子展开就是把各种各样的电荷结构变成很多种多极子吧。

八六：没错！就像傅里叶展开可以把任意函数展开成用不同频率的正弦函数的级数一样，多极子展开可以把任意电荷分布展开成用不同极数的多极子的级数呢。而刚刚evil说的电偶极子，就像是正弦函数一样是纯粹的呢，在多极子展开里就像正弦函数在傅里叶展开里一样都只有一个项哦。

Evil：哼哼，我真厉害！那……难道就像函数近似一样，因为越高次的极子对远场的影响的距离上的次方越小，比如单极也就是一个点电荷的电势正比于$\frac{1}{r}$，双极正比于$\frac{1}{r^2}$，四极正比于$\frac{1}{r^3}$等等，也可以按精度在级数的某一级截断吧。

八六：没错，而且误差项也是很明确定义的！

八六：总的来说，ME就是将一组椰蓉对远处观测点的贡献表示为关于这组椰蓉源中心的级数展开的操作啦。比喻一下的话，就是可以像“发射器一样”向自己周围广播自己对外面椰蓉的影响吧~

八六：嗯……还有一点想补充说明的。evil同学知道的吧，在某些力的计算里面，力是矢量所以计算起来不仅要考虑大小还要考虑方向很麻烦呢，而势是标量可以直接简单相加呢。到时候可以再计算势的梯度来求得力。就像evil之前说的：

Evil：是呢……不过要是有一团比较紧密的椰蓉和一个离得很远的椰蓉，那这一团的椰蓉对这个离得很远的椰蓉的引力也要一对一对计算呢，明明感觉可以当成差不多的力……跳转

八六：要是是势的话就可以很简单的直接让远场的各个椰蓉团的ME的势加和在一起求出整体的势了吧？嘛，虽然按照同样的思想，也可以先求出整体的力场之类的……这种细节就不要在意了吧？

Evil：很在意啊！

八六：最后总结的时候会讲到的哦所以就先忍耐一下吧~

八六：好啦，终于画完啦，有彩色粉笔真是不错呢。

Evil：标记起来很方便呢~

八六：接下来稍微说一下FMM到现在的流程吧~第一遍的ME将所有在叶节点的小团，也就是黑板上绿色的这几个框框啦，里面的椰蓉的势函数级数求出来。然后，和Barnes-Hut类似的，虽然数学方面更加复杂一点，可以把这几个绿色的小团的势函数平移和合并，变成红色框框的势函数呢。最后一直向上传递势函数，直到获得整个模拟里最高层的框框的势函数为止呢。顺便一提，这两个步骤分别叫做粒子到多极子（Particle to Multipole，P2M）和多极子到多极子（Multipole to Multipole，M2M）呢，而这两个步骤合在一起就是FMM的向上传播的过程了。

Evil：那个彩色的框框是什么呢？

八六：啊，那个只是无聊才试着画画渐变啦，想表达的是要是是引力势的话热图看起来大概长什么样哦。

Evil：呐鲁侯多。

八六：经过P2M和M2M之后，所有的节点，也就是所有的框框，都计算和储存了一份自己的ME展开的级数的系数呢，就像$M_1,M_2,…$。这个向上传播的过程是对所有框框进行的，而框框的数量和总共椰蓉的数量是成线性关系的呢，所以这一步的时间复杂度是$O(N)$哦。

Evil：（总感觉怪怪的……真的是线性吗？）

Evil：不过这样的操作之后，接下来就是用总的ME展开的场对所有的椰蓉进行计算了吧？就像已经有了势场，只需要对每个椰蓉找到势的梯度就解决了！那确实也是$O(N)$的一次pass就解决了呢~

八六：不要这么快嘛，八六同学~刚刚的讨论里面，我们都忽略了多极子展开得到的函数的适用域哦？就像泰勒级数一样，多极子展开得到的多项式有可能会有收敛范围呢，毕竟一开始就是假设源椰蓉对远场的作用才想到多极子展开的嘛~所以，要保险起见的话，自己同级的框框肯定是要直接里面的椰蓉成对计算了吧，甚至上一级和上两级的也要呢~

Evil：呜……那这么说的话还是感觉没什么用啊，引入多极子展开这种事。虽然可以对一小团椰蓉求出它们的势再用在远处的椰蓉上求力，但是感觉和Barnes-Hut算法的想法大差不差呢，还不是需要分级运算之后各个子叶成对计算，只是把最末尾的一级从单个椰蓉变成了单群靠的很近的椰蓉，在八六画的图上就是还不是要进行橙色框框一对一对的对红色框框进行计算，然后之后这些橙色框框又要对绿色框框进行同样重复的计算……在大尺度的模拟下也只是对总运算量除以了一个常数，比如每格最少要10个椰蓉以上才再细分，考虑到椰蓉不是完全均匀分布，就只是对总运算量除以了不到10吧，还是没能想到为什么会比Barnes-Hut有质的变化呢。

八六：确实呢，要是单看ME的话。不过计算误差有了明确的定义，还是有进步的呢。

Evil：唉……结果多极子展开也没能解决问题吗？要是这些计算像向上传播的过程一样，每层都能存下来一个信息给自己的下一层用就好了，不过下面层的节点又应该怎么……

八六：诶？！evil同学难道会读心吗？快速多极子方法确实还不止有ME这个操作，还有另一个操作，叫做局部展开（Local Expansion，LE）呢。

Evil：诶，那又是什么？

八六：是一种将远处的椰蓉团对某区域内部所有椰蓉的贡献表示为关于该区域中心的级数展开的操作哦。也就是说，能像“接收器”一样把各个框框的ME“发射器”对自己格的影响转换成以自己为中心的表现形式呢。

Evil：不就是换了个坐标原点吗……像是把泰勒展开从$(x-x_1)$中心变成$(x-x_2)$一样，感觉没什么用呢……

八六：呵呵……单看一次LE操作确实会有这种感觉呢，拿黑板上的图来说，就只是把橙色框框为中心的展开的原点移到了红色框框里嘛。但是，接下来的一步就是连接FMM向上传递和之后向下传递的桥梁，也就是多极子到局部（Multipole to Local，M2L）操作。

Evil：啊，我更在意向下传递呢，跳过吧~（微笑）

八六：啊啊！evil同学也会使坏了！

Evil：我本来就是evil啊……嘛，neuro完全是另外一个层级的evil就是了……八六同学，继续M2L吧~

八六：好，好吧……还不知道evil的姐姐是什么样的人呢……

八六：M2L，就是对与它“足够远”、满足多极展开收敛条件的格子，将其他同级的格子的ME展开利用LE转换成以自己为原点的展开形式，之后再汇总起来，记录成LE的系数，比如$L_1,L_2,…$，就是为了接下来的向下传播把同级的ME信息都收集起来了呢。所以这也是FMM很核心的一步呢，就是evil刚刚想到的：

Evil：唉……结果多极子展开也没能解决问题吗？要是这些计算像向上传播的过程一样，每层都能存下来一个信息给自己的下一层用就好了，不过下面层的节点又应该怎么……跳转

Evil：哇，难道，我真的是天才？

八六：当然是哦，evil~读到这里还一遍都没有倒回去的读者桑，也是天才呢~

Evil：读者桑？啊，就是刚刚一直感觉到的视线吗？多摩多摩~

八六：接下来，就是向下传播的过程呢。用刚刚得到的LE的系数，再利用局部到局部（Local to Local，L2L）向下让每个框框都有自己的LE级数系数。

Evil：到这里，每个框框都有了一组ME的系数和一组LE的系数呢。

八六：嗯~最后，只要对框框里的椰蓉使用自己的LE展开的系数还原LE展开的级数，也就得到了势函数了，力也就顺势得到了。这步就是局部到粒子（Local to Particle，L2P）。哦，当然势函数里面还要记得加上近场那些不能用级数近似的椰蓉的影响哦~

Evil：啊，从椰蓉开始转了一圈回到椰蓉了呢，其间就是层层更新了各种信息，对原本两两组对计算引力的纠缠不清的问题建立了结构，而每个操作都是对和椰蓉个数成正比的框数进行遍历，时间复杂度也还是$O(N)$呢，FMM总体就一定是$O(N)$的时间复杂度了~虽然感觉这个$O(N)$前面藏起来的系数和常数项一定大的离谱……

八六：不过，实际使用起来确实会比Barnes-Hut快很多哦？在计算二维引力的时候，好像只需要大约100k以上的椰蓉FMM就能快过B-H了，当然，不同的实现也会有运行速度上的差异啦。

Evil：嗯……虽然感觉讲完了，但是感觉心里还是沉沉的，感觉有好多跳过的部分或者细节不清楚的部分呢……

八六：嘿嘿，其实有部分细节八六我是故意跳过的呢，比如ME展开和LE展开和这些M2M，L2L的上下传递的具体形式。因为！不只是引力，只要是能有“远处的作用没有近处的作用强，可以适当忽略”的性质的操作，就可以用FMM来分层加速计算呢！比如多极子诞生的领域——静电力，还有声学，分子动力学，流体力学，甚至是机器学习和量子化学都能看到FMM的身影呢。

Evil：呜哇，这么一想确实呢。

八六：这些不同的过程的多极子展开和局部展开以及平移操作的形式都完全不同，当然，收敛域也会随之变化，从而近场和远场的定义也会跟着变化呢。（看向evil）看着evil同学皱着的眉头，想必evil也想到了吧，要是有些影响随距离下降很慢的，又或者频率/信息密度很高的场，FMM又该如何应对呢？

八六：这需要涉及到我们对精度的要求了呢。之前说的ME，LE等等操作，实际上都是假设了展开的级数展开阶数$p$是一个定值，才得出了$O(N)$的时间复杂度的结论呢，实际上的时间复杂度在二维是$O(pN)$，三维至少是$O(p^2N)$呢，因为要用球谐展开来替代多极子展开了。要是展开阶数$p$和$N$有函数关系的话，也就是说不是固定展开阶数的时候，FMM的时间复杂度一般会退化成$O(N\log^\epsilon N)$呢，$\epsilon$具体是多少也要看是对什么场进行计算了，说不定会比Barnes-Hut算法还慢呢。当然，要是碰到非常极端的情况，比如所有椰蓉都聚在一团，那都只能考虑成近场处理，就退化成朴素方法的$O(N^2)$了。

Evil：不过，想想椰蓉面包上只有一个椰蓉向上堆成的高塔，也是挺好笑的呢~

八六：呵呵，是啊，那老板也是费心了呢~

八六：还有前面提到的问题呢：

八六：要是是势的话就可以很简单的直接让远场的各个椰蓉团的ME的势加和在一起求出整体的势了吧？嘛，虽然按照同样的思想，也可以先求出整体的力场之类的……这种细节就不要在意了吧？

Evil：很在意啊！跳转

八六：其实是可以的，但是用矢量场的ME，LE等等操作的展开会很麻烦，也会增加空间复杂度呢。但是对非保守力，比如有摩擦力等等，因为可以经过一个回到原点的位移之后能量却不一样，所以定义不了势函数这么便利的东西，也只能硬着头皮用力场这种矢量场当展开的对象了。

八六：还有Evil之前想的：

八六：经过P2M和M2M之后，所有的节点，也就是所有的框框，都计算和储存了一份自己的ME展开的级数的系数呢，就像$M_1,M_2,…$。这个向上传播的过程是对所有框框进行的，而框框的数量和总共椰蓉的数量是成线性关系的呢，所以这一步的时间复杂度是$O(N)$哦。

Evil：（总感觉怪怪的……真的是线性吗？）跳转

Evil：不是，八六同学怎么会知道我在想什么啊！

八六：毕竟evil同学当时的表情就是这种表情嘛！

Evil：我真的有这么好懂吗……

八六：关于FMM本身在“力的计算”这部分是线性时间复杂度的，除了刚刚对展开级数$p$的影响之外，基本是没有异议的。但是关于建立四叉/八叉树，也就是框框的数量和总共椰蓉的数量是不是线性操作呢？也是比较复杂的问题呢，不过简单来说，在工程上我们会认为是线性的。因为，其实，evil酱知道吗，所有的算法都不可能是线性的~

Evil：诶？诶？！

八六：因为随着$N$的增大，计算机里面要用更多bit来表示$N$或者$N$的地址呢，那么只要涉及到要存取数的情况，就需要考虑这个超级小，但是确实和$\log N$正相关的变量，不是吗~

Evil：确，确实……但是也太钻牛角尖了感觉……

八六：嗯，所以涉及到现实的，其实很多东西都有很多隐性的前提呢~

Evil：看费恩曼物理学讲义和其他的物理知识的时候深有体会啊……即使广义相对论说是只需要等效原理和广义相对性原理，推理前还需要相信物理定律的广义协变性，实际应用的时候还需要考虑能量条件等等呢……

八六：不过，为了解决这些问题，有很多种改进版本的FMM呢，比如自适应FMM让网格划分更高效，还有Kernel-independent FMM (KIFMM)，就不需要考虑运算的核是引力还是什么其他的了，适用于任意平滑核函数。还有……

八六：（肚子咕咕叫）（脸红）

Evil：啊，才注意到时间，怎么感觉外面都快晚上了！对不起啊八六酱！

八六：嗯，没关系哦，倒不如说很开心呢！之后的一些琐碎的内容，就用邮件发给你吧~

Evil：好啊。让我来带八六同学吃点除了椰蓉面包的其他的吧！章鱼小丸子怎么样？

八六：想要物理口味的呢~

Evil：好~

参考资料

---

相关链接：

介绍FMM的文章：

https://zhuanlan.zhihu.com/p/2061732285787936506

https://blog.csdn.net/weixin_42849849/article/details/153407237

https://www.bohrium.com/sciencepedia/feynman/keyword/multilevel_fast_multipole_algorithm_(mlfma)

实现FMM的github库：

https://github.com/davidson16807/fast-multipole-method

介绍和实现FMM的视频：

https://github.com/keyframe41/Videos/tree/main/Fast-Multipole-Method

（顺便这是作者Keyframe Codes参加#SoME5的视频呢）

后日谈

---

八六发给Evil的邮件：

亲爱的Evil同学：

近来可好？相信你还是想着今天下午我说的“琐碎的内容”吧~这就把讨论用附件发给你哦。

顺便，想让Evil同学下次给我讲些物理上好玩的问题呢~

祝：

学业进步，更爱八六~

附件：

heartheartheart♥