Scalable parallel 3-D TEM inversion via rational approximation of the matrix exponential

技術的概要：行列指数の有理近似によるスケーラブルな並列 3 次元 TEM 反転

問題提起
過渡電磁気法（TEM）は、鉱物探査、水文地質学、環境研究における応用地質学において重要な手法である。しかし、非自明な環境における 3 次元導電率モデルの反転は、計算集約的である。従来のアプローチは、時間領域でマクスウェル方程式を解くために、陰的な時間ステップ法（例：後退オイラー法）に依存している。これらの手法は時間進行を中心とした計算ワークフローを組織化しており、順問題および反転のコストは、過渡応答を解像するために必要な時間ステップの数と密接に結合している。ムーアの法則の限界によりシングルコアの性能向上が頭打ちとなる中、これらの本質的に逐次的なアルゴリズムは、現代の並列ハードウェア上で効率的にスケーリングすることに苦労している。

手法
著者らは、行列指数の作用を用いて順問題を再定式化することにより、計算コストを観測時間の数から切り離す新しい反転フレームワークを提案する。

1. 行列指数の有理近似:
   時間ステップ法の代わりに、過渡解 $u(t)$ を、初期ベクトルに対する行列指数 $\exp(-t M^{-1}K)$ の作用として表現する。著者らは、RKFIT アルゴリズムを用いて構築された、この指数関数の有理近似（RBA）を採用する。この近似は、解をシフトされた線形方程式系の線形結合として表す：
   $$u(t_j) \approx 2 \text{Re} \left( \sum_{i=1}^m \alpha_i^{(j)} (K - \xi_i M)^{-1} f \right)$$
   ここで、$\{ \xi_i \}$ はすべての時間点にわたって共有される極であり、$\{ \alpha_i^{(j)} \}$ は事前に計算された係数である。重要なのは、シフトされた行列 $(K - \xi_i M)$ が時間に依存しないことである。

2. 並列順問題および感度計算:
   この再定式化により、$m$ 個のシフトされた線形方程式系を独立して、かつ並列に解くことが可能になる。時間依存性は完全に残差（係数 $\alpha$）に含まれるため、順応答および感度の計算は、必要な観測時間数 $K$ に事実上依存しない。

   • 順ソルバ: MUMPS を用いた並列直接ソルバを使用して、シフトされた系を因数分解および解く。
   • 感度（ヤコビアン）: ヤコビアンは、同じシフトされた系を用いて暗黙的に計算される。質量行列の微分は、事前に計算された解ベクトル $g_i = (K - \xi_i M)^{-1}f$ に適用される。これにより、時間ステップ法に典型的な逐次的な順/随伴再帰を回避する。

3. 反転フレームワーク:

   • 最適化: 逆問題はガウス・ニュートン法フレームワークを用いて解かれる。線形化された最小二乗問題は、ヤコビアンを暗黙的演算子（行列フリー）として扱う LSQR アルゴリズムによって解かれる。
   • 正則化: 滑らかさの制約は、非構造化メッシュ上の区画定数モデルパラメータを処理するためにラビアルト・トーマス（RT0）要素を用いて定式化され、勾配の一貫した離散化を確保する。
   • 実装: コードは Julia で記述され、有限要素アセンブリには Gridap.jl フレームワークを、共有メモリアーキテクチャ上の並列因数分解の管理には MUMPS 用のカスタム MPI インターフェースを利用している。

主要な貢献

• アルゴリズム設計のパラダイムシフト: 本論文は、逐次的な時間進行法を並列ハードウェアに適応させることから、本質的に並行性を露出させる数学的定式化（RBA）を設計することへの転換を実証している。
• 時間チャネルの切り離し: 計算コストは、時間チャネル数 $K$（しばしば数百）ではなく、極の数 $m$（典型的には小さく、例：16〜21）によって決定される。これにより、線形代数の作業量を増やすことなく、高解像度の時間サンプリングが可能になる。
• スケーラブルな反転ワークフロー: 数学モデルのレベルで並列性を組み込むことにより、著者らはスケーラブルな反転を実現した。支配的なコスト（シフトされた系の解）は MPI ランク間で分散可能である。
• オープンソース化: 実装は、地質物理学コミュニティにおける再現可能な研究を支援するために、オープンソースソフトウェアとしてリリースされた。

結果
この手法は、半空間内に導電性および抵抗性の異常を有する、約 70 万の自由度を持つ合成 3 次元モデルで検証された。

• 再構成精度: 反転は、異常の位置と全体的な形状を正常に回復した。導電性ブロックはよく解像されたが、抵抗性ブロックは、TEM 感度および滑らかさ正則化の既知の限界と一致して、コントラストが低下した。
• 性能:
  • スケーリング: 強スケーリング実験は、因数分解および後退代入の両方において、8 MPI ランクまでほぼ理想的な性能を示した。共有メモリの NUMA（非一様メモリアクセス）効果およびメモリ帯域幅の競合により、16 ランクでは効率がわずかに低下した。
  • メモリ: ピークメモリ消費量は約 620 GB で、21 個のシフトされた系の因数分解を同時に保存するために必要であった。
  • 比較: 陰的オイラー時間ステップ法（同じ時間窓に対して 31 個の実数値因数分解を必要とする）と比較して、RBA 手法は 21 個の複素数値因数分解を必要とする。複素数演算はより高価であるが、RBA 手法は逐次的な時間進行のクリティカルパスを排除し、分散アーキテクチャにおけるウォールクロック時間の削減において優れた潜在能力を提供する。

意義と主張
著者らは、この研究が時間再帰のボトルネックに特に対処する、標準的な 3 次元 TEM 反転に対するスケーラブルな代替手段を提供すると主張している。主な利点は、必ずしも直列演算作業の削減ではなく、並列クリティカルパスの削減およびコストの時間チャネル数からの独立性にある。

本論文は、複数の疎な因数分解の同時保存が大きなメモリフットプリントを生み出し、共有メモリアーキテクチャが帯域幅のボトルネックに直面しているという現在の限界を認めている。著者らは、将来の研究は、スケーラビリティのさらなる向上とメモリ使用量の削減のために、分散型高性能計算（HPC）クラスターへの移行、および直接ソルバから前処理付き反復ソルバへの置き換えに焦点を当てるべきであると示唆している。本研究は、提案されたフレームワークが 3 次元構造を正常に回復し、より効率的で大規模な電磁気反転の基盤を確立したと結論付けている。