No 3D Matrices: A Unified Tensor-Product View of Matrix-Free Cartesian PDE Solvers

技術要約：3D行列は存在しない：行列フリーな直交座標系PDEソルバーの統一的テンソル積による視点

問題提起
実用的な計算流体力学（CFD）コードやスペクトラル要素法ソルバーは、長年、ある構造的な「秘密」に依拠してきた。それは、直交座標系の偏微分方程式（PDE）（ラプラシアン、質量行列、コンパクト微分スキームなど）を支配する3次元演算子は、決して $N \times N$ のモノリシックな行列（ここで $N = N_x N_y N_z$）として組み立てられたり、因数分解されたり、格納されたりはしない、ということである。教科書ではこれらの演算子が、高密度または疎な3D行列として提示されることが多いが、そのような手法は計算的に極めて困難である。例えば、1辺あたり200点を持つ格子の場合、高密度なラプラシアンは約0.5ペタバイトのストレージを必要とし、たとえ疎な7点ステンシルであっても、直接法による因数分解はワークステーションのメモリを使い果たすことになる。本論文が取り組む問題は、これらの演算子の理論的な表現と、効率的な行列フリーの実装との間にある乖離である。本論文は、これらの3D演算子が1D演算へと「崩壊」する代数的な正当化を統一し、このアプローチをハードウェア効率的にするための具体的な実装テクニックを体系化することを目的としている。

手法
本論文は、クロネッカー積代数に基づく統一的なフレームワークを導出している。核心となる洞察は、任意の3次元直交座標系演算子 $L$ は、1次元演算子のテンソル積の和として表現できるということである：
$$L = \sum_{r} C_z^{(r)} \otimes C_y^{(r)} \otimes C_x^{(r)}$$
ここで、$C_\xi$ は各座標軸方向に作用する小さな帯行列または対角行列である。本手法は、以下の2つの基本的な恒等式に基づいている：

1. 混合積の法則（Mixed-Product Rule）: $(A \otimes B)(C \otimes D) = (AC) \otimes (BD)$。これにより、方向ごとの演算子を、その積を形成することなく合成できる。
2. 逆行列の法則（Inverse Rule）: $(A \otimes B)^{-1} = A^{-1} \otimes B^{-1}$。これにより、テンソル積形式の質量行列やコンパクト因子の逆転が、独立した1D解法へと簡約されることが説明される。

本論文では、この代数を以下の対象に体系的に適用している：

• 差分法: 中央差分および高次ステンシル。
• コンパクト（Padé）スキーム: $A_x u' = R_x u$ のような、高密度な積 $A_x^{-1}R_x$ の形成を回避する暗黙的関係。
• ガラーキン法およびスペクトラル法: B-スプライン、等幾何解析（IGA）、およびコロケーションを含み、質量行列と剛性行列がテンソル積基底から継承するクロネッカー構造を示す。
• 時間ステップ法: 代替方向陰解法（ADI）（Peaceman–Rachford、Douglas–Gunn）は、結合された3D暗黙ステップを、3つの連続する1D帯行列解法へと分解できることを示す。
• 直接ソルバー: ポアソン方程式およびヘルムホルツ方程式に対する高速対角化法（Lynch–Rice–Thomas）において、3D問題が1D演算子の固有基底におけるスカラー除算へとデカップリングされることを示す。

主な貢献
本論文は、単なる代数的な簡約を超えて、ハイパフォーマンス・コンピューティングに必要な「実用的な実装テクニック」へと踏み込み、主に3つの貢献を行っている：

1. 統一的な代数的導出: 有限差分、コンパクトスキーム、ガラーキン法、B-スプライン、コロケーションのすべてに、1D演算への簡約が普遍的に適用できることを示す、クリーンで単一の導出を提供している。これは、「3D演算子」がメモリ上にインスタンス化される必要のない数学的な抽象概念であることを明確にしている。
2. 3つの実用的な実装テクニック:
   • マルチ右辺（Multi-RHS）再整形: $N_y N_z$ 個の独立した1Dラインをループ処理する代わりに、フィールドを $N_x \times (N_y N_z)$ の行列へと再整形する。これにより、操作を単一のバッチカーネルとして露出し、高密度/要素ローカルな因子に対する最適化されたBLAS-3 GEMM呼び出しや、ステンシルに対するバッチ化された帯行列ソルバーの使用を可能にする。
   • 和因子化（Sum Factorization）: 高次スペクトラル要素法において、全3Dテンソル積を行うのではなく、逐次的な1D縮退を行うことで、積分のコストを $O(n^{2d})$ から $O(d n^{d+1})$ （$n$ は方向あたりの点数）へと削減する方法を詳述する。
   • ペンシル分解（Pencil Decomposition）: 分散メモリ（MPI）環境において、格子を一方の方向に整列した「ペンシル（鉛筆状の領域）」に分解する。コレクティブ転置を用いてペンシルの向きを回転させることで、各方向のスイープにおいて、ローカルランク内でデータがメモリ上で連続することを保証する。
3. 計算量分析: 固定されたステンシル幅または多項式次数に対して、計算コストは $O(N)$ であり、演算子のストレージは $O(N_x + N_y + N_z)$ であることを厳密に確立している。これは、$O(N)$ のストレージと $O(N^2)$ または $O(N^{4/3})$ の因数分解コストを要する組み立てられた3D行列とは対照的である。

結果
本論文は、単位立方体上の3Dポアソン方程式の解法に関する数値比較を通じて、その主張を検証している：

• ストレージ: 未知数が $N \approx 1.7 \times 10^7$ ($256^3$) のとき、組み立てられた疎直接ソルバーはギガバイト単位のフィルイン（および標準的なワークステーションでの $48^3$ を超える実行不能状態）を必要とするが、行列フリーの高速対角化法は、キロバイト単位のストレージ（1Dの固有値）のみを必要とする。
• 実行時間: 行列フリーソルバーは $256^3$ の問題を約2秒で完了させるが、組み立てられた直接ソルバーは $48^3$ に対してさえ約1分を要し、より大きな格子では実行不可能になる。
• 精度: 両手法は同一の離散解を与え（誤差 $\approx 10^{-14}$ で一致）、行列フリーのアプローチが精度の低下を招かないことを確認している。

意義
本論文の意義は、半世紀もの間、専門的なコードにおける「明文化されていない習慣」であった構造的な秘密を「公にする」ことにある。クロネッカー積を組織化の原理として枠付けることで、性質の異なる様々な離散化手法（有限差分、スプライン、スペクトラル要素法）を単一のアルゴリズム・パラダイムの下に統合している。本論文は、「3D行列」は表記上の産物に過ぎず、計算に不可欠な要素ではないと主張している。本論文は、これらの簡約を実装するためのマニュアルとして機能し、将来の行列フリー、低ランク、あるいは学習ベースのアクセラレータを評価するための鋭い $O(N)$ の基準を提供している。現代の直交座標系ソルバーの効率性は、複雑な新しいアルゴリズムからではなく、特定のデータレイアウトとカーネル戦略を通じて、グリッドと演算子の分離可能な構造を正しく活用することから生まれるということを強調している。