A persistent-homology-based Bayesian prior for potential coefficient reconstruction in an elliptic PDE

技術的要約：楕円型偏微分方程式におけるポテンシャル係数再構成のためのパーシステント・ホモロジーに基づくベイズ事前分布

問題設定

本論文は、分布されたノイズを含む解 $u$ の観測値から、楕円型偏微分方程式（PDE）におけるポテンシャル係数 $q$ を再構成する逆問題を扱う。支配方程式は以下の通りである：
$$\begin{cases} -\Delta u + qu = f & \text{in } \Omega, \\ u = 0 & \text{on } \partial\Omega, \end{cases}$$
ここで、$\Omega \subset \mathbb{R}^s$ ($s=1,2$) は有界領域であり、$f$ は既知のソース、$q$ は有界な関数の許容集合に属する。観測値は $u^\eta(x) = u(x) + \eta(x)$ で与えられ、$\eta$ は測定ノイズを表す。

核心となる課題は、未知の関数 $q$ の性質にある。多くの実用的なアプリケーション（例：イメージング、エラストグラフィ）において、$q$ は鋭い不連続性や非滑らかな特徴を示すことがある。ガウス事前分布に依存する従来のベイズ的手法は、過度な滑らかさを強制するため、こうしたシナリオでは失敗することが多い。その結果、急激なジャンプやステップ不連続を捉えることができず、ぼやけた再構成をもたらしてしまう。全変動（Total Variation; TV）事前分布は、区分的に定数である解を促進することでこの問題に対処するが、一般的な位相空間上の分布を定義することや、複雑な構造情報を捉えることには限界がある。

手法

著者らは、ハイブリッド・フレームワーク内における新しいパーシステント・ホモロジー（PH）に基づくベイズ事前分布を提案している。この手法は、以下のコンポーネントで構成される：

1. ベイズ・フレームワーク:
   問題はベイズの定理を用いて定式化される。事後分布 $\mu_{post}$ は、尤度関数（ガウスノイズに基づく）と事前分布 $\mu_{pr}$ から導出される。無限次元空間において事前分布が適切に定義されることを保証するため、著者らはガウス基準測度 $\mu_0$（通常は $L^2(\Omega)$ 上で定義される）に関する事前分布を構築する。事後密度は、ラドン＝ニコディム微分を用いて次のように表現される：
   $$\frac{d\mu_{post}}{d\mu_0}(q) \propto \exp(-\Phi(q; u^\eta) - R(q))$$
   ここで、$\Phi$ はデータ不適合度、$R(q)$ は正則化項である。

2. パーシステント・ホモロジー（PH）の構築:

   • 位相的特徴: 著者らは、候補関数の位相的特徴を定量化するためにパーシステント・ホモロジーを利用する。領域を単体複体（simplicial complex）に離散化し、劣レベル集合（sublevel sets）を定義することで、関数値の閾値が増加するにつれて、ホモロジー類（1次元における連結成分、2次元におけるループ）の「誕生（birth）」と「消滅（death）」を追跡する。
   • パーシステンス・ペア: これらの特徴は、誕生単体 $\sigma$ と消滅単体 $\tau$ のペア $(\sigma, \tau)$ としてエンコードされる。パーシステンスは $|q(\tau) - q(\sigma)|$ と定義される。
   • パーシステンス距離: メトリック $\|q\|_{per}$ は、すべてのペアのパーシステンスの総和として定義される。1次元においては、これが全変動（TV）と密接に関連していることが示されている。具体的には、$\|q\|_{per} + \text{global range} = TV(q)$ である。

3. PH-ガウス事前分布:
   著者らは、パーシステンス距離に基づく正則化項 $R_{PH}(q)$ を導入する。事前分布密度は次のように定義される：
   $$\frac{d\mu_{pr}}{d\mu_0}(q) \propto \exp\left(-\lambda \sum_{(x_j, \tilde{x}_j) \in P(q)} \alpha_j(q) |q(x_j) - q(\tilde{x}_j)|\right)$$

   • 重み付け戦略: 項 $\alpha_j(q)$ は、パーシステンスの大きさに依存する重み関数である。小さなパーシステンス（ノイズである可能性が高いもの）は厳しく罰せられ（大きな $\alpha_j$）、大きなパーシステンス（重要な構造的特徴）は罰せられにくくなる（小さな $\alpha_j$）。
   • 一般化: $\alpha_j \equiv 1$ の場合、この事前分布は標準的なTV事前分布に帰着する。したがって、PH事前分布は、位相的な重み付けを組み込んだTV事前分布の一般化として提示されている。
   • 2次元への拡張: 2次元の正方形領域に対しては、離散化された関数行列の各行および各列に対して1次元のPH事前分布を適用し、その結果を合計することで事前分布を構築する。

4. サンプリング・アルゴリズム:
   事後分布からサンプリングするために、著者らは**遅延受理メトロポリス・ヘイスティングス（Delayed Acceptance MH）**アルゴリズムを採用している。ポテンシャル関数 $V(q)$ は、計算コストの高いデータ不適合項と、より低コストな正則化項の2つの部分に分解される。アルゴリズムは、以下の2段階の受理メカニズムを使用する：

   • ステージ1: 提案をフィルタリングするための、正則化項 $R(q)$ を用いた事前スクリーニング。
   • ステージ2: 提案が第1ステージを通過した場合に、完全な事後ポテンシャルを用いて行う補正ステップ。
     このアプローチは、フォワード演算子の評価が高コストである場合に、特に計算効率を向上させる。

主な貢献

• 新しい事前分布の定式化: 本論文は、ガウス測度の確率的フレームワークと、パーシステント・ホモロジーから導出される位相的制約を組み合わせたハイブリッドな事前分布を導入している。
• 位相的正則化: 標準的なTV事前分布とは異なり、PH事前分布は未知の関数の位相的構造を明示的に符号化し、重み付きパーシステンス距離を通じて、ノイズ（短命な位相的特徴）と重要な特徴（長命な特徴）を区別する。
• TVの一般化: 著者らは、PH事前分布がTV事前分布の一般化であり、より一般的な位相空間上で動作し、より豊かな構造情報を保持できることを示している。
• 効率的なサンプリング: 高次元の逆問題において効率的なサンプリングを可能にする、遅延受理アルゴリズムの開発。

数値実験の結果

著者らは、1次元および2次元の設定において、PH-ガウス事前分布を以下の手法と比較することで、手法の妥当性を検証した：

1. 標準的なガウス事前分布（様々な共分散演算子を使用）。
2. TV-ガウス事前分布。

知見:

• 非滑らかな関数: 区分的に定数である関数や鋭いジャンプを伴うケース（例3, 4, 5, 6, 7）において、標準的なガウス事前分布は不連続性を捉えることができず、振動的または過度に滑らかな結果を生じさせた。TV-ガウスおよびPH-ガウス事前分布は、いずれも標準的なガウス事前分布よりも大幅に優れた性能を示した。
• 振動関数: 例1（ワイエルシュトラス型の関数）において、PH-ガウス事前分布はTV-ガウス事前分布と比較して優れた精度を示し、誤差の軌跡も明らかに低かった。PH事前分布は、複雑な振動構造を保持しつつ、高周波ノイズを効果的に抑制した。
• 2次元再構成: 円形の不連続を含む2次元の例において、PH-ガラス事前分布は、特に複雑な位相構造を扱う場合（例7）において、TV-ガウス事前分布と同等か、あるいはわずかに優れた再構成を提供した。
• 計算コスト: 遅延受理アルゴリズムは、標準的なpCNスキームよりもわずかに高速であることが判明した。PH事前分布自体は、TV事前分布よりも評価に計算コストがかかるが、特定の非滑らかなシナリオにおける精度の向上がそのコストを正当化している。

意義と主張

本論文は、提案されたPHベースの事前分布が、複雑な位相的特徴や高周波振動を持つ関数に対して、古典的なTV事前分布に対する控えめながらも一貫した改善を提供すると主張している。

著者らは、本手法の主な利点はその適用可能性にあると強調している。未知の関数が適切な位相空間に属することのみを要求し、ガウス基準測量に関して事前分布を定義することで、この手法は古典的なアプローチの厳格な正則性制約を回避している。この柔軟性により、PH事前分布は、標準的な滑らかさの事前分布では捉えられない、非滑らかで複雑な構造情報を持つ逆問題に対する堅牢なツールとなる。本研究は、PH事前分布を、あらゆる既存手法の代替としてではなく、困難な潜在的逆問題に対する、位相的な知見に基づいた専門的なTV事前分布の拡張として位置づけている。