[论文解读] 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 该方法在2D与3D中随问题规模增长的可扩展性如何，特别是关于界面大小与矩阵复杂度？

主要发现

• HDD方法结合$\mathcal{H}$-矩阵压缩，将计算成本从$\mathcal{O}(n^3)$降低至$\mathcal{O}(n \log^3 n)$，存储需求从$\mathcal{O}(n^2)$降低至$\mathcal{O}(n \log^2 n)$。
• 该方法具有高精度，在数值测试中误差低于$10^{-5}$，且解的质量与PCG结果一致。
• 对于$ n = 513^2 $自由度的情况，HDD所需的内存少于使用$\mathcal{H}$-Cholesky预条件的PCG，后者因内存溢出而失败。
• 所有$\Phi^{f}_{\omega}, \Phi^{g}_{\omega}$算子的总内存需求几乎随$n$线性增长，证实了其数据稀疏效率。
• HDD方法在每一棵树层级上计算完整的解算子（$\Phi^{g}, \Phi^{f}$），实现对多个泛函的重复使用而无需重新计算。
• 数值结果表明，HDD在效率上优于直接的$\mathcal{H}$-矩阵求逆，且在解的细节精度上优于PCG，尽管其资源消耗高于PCG。