高速多重極法（FMM）とN体問題

なぜN体問題は総当たりだと破綻するのか

N体問題とは、N個の粒子が互いに及ぼし合う長距離相互作用（重力や静電クーロン力など）を計算し、系の時間発展を追う問題です。重力多体シミュレーションによる銀河形成、分子動力学における静電力、渦法による流体計算などが代表例です。

長距離力の厄介さは、相互作用が距離とともにゆっくり（1/rや1/r²で）しか減衰しないため、原理的にはすべての粒子ペアを無視できない点にあります。ある1粒子が受ける力を正確に求めるには、残りN−1粒子すべてからの寄与を足し合わせる必要があります。これを全粒子について行えば、計算すべきペア数はおよそN×(N−1)/2、すなわちN²のオーダーになります。

直接計算（総当たり）のコスト

  1粒子が受ける力 = Σ_j (j≠i) 粒子jからの寄与
  全粒子ぶん繰り返す → 演算回数 ∝ N²

  例）N = 1,000        → 約 10^6 回
      N = 1,000,000    → 約 10^12 回
      N = 1,000,000,000 → 約 10^18 回

粒子数が1桁増えるたびに演算量は2桁増えます。近距離力（カットオフ距離を超えると無視できる相互作用）ならセルリストや近傍探索で線形時間に落とせますが、これは近距離力を扱う格子ボルツマン法と粒子法と同じ発想であり、遠くまで効く長距離力には使えません。この「遠方も無視できない」という性質こそが、FMMのような階層的近似アルゴリズムを必要とする根本理由です。

核心となる近似 ── 遠い群れはまとめて1個とみなす

FMMとBarnes-Hutの出発点は共通で、素朴だが強力な観測に基づきます。「遠くに固まって存在する粒子の集団は、遠くから見れば1つのかたまりとして扱ってよい」という近似です。

観測点から十分遠い粒子群があるとき、その群れが観測点に及ぼす場（ポテンシャルや力）は、群れ全体の質量・電荷とその分布のモーメントだけで高精度に表現できます。この表現が多重極展開（multipole expansion）です。

FMMではさらにもう一種類の展開を導入します。多重極展開が「遠方の群れが作る場を、群れの側から表す」のに対し、局所展開（local expansion）は「その場を、観測する側の近傍で成り立つ級数として表す」ものです。ある領域に届く遠方からの寄与をいったん1つの局所展開にまとめてしまえば、その領域内のどの点についても、級数を評価するだけで力が求まります。

2種類の展開の役割

  多重極展開(M):  発生源となる粒子群の側に置く
                  「この群れが遠方に作る場」を係数で表現

  局所展開(L):    影響を受ける観測領域の側に置く
                  「遠方から届く場」をその領域内で評価できる級数に変換

  → M から L への変換(M2L)が、FMM を O(N) にする心臓部

木構造による階層化

「どの群れが観測点から十分遠いか」は粒子の配置によって変わります。これを系統的に扱うため、空間を再帰的に分割した木構造を使います。2次元なら四分木（quadtree）、3次元なら八分木（octree）が標準で、1つのセル（ノード）に含まれる粒子数がしきい値を下回るまで、各セルを4分割・8分割し続けます。

八分木による空間分割（3次元）

  ルート = 全空間を含む1つの立方体セル
    ├─ 8個の子セルに分割
    │    ├─ さらに8分割 … （粒子が疎な領域は浅く）
    │    └─ …            （粒子が密な領域は深く）
    └─ 葉セル = 含む粒子数がしきい値未満で分割停止

この木があると、あるセルから見て「遠い／近い」を階層のレベルごとに判定できます。粒子が一様でなく偏って分布する系では木の深さがばらつき、これが並列化での科学シミュレーションの負荷分散を難しくする要因になります。空間充填曲線で木のノードをランクへ割り付ける手法が広く使われます。

Barnes-Hut法 ── 各点から木をたどる O(N log N)

Barnes-Hut法は、木の各セルにそのセル内粒子の重心と総質量（多重極の最低次に相当）を持たせておき、1粒子ずつ力を評価します。評価対象の粒子について、ルートセルから木を下向きにたどり、各セルで開き角度の判定を行います。判定を満たせばそのセルで打ち切り、満たさなければ子セルへ降りて再帰します。

Barnes-Hut の判定（セルサイズ s、粒子までの距離 d）

  s / d < θ  なら … そのセルは「十分遠い」
                     → セルの重心1点として寄与を加え、子は見ない
  s / d ≥ θ  なら … 「近すぎる」
                     → そのセルを開き、子セルを個別に再帰的に調べる

  θ（開き角パラメータ）が小さいほど高精度・低速

遠いセルは中身を展開せず1点で済ませるため、1粒子あたりに調べるセル数は木の深さ、すなわちおよそlog Nのオーダーに収まります。これを全粒子について行うので全体はO(N log N)です。総当たりのN²から劇的な改善であり、実装も比較的単純なため天体物理などで広く使われます。ただしBarnes-Hutはあくまで「粒子ごとに木をたどる」構造であり、対数因子が残ります。

FMMが O(N) を達成する仕組み

FMMがBarnes-Hutを超えて完全な線形時間に届くのは、粒子どうしではなくセルどうしで相互作用を処理し、展開係数を木の上で使い回すからです。FMMは木を2回走査します。

横にスクロール

FMM の2パス構造

  上りパス（葉 → 根、P2M と M2M）
    P2M: 各葉セルで、内部粒子から多重極展開 M を生成
    M2M: 子セルの M を親セルの M へ統合（重心をずらして足し込む）
         → 上位セルほど広い群れの多重極を1個持つ

  下りパス（根 → 葉、M2L と L2L）
    M2L: 「十分離れたセル」の多重極 M を、
         観測側セルの局所展開 L へ変換して足し込む（最重要）
    L2L: 親セルの L を子セルの L へ受け渡す（中心をずらして伝播）

  仕上げ（葉セル）
    L2P: 葉の L を内部各粒子で評価し、遠方からの寄与を得る
    P2P: 隣接する近傍セルの粒子とだけ直接（総当たり）計算

鍵はM2L変換です。あるセルに届く遠方からの寄与は、多数の遠方セルの多重極を1つの局所展開にまとめて蓄えられます。その局所展開は上位から下位へL2Lで伝播し、最後に葉で各粒子に配られます。こうして「遠方からの寄与」は粒子ペア単位ではなくセル単位で一度だけ計算され、木の各レベルで扱うセル数の総和が粒子数に比例するため、全体がO(N)に収まります。近傍のセルどうしだけは近似が効かないので、そこだけP2Pで直接計算します。

どちらを使うか、そして他の選択肢

精度をあまり要求せず実装を単純に保ちたい場面、粒子数が中規模なら、Barnes-HutのO(N log N)で十分なことが多いです。一方、超大規模かつ高精度が要る分子動力学や重力多体では、FMMのO(N)と任意精度制御が効いてきます。

なお、長距離力を階層近似で扱うのはFMM系だけではありません。周期境界を持つ系ではEwald法系のPPPM（particle-particle particle-mesh）が広く使われ、遠方成分を格子上のポアソン方程式に落として並列FFTアルゴリズムで解きます。木ベースのFMMは通信が局所的で任意配置に強い一方、FFTベースの手法は一様格子と周期境界に強く、系の性質と並列規模に応じて使い分けられます。

まとめ

• N体長距離相互作用の直接計算はペア数がN²のオーダーで増え、大規模系では実時間で終わらない。近傍探索での線形化は近距離力にしか使えない。
• 「遠くの粒子群はまとめて多重極展開で近似できる」という観測が出発点で、比 d/r が小さいほど低次の展開で高精度になる。空間を四分木・八分木で階層化し、遠近の判定を系統的に行う。
• Barnes-Hut は各粒子から木を降りて点とセルの相互作用を評価しO(N log N)。FMM は多重極→局所展開の変換（M2L）を導入してセルとセルの相互作用を処理し、展開係数を木上で共有してO(N)を達成する。
• FMM は任意精度に届く一方で定数項が大きく、規模・精度・並列アーキテクチャに応じてBarnes-Hut、FMM、FFTベースのPPPMを使い分けるのが実務上の判断になる。