本文是快速多极子算法FMM入门的续集。在前一篇文章中，我们已经了解了如何使用FMM解决一个一维的简单求和问题。在本文中，我们将推广这一思想，用以解决更广泛的高维问题。虽然我们仍将围绕一个具体的示例问题来讨论，但与前一篇文章最大的区别是：本文中的叙述将适用于相当广泛的一大类问题，而非针对示例问题中的独特性质来设计。
问题的背景与其它解法
这个示例问题仍然是OI背景的；具体来说，是我从2022年的集训队论文中看到的，来自邓明扬命制的模拟赛题目《差不多得了》。题面如下：
给定 $n$个二维平面上的整点，对每个点求到所有其他点的距离之和，允许相对误差范围 $\epsilon=10^{-4}$，保证 $n\le 2\times 10^5$，坐标绝对值不超过 $10^9$。
我的解法
首先讲一下我第一次看到这道题时自己的思考吧。我在2022年看过一个YouTube视频：A Ridiculous Approximation，视频的主题是如下的近似公式：
\[\sqrt{x^2+y^2}\approx 0.96x+0.4y,\qquad x\ge y\ge 0.\]第一次看到这个公式时，我相当惊讶；然而冷静思考之后我意识到这并没有什么神秘之处，因为这就是在用一条直线近似圆弧；想象一下等式两边函数各自的等高线就明白了。视频中的近似公式就是在长为$\frac{\pi}{4}$的圆弧上的近似，而一般来讲对于长度$x\to 0$的圆弧，对应的弦到原点的距离为$\cos \frac x2$，我们可以选择位于这条弦和与其平行的切线正中的直线作为近似。这样的直线近似带来的误差是
\[\frac12(1-\cos \frac x2)=\frac{x^2}{16}+O(x^3).\]代入$x=\frac{\pi}{4}$，得到误差为3.8%，确实与视频中描述的相同。回到例题中，我们用这个近似公式就很容易得出一个做法了：把圆周等分为$N$份，在每一个扇形中使用上述的近似公式将距离线性化来求和。具体来讲，把扇形的两边使用线性变换变到$x$轴和$y$轴上，然后简单地沿$x$排序后扫描线就可求出线性函数的和。由于误差是$O(x^2)$的，需要选择$N=O(\epsilon^{-1/2})$达到精度要求，故复杂度为$O(n\epsilon^{-1/2}\log n)$。实际上由于误差是$\frac{x^2}{16}$，应该跑起来还算快。
标算
标算的简洁和优美给我留下了深刻的印象。二维中长度为1的线段在随机方向上投影长度的期望为$\frac{2}{\pi}$，因此我们取$N$个角度的投影，将其转换为一维问题，排序后简单统计出每个点到其它点的距离和；最后输出所有投影方向上答案的平均值乘上$\frac{\pi}{2}$即可。由于用$N$个等距节点近似一个光滑函数的积分误差为$O(N^{-2})$，该算法也是选择$N=O(\epsilon^{-1/2})$，故时间复杂度也是$O(n\epsilon^{-1/2}\log n)$，和我自己想出的算法相同，但是实现要简单太多了：见dist_sum_2d_by_projection.cpp。
其实这一算法的时间复杂度仍可进一步减小：我们并不需要精确地求解一维问题，只需要让其相对误差不超过$\epsilon$即可。因此我们将每个点投影后的坐标四舍五入近似到$O(\epsilon^{-1})$个离散值之中，就能用桶排序在$O(n)$的时间内完成排序，将整个算法的时间复杂度降低为$O(n\epsilon^{-1/2})$。这一思想在近似算法中相当常用。
FMM
虽然通过投影期望计算距离和的方法相当优美，但它本质上甚至不是一个多项式时间算法，在$\epsilon$减小时（如$10^{-10}$）将会完全失效。并且距离和问题的一维类比非常容易计算，所以才能使用投影来降维，然而更广泛的问题并不具备这一性质。下面我将展示FMM是如何处理普遍的高维求和问题的，从而将距离和问题作为一个特例来解决。这一节所有的东西都是我独立想出来的，所以可能和常见的FMM算法并不完全相同，但思想是同样的。
问题定义
设$x_1\dots x_n$是$\mathbb{R}^d$中的$n$个点，分别带有权重$w_1\dots w_n$，而$f$是一个$\mathbb{R}^d$中除原点之外解析的函数，我们需要计算如下的求和：
\[g(x_i)=\sum_{j\ne i}f(x_i-x_j)w_j.\]并且假设$f$的解析式已知，从而可以任意地求导。我们需要让计算结果的相对误差达到某个给定的$\epsilon>0$。对于这个距离和问题，有$d=2,f(x)=|x|,w_j=1$。对于前一篇文章中的问题，则是$d=1,f(x)=\frac{\operatorname{sgn}(x)}{x^2},w_j=q_j$。
前置知识：多元泰勒展开
类似一维的情形，对于一个$\mathbb{R}^d$中的解析函数$f$，有如下的展开：
\[f(x)=\sum_{\mathbf{k}}\frac{1}{\mathbf{k}!}\partial_x^\mathbf{k}f(x_0)(x-x_0)^{\mathbf{k}},\]其中$\mathbf{k}=(k_1,k_2,\dots,k_d)$是一个自然数序列，记$\mathbf{k}!=k_1!k_2!\dots k_d!$, $\partial_x^\mathbf{k}=\frac{\partial^{k_1+\dots+k_d}}{\partial x_1^{k_1}\dots \partial x_d^{k_d}}$, $x^\mathbf{k}=x_1^{k_1}x_2^{k_2}\dots x_d^{k_d}$, $|\mathbf{k}|=k_1+k_2+\dots+k_d$. $P-1$阶泰勒展开给出的近似就是仅保留所有$|\mathbf{k}|<P$的项给出的部分和。
在FMM中，我们的卷积核$f$常常来自Green函数，因此在对角线部分（即$x_i=x_j$处）不可导，而在非对角线部分是光滑的。虽然距离和问题中的$|x_i-x_j|$不是一个Green函数，但它也满足类似的性质。根据一维的情形类比，一个在原点处不可导、其余地方解析的函数在$x$点的泰勒级数收敛半径是$|x|$，故在以$x$为圆心、半径$\frac{|x|}{2}$的圆内泰勒展开是线性收敛的，取$P=\log_2\epsilon^{-1}$即可。（这一点我其实没证明过，不敢保证是对的，但感性理解觉得还是挺对的，实际上也没出问题）
FMM的核心公式
在前文中，我们已经知道了FMM所需统计的两项数据：$M$和$L$。其中$M_p$表示的是节点$p$内部的求和结果，$L_{p’}$表示的是节点$p’$处的远场求和结果，这里的“远场”是指除了节点本身和相邻节点之外的全部位置。我们仍然用$x$表示求和的位置，用$y$表示产生贡献的点。在节点$p$的中心$c_p$处我们展开$f(x-y)$得到
\[\begin{aligned} f(x-y)&=\sum_{\mathbf{k}}\frac{1}{\mathbf{k}!}\partial_y^\mathbf{k}(f(x-y))\vert_{y=c_p}(y-c_p)^{\mathbf{k}}\\ &=\sum_{\mathbf{k}}\frac{1}{\mathbf{k}!}\partial_x^\mathbf{k}f(x-c_p)(-1)^{|\mathbf{k}|}(y-c_p)^{\mathbf{k}}. \end{aligned}\]因此对于足够远的$x$（即$|x-c_p|$超过节点$p$代表区域的半径的两倍），在精度范围内就有
\[\sum_{j\in p}f(x-x_j)w_j=\sum_{|\mathbf{k}|<P}\frac{(-1)^{|\mathbf{k}|}}{\mathbf{k}!}\partial_x^\mathbf{k}f(x-c_p)\sum_{j\in p} w_j(x_j-c_p)^{\mathbf{k}}.\label{def_m}\tag{1}\]因此我们维护$M_{p,\mathbf{k}}=\sum_{j\in p} w_j(x_j-c_p)^{\mathbf{k}}$. 在叶节点处暴力求和算出$M$，而M2M过程就是一个多项式平移，从子节点的$c’$平移到父节点的$c$。这可以用$d$步完成，第$i$步做$P^{d-1}$个一维多项式在第$i$维度上的平移。
下面我们转换关注点，考虑用$p’$节点的中心$c_{p’}$上的求和结果表示$x$上的求和结果。为了做到这一点，在$c_{p’}$点展开$\eqref{def_m}$式即可：
\[\begin{aligned} \sum_{j\in p}f(x-x_j)w_j&=\sum_{|\mathbf{k}|<P}\frac{(-1)^{|\mathbf{k}|}}{\mathbf{k}!}(\sum_{|\mathbf{k}'|<P}\frac{1}{\mathbf{k'}!}\partial_{x}^{\mathbf{k}'}(\partial_x^\mathbf{k}f(x-c_p))\vert_{x=c_{p'}}(x-c_{p'})^{\mathbf{k}'})M_{p,\mathbf{k}}\\ &=\sum_{|\mathbf{k}'|<P}\left(\frac{1}{\mathbf{k'}!}\sum_{|\mathbf{k}|<P}\frac{(-1)^{|\mathbf{k}|}}{\mathbf{k}!}\partial_x^\mathbf{k+\mathbf{k}'}f(c_{p'}-c_p)M_{p,\mathbf{k}}\right)(x-c_{p'})^{\mathbf{k}'}. \end{aligned}\]由此可见，我们需要维护
\[L_{p',\mathbf{k}'}=\sum_{p\textrm{ is far from }p'}\frac{1}{\mathbf{k'}!}\sum_{|\mathbf{k}|<P}\frac{(-1)^{|\mathbf{k}|}}{\mathbf{k}!}\partial_x^\mathbf{k+\mathbf{k}'}f(c_{p'}-c_p)M_{p,\mathbf{k}}.\]该式即代表了M2L过程；远场对于位于节点$p’$中的$x$点的贡献即为$\sum_{|\mathbf{k}’|<P}L_{p’,\mathbf{k}’}(x-c_{p’})^{\mathbf{k}’}$。我们还需要一个L2L过程，把父节点的$L$传给子节点。这可以用$d$步完成，每一步把中心$c$的某一维坐标从父节点调整到子节点的位置。具体来讲，假设要把$c$的第$i$维调整到$c’$的第$i$维，则有
因此给子节点的$L$的$j$位置加上父节点的$L$的$k_i$位置乘上$\binom{k_i}{j}(c’_i-c_i)^{k_i-j}$即可。还有一种思考方式就是把$\sum_{|\mathbf{k}’|<P}L_{p’,\mathbf{k}’}(x-c_{p’})^{\mathbf{k}’}$直接在子节点的$c’$处泰勒展开，可以得出等价的L2L公式，就是推导稍微麻烦一点，这里不展示了。
程序实现细节
由于FMM要求每一个节点上的级数在该节点所代表的区间内足够快速地收敛，我们必须让节点都是正方形才能得到最高的效率，因此FMM一般使用的数据结构是四叉树、八叉树等。与一般的仅有沿树边的访问需求的数据结构不同，FMM里还需要获取一个节点相邻的全部节点，因此在程序中我采用的策略是把每一层的节点组织成一个$d$维数组，比如二维的情况就是使用$(dep,p_x,p_y)$来代表一个节点；其父节点位于$(dep-1,p_x/2,p_y/2)$，邻居们位于$(dep,p_x+\{-1,0,1\},p_y+\{-1,0,1\})$，四个孩子位于$(dep+1,2p_x+\{0,1\},2p_y+\{0,1\})$。在叶节点上则使用一个std::vector来存储其中的全部$x_i$。
为了方便起见，我假设输入的$x_i$都是均匀分布的，建立固定高度的满四叉树，这样所有叶节点就能位于同一深度了；否则代码会相当复杂，不能采用上文所述的简单存储方式了，查找每个方向上的邻居节点不再是平凡的事情；M2L时也需要特判父亲的邻居是叶节点的情况。但就算是非均匀分布的，下文的复杂度分析也不会有任何影响，就普通地动态开点建树，然后在递归时维护一下邻居的情况，即可用同样的复杂度完成均匀分布下能完成的所有事，并不是我不会写，只是我很懒。
我们可以选择恰当的层数，使得叶节点中的$x_i$个数为$O(P^d)$。之前描述的M2M和L2L过程由于是每个维度分开卷积的，时间复杂度都是$O(P^{d+1})$；但M2L一般来讲是不可分开的卷积，故复杂度为$O(P^{2d})$。每个叶节点上的暴力都需要$O(P^{2d})$的计算量，因此总复杂度是$O(nP^{-d}P^{2d})=O(nP^d)=O(n\log^d\epsilon^{-1})$。与一维的情形类似，我们不可能用其他的卷积策略达到更优的复杂度，因为对于$n$个点中的每一个都需要$O(P^d)$的计算量来给泰勒级数求值。并且一般的问题规模中$P$非常小，暴力卷积总是最快的卷积策略。
值得注意的是，这一复杂度与普通的物理问题FMM的二维$O(nP)$，三维$O(nP^{3/2})$无法比较。因为物理问题中的核总是调和函数，在一个区间内的值仅由区间的边界决定，所以自由度是和区间边界的测度成正比的；而我们所解决的是所谓的kernel-independent FMM问题，对于一般的核（如本问题中的$|x-y|$）只在区间的边缘上设置代理点是无论设多少个都不够的。此外，虽然所有的FMM文献都号称自己在$P=O(1)$时是$O(n)$的复杂度，但这是不可能做到的，光是建成四叉树/八叉树就至少要$O(n\log n)$的代价了。
在M2L过程中我们需要计算$f$的任意阶偏导。对于本问题中的$f(x)=|x|$，归纳易得
其中的系数$coeff_{k_1,k_2,k}$可以简单地递推算出。我们在M2L中总共需要查询$O(P^2)$个偏导，因此用$O(P^3)$的时间暴力计算就可以了；但就算偏导数计算代价很大也没有关系，在FMM过程中只要保持均匀地分割空间，两个邻居节点间本质不同的位置关系就只有$O(1)$种，预处理出每种情况的各阶偏导即可。
完整的程序见fmm_2d.cpp。
变种：使用其它的基
这里使用了多元泰勒展开来近似核函数$f$在远离原点处的值。就像多项式有系数和点值两种表示，除了多元泰勒展开之外，我们也可以使用$f$在一系列点的取值来近似$f$。根据数值分析，如果我们选择Chebyshev节点上插值出的多项式来近似$f$，在$f$的各阶导数有界的前提下误差会随着节点数线性收敛到0，因此也是只用$O(P^d)$个点就能达到所需的精度。（这是我最喜欢的数值方法之一，我的本科毕业设计也与此密切相关）对于点值表示，我们也能计算出相应的M2M、M2L和L2L系数矩阵，当然比这里的多元泰勒展开给出的公式更加复杂，但整个算法的流程还是完全相同的。相比于多元泰勒展开，点值表示的主要优势是无需计算高阶导，因为并非所有的$f$都是知道解析式的，而且就算知道解析式求高阶导可能也相当复杂。
对于物理问题里一些特殊的核，科研中会针对这些核的特性设置一些独特的展开（如将$\mathbb{R}^2$看成$\mathbb{C}$用复数泰勒级数、$\mathbb{R}^3$中的球谐展开），使得用较低的展开阶数就能达到足够精度，且转移矩阵还存在简单的封闭形式，这里就不展开讨论了。
理论时间复杂度
从理论计算机科学的角度来看，其实理论上我们可以得到更优的复杂度。首先，数值计算里是没有“多项式多点求值”这种科技的，虽然确实存在比$O(P^2)$更快的算法，但要么是需要用高精度来运算（source）、要么是需要用Chinese Remainder / rational reconstruction等技巧转化为有限域运算（source），具体的时间复杂度无从得知。这里我们索性做出最好的假设，即$T$个点上的$T$次多项式求值要$O(T\log^2 T)$的时间，前一篇文章所涉及到的一维FMM问题理论复杂度就将变为$O((n/T)(P\log P+P\log ^2 T+T^2))$，这个$T$的最小值我不会求，但易得$T=\Omega(\sqrt{P\log P}),T=o(\sqrt{P}\log P)$。总之，理论上讲一维FMM的复杂度应该是$\tilde{O}(n\sqrt{P})$的。类似地，在假设多项式多点求值真的有这么完美的做法的前提下二维的FMM是$\tilde{O}(nP)$的。
与一维情形不同的是，在二维同时计算$P$个点的泰勒级数求和可以做到$O(P^\omega)$，因为本质上是在同时计算$P$个$u^T Av$，拼成一个大矩阵去求$U^TAV$的对角线即可做到矩阵乘法的复杂度。由于$M$和$L$之间的那些转换过程都是卷积，而偏导数可以预处理，单个节点的计算时间理论上是$O(P^2\log P)$的，故总所需时间为$O((n/T)(P^2\log P+TP^{\omega-1}+T^2))$，易得最小值的$T$为$P\sqrt{\log P}$，故时间复杂度为$O(nP^{\omega-1})$，这是无需任何假设的二维FMM理论时间复杂度。
还有一种神奇的方法是抛开泰勒展开/Chebyshev插值等有固定过程的近似算法不管，直接使用低秩近似作为$M$和$L$。与前面那几种galactic algorithm不同，这确实是一种现实中可行而且广泛应用的FMM算法，美中不足的是理论复杂度难以分析：低秩近似的rank到底是多少？我画出了精度要求和所需rank的关系图，结果显示对于$|x-y|$这个核，rank随着$\log(1/\epsilon)$的增长快于线性，但也不敢说就是泰勒展开/Chebyshev插值给出的平方级。我还观察到，无论是泰勒展开还是Chebyshev插值，达到精度所需的展开项数都超过数值rank，但也没有很大的差距。在这种差距下，SVD所带来的巨大常数就不可忽视了，因此我想低秩近似可能无法快过解析求解的FMM。
结果
我首先比较了原题目的相对误差要求$10^{-4}$下$O(n\epsilon^{-1/2}\log n)$的标算和$O(n\log^2\epsilon^{-1})$的FMM的用时。标算需要选择$N=13$个投影，用时为0.240s；而FMM需要展开$P=5$阶，四叉树深度为$6$，用时0.168s。由此可见，就算精度要求极低时，FMM也要比标算更快。而如果提升精度要求到$10^{-6}$，标算所需的投影数就增长到$65$个，用时0.786s；而FMM仅需增加展开到$11$阶，四叉树深度降到$5$，用时0.305s。精度要求$10^{-10}$时，标算需要多达七千个投影，跑了42.9s，已经比暴力还要慢4倍了；而FMM需要$29$阶展开，四叉树深度降到$4$，用时仅0.834s。这与时间复杂度给出的预测一致，精度要求更高时FMM的优势将显著扩大。在叶节点大小选择恰当时，FMM的运行时间中M2L部分占了总时间的约$1/3$，而叶节点上的求值和近场贡献统计部分占了总时间的约$1/2$。
值得注意的是，如果对于$10^{9}$范围的输入数据直接运行$P=16$及以上的FMM，会超出double类型所能表示的范围，让程序失效。解决方案很简单，先把输入数据归一化到$[-1,1]^2$里再求解，把答案乘上$10^9$即可。
虽然速度上被FMM完爆，而且无法处理更广泛的问题，这个投影长度求期望的算法也有自己的长处，其中最明显的优点当然是实现比FMM要简单几十倍。对于足够高维的数据（如100维），FMM将完全失效变成指数级的复杂度，而投影算法仍然会比暴力更加有效。对于噪声相当大的数据（假设每次测量点的位置时都会有一个独立的巨大噪声干扰），FMM也是完全失效，而只要有足够多的投影，哪怕信噪比是$1:100$，投影算法也能有效计算。
