空间层次压缩远场
QuadTree/OctTree 让近场精确、远场聚合，把直接 (O(N^2)) 计算降到 (O(N\log N)) 或 (O(N))。
CS267 · 第 26 课
N-body 问题把“每个粒子与所有其他粒子相互作用”的 (O(N^2)) 计算，转化为一组可控误差的层次化近似。本课从 quad/octree 出发，推导 Barnes–Hut 与 Fast Multipole Method，讨论并行树划分、locally essential tree、性能调优，并用 1.5D 复制解释如何用额外内存换取更少通信。
直接 N-body 对每一对粒子计算引力、静电力或其他成对作用。若有 (N) 个粒子，交互数量为 (N(N-1)/2)，计算复杂度是 (O(N^2))。它出现在天体物理、分子动力学、等离子体模拟、边界元方法、计算机图形学和电子束模拟中。
for i in range(N):
force[i] = 0
for j in range(N):
if i != j:
force[i] += pair_force(position[i], position[j])近处粒子必须精确区分；远处的一簇粒子，从当前观察点看几乎位于同一方向，可以用总质量、质心或更高阶矩来近似。层次化方法就是不断回答：“这簇粒子离我足够远吗？”
二维空间递归分成 4 个象限得到 QuadTree；三维空间递归分成 8 个子立方体得到 OctTree。若粒子分布不均匀，应只细分粒子密集的区域，形成自适应树。
逐粒子插入的简单构建成本通常为 (O(N\log N))。若粒子近似均匀，也可以通过排序、空间编码或批量构建把常数显著降低。
| 参数 | 太小 | 太大 |
|---|---|---|
| 每叶最大粒子数 (q) | 树更深、节点更多、遍历开销大 | 叶内直接相互作用变多 |
| 接受阈值 ( heta) | 精度高但访问更多 box | 速度快但近似误差增大 |
| 展开阶数 (r) | 单次展开便宜但精度低 | 精度高但系数和转换成本增加 |
Barnes–Hut 为每个 box 保存总质量 (M) 和质心 (c)。设 box 尺寸为 (D)，目标粒子到 box 中心的距离为 (r)。当：
\[\frac{D}{r} < \theta\]
就把整个 box 当作位于质心 (c) 的单个“超级粒子”；否则继续访问子节点。典型 ( heta) 在 0.5 到 1 附近，但实际值应由误差要求决定。
def force_from_node(particle, node, theta):
if node.is_leaf:
return direct_force(particle, node.particles)
if node.size / distance(particle, node.center_of_mass) < theta:
return monopole_force(particle, node.total_mass, node.center_of_mass)
return sum(force_from_node(particle, child, theta)
for child in node.children)在常见粒子分布上，Barnes–Hut 的工作量约为 (O(N\log N))。极端聚簇或很严格的精度要求会增加访问 box 的数量，因此复杂度和负载都与数据分布有关。
FMM 不只用总质量和质心，而是保存截断的 multipole expansion。它区分两类展开：
二维中可用复数 (z=x+iy) 表示位置。远离源点时，(log(z-z_k)) 可在 (1/z) 上做 Taylor 展开：
\[f(z) = M\log z + \sum_{\ell=1}^{r} a_\ell z^{-\ell} + \text{error}(r)\]
若所有源粒子都在半径 (R) 内，而目标点满足 (|z|>R)，令 (c=R/|z|<1)，截断误差满足：
\[\text{error}(r) = O(c^{r+1})\]
因此每多保留一阶展开，误差按几何级数下降。三维实现通常使用球谐函数或等价展开。
| 算法 | box 信息 | 访问模式 | 典型复杂度 |
|---|---|---|---|
| Direct | 粒子列表 | 所有粒子对 | (O(N^2)) |
| Barnes–Hut | 总质量 + 质心 | 精度越高，访问 box 越多 | (O(N\log N)) |
| FMM | 多极与局部展开 | 每层固定邻域交互 | (O(N)) 或 (O(bN)) |
层次化 N-body 同时包含规则工作与非规则工作：建树、向上聚合和向下传播是树遍历；逐粒子求力的访问路径随粒子位置变化。并行策略必须同时考虑计算负载、树依赖和远程 box 数据。
每个进程不需要复制整棵全局树，只需要“计算本地粒子时可能访问的远程节点”。这部分远程子树称为 locally essential tree。实际实现通常先交换粗层摘要，再按需请求更细节点。
若每个进程把所有树节点广播给所有其他进程，计算虽然减少，通信却会主导性能。高效实现必须利用空间层次、邻域关系和批量消息，限制远程节点集合。
先看不做近似的直接方法。设 (p) 个处理器均分 (n^2) 次相互作用，每个处理器本地快速内存可容纳 (M) 个粒子或中间量。课件给出的下界是：
\[\#\text{words} = \Omega\!\left(\frac{n^2/p}{M}\right)\]
\[\#\text{messages} = \Omega\!\left(\frac{n^2/p}{M^2}\right)\]
它表达了一个普遍规律：更大的本地内存允许重用更多粒子数据，从而同时减少带宽和消息启动次数。
每个 rank 保存 (n/p) 个粒子，并把一份粒子块在环上移动 (p) 次。计算量为 (O(n^2/p))，通信带宽为 (O(n))，消息数为 (O(p))。当 (M=O(n/p)) 时，它达到了上述下界。
将处理器组织为 (c) 个复制层，每个 rank 的内存增至：
\[M = O\!\left(\frac{nc}{p}\right)\]
主循环步数从 (p) 降到 (p/c^2)，总带宽与延迟近似为：
\[\text{bandwidth} = O\!\left(\frac{n}{c}\right), \qquad \text{latency} = O\!\left(\frac{p}{c^2}\right)\]
当 (1 \le c \le \sqrt p) 时可以连续调节内存与通信。极限 (c=\sqrt p) 对应 2D 分解：主循环只剩一步，带宽降到 (O(n/\sqrt p))。
直接 N-body 的两重循环和矩阵乘法的三重循环都可以通过 Hölder–Brascamp–Lieb 线性规划分析数据复用。循环维度、数组下标投影和可用内存共同决定最优 tile 形状与通信下界。
Barnes–Hut 与 FMM 的答案是近似的，因此通信下界会依赖精度、粒子分布、树形状和接受准则。减少算术通常也会减少通信，但“达到某一误差所必需的最少通信量”仍比直接 (N^2) 方法复杂。
课件展示了同一 FMM 实现经过手工优化、自动调优、OpenMP 和 GPU 移植后，性能曲线会显著改变。最重要的参数之一是每个叶 box 的最大粒子数 (q)：它决定树遍历和叶内直接计算的平衡。
| 优化层 | 关注点 |
|---|---|
| 算法 | ( heta)、(r)、(q)、树深、交互列表 |
| 单核 | 数据布局、SIMD、循环展开、减少分支 |
| 多核 | OpenMP task、NUMA first-touch、细粒度与调度开销 |
| GPU | 批量 box 交互、结构化数组、kernel 合并与占用率 |
| 分布式 | LET、消息聚合、重叠通信、拓扑感知划分 |
Q1：Barnes–Hut 使用 (D/r < \theta) 判断什么？
Q2：FMM 相比 Barnes–Hut 的关键增强是什么？
Q3：1.5D 直接 N-body 中增加复制因子 (c) 的主要效果是什么？
Q4：Locally Essential Tree 的目的是什么？
QuadTree/OctTree 让近场精确、远场聚合，把直接 (O(N^2)) 计算降到 (O(N\log N)) 或 (O(N))。
Barnes–Hut 使用低阶聚合并按精度访问不同数量的 box；FMM 使用固定交互列表与高阶展开。
好的树划分不仅要均衡粒子数，还要控制远程必要节点、消息数量和空间局部性。
1.5D 复制把带宽降到 (O(n/c))、延迟降到 (O(p/c^2))，体现 communication avoiding 的通用思想。