[論文レビュー] Partial inversion of the elliptic operator to speed up computation of likelihood in Bayesian inference

Alexander Litvinenko|arXiv (Cornell University)|Aug 7, 2017

Seismic Imaging and Inversion Techniques参考文献 57被引用数 4

ひとこと要約

本稿では、楕円型作用素の部分的逆行列を用いて、全システムを逆行列にしなくても、点値、界面解、または部分領域の平均値などの必要な解関数の機能のみを計算できる階層的ドメイン分割（HDD）手法を提案する。有限要素法（FEM）と$χ$-行列圧縮を組み合わせることで、計算コストを$χ(n^3)$から$χ(n \log^2 n)$に削減し、ベイズ推論およびデータ同化における尤度評価を著しく高速化しながら、高い精度を維持する。

ABSTRACT

Often, when solving forward, inverse or data assimilation problems, only a part of the solution is needed. As a model, we consider the stationary diffusion problem. We demonstrate an algorithm that can compute only a part or a functional of the solution, without calculating the full inversion operator and the complete solution. It is a well-known fact about partial differential equations that the solution at each discretisation point depends on the solutions at all other discretisation points. Therefore, it is impossible to compute the solution only at one point, without calculating the solution at all other points. The standard numerical methods like a conjugate gradient or Gauss elimination compute the whole solution and/or the complete inverse operator. We suggest a method which can compute the solution of the given partial differential equation 1) at a point; 2) at few points; 3) on an interface; or a functional of the solution, without computing the solution at all points. The required storage cost and computational resources will be lower as in the standard approach. With this new method, we can speed up, for instance, computation of the innovation in filtering or the likelihood distribution, which measures the data misfit (mismatch). Further, we can speed up the solution of the regression, Bayesian inversion, data assimilation, and Kalman filter update problems. Applying additionally the hierarchical matrix approximation, we reduce the cubic computational cost to almost linear $\mathcal{O}(k^2n \log^2 n)$, where $k\ll n$ and $n$ is the number of degrees of freedom. Up to the hierarchical matrix approximation error, the computed solution is exact. One of the disadvantages of this method is the need to modify the existing deterministic solver.

研究の動機と目的

• 楕円型PDEの全解を解かずに、ベイズ推論における尤度計算を高速化すること。
• 解関数の機能（例：点値、界面データ、部分領域の平均）を全解を計算せずに効率的に計算できること。
• 逆問題、データ同化、カルマンフィルタリングにおいて、計算コストとメモリ使用量を削減し、観測可能または関連する部分の解にのみ焦点を当てる。
• 階層的行列（$χ$-行列）圧縮とドメイン分割を統合し、近似的線形の複雑度を達成すること。

提案手法

• 計算ドメインを再帰的に部分領域に分割する階層的ドメイン分割（HDD）を用い、局所的解作用素を計算する。
• シュール補完技術を用いて、部分領域データと界面・境界値を関連付ける界面および境界写像$\Phi^{f}_{\omega}, \Phi^{g}_{\omega}$と$\Psi^{f}_{\omega}, \Psi^{g}_{\omega}$を定義する。
• $\mathcal{H}$-行列近似を中間行列に適用し、$\mathcal{O}(n \log^2 n)$のストレージと$\mathcal{O}(n \log^3 n)$の複雑度でデータスパース表現を達成する。
• 二段階アルゴリズムを採用：『葉から根へ』で局所的作用素を計算し、『根から葉へ』で解関数の機能を再構築する。
• 右辺と剛性行列のための別々の細かいメッシュと粗いメッシュを用いることで、マルチスケール効率を実現する。
• 全解の逆行列を避けることで、ベイズ推論に統合し、全解の逆行列を計算しないで必要なデータ不適合（尤度）関数の機能のみを計算する。

実験結果

リサーチクエスチョン

• RQ1 ベイズ推論における尤度評価に必要な解の一部またはその関数的量を、全PDEを解かずに計算できるか？
• RQ2 階層的ドメイン分割と$\mathcal{H}$-行列圧縮をどのように組み合わせて計算コストを低減しつつ精度を保持できるか？
• RQ3 標準的な反復解法（例：PCG）と比較して、部分的逆作用素の計算にかかる計算およびストレージのオーバーヘッドはどの程度か？
• RQ4 異なる離散化レベルを右辺と剛性行列に用いるマルチスケール問題に、この手法は効率的に適用可能か？
• RQ5 2次元および3次元において、問題サイズの増大に伴い、界面サイズと行列の複雑度に関して、この手法はどのようにスケーリングするか？

主な発見

• HDD手法により、$\mathcal{O}(n^3)$から$\mathcal{O}(n \log^3 n)$への計算コストの削減と、$\mathcal{O}(n^2)$から$\mathcal{O}(n \log^2 n)$へのストレージ削減が$\mathcal{H}$-行列圧縮によって達成された。
• 数値実験では、誤差が$10^{-5}$未満であり、解の品質においてPCGの結果と一致した。
• $n = 513^2$の自由度の場合、HDDは$\mathcal{H}$-コレスキー前処理を施したPCGよりもメモリ使用量が少なく、PCGはメモリオーバーフローにより失敗した。
• すべての$\Phi^{f}_{\omega}, \Phi^{g}_{\omega}$作用素の合計メモリ要件は、$n$に対してほぼ線形にスケーリングされ、データスパース効率を裏付けた。
• HDD手法は、各ツリー段階で全解作用素（$\Phi^{g}, \Phi^{f}$）を計算し、再計算なしに複数の関数的量に再利用可能である。
• 数値結果から、HDDは直接$\mathcal{H}$-行列逆行列よりも効率的であり、PCGよりも解の詳細度において高い精度を達成したが、PCGよりリソースを多く使用した。