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

Problem Statement

The paper addresses the inverse problem of reconstructing a potential coefficient $q$ in an elliptic partial differential equation (PDE) from distributed, noisy observations of the solution $u$. The governing equation is:
$$\begin{cases} -\Delta u + qu = f & \text{in } \Omega, \\ u = 0 & \text{on } \partial\Omega, \end{cases}$$
where $\Omega \subset \mathbb{R}^s$ ($s=1,2$) is a bounded domain, $f$ is a known source, and the potential $q$ belongs to an admissible set of bounded functions. The observations are given by $u^\eta(x) = u(x) + \eta(x)$, where $\eta$ represents measurement noise.

The core challenge lies in the nature of the unknown function $q$. In many practical applications (e.g., imaging, elastography), $q$ may exhibit sharp discontinuities or non-smooth features. Conventional Bayesian approaches relying on Gaussian priors often fail in these scenarios because they enforce excessive smoothness, leading to blurred reconstructions that cannot capture abrupt jumps or step discontinuities. While Total Variation (TV) priors address this by promoting piecewise constant solutions, they have limitations in defining distributions on general topological spaces and capturing complex structural information.

Methodology

The authors propose a novel Persistent Homology (PH)-based Bayesian prior within a hybrid framework. The methodology consists of the following components:

1. Bayesian Framework:
   The problem is formulated using Bayes' theorem. The posterior distribution $\mu_{post}$ is derived from a likelihood function (based on Gaussian noise) and a prior measure $\mu_{pr}$. To ensure the prior is well-defined in infinite-dimensional spaces, the authors construct the prior with respect to a Gaussian reference measure $\mu_0$ (typically defined on $L^2(\Omega)$). The posterior density is expressed via the Radon-Nikodym derivative:
   $$\frac{d\mu_{post}}{d\mu_0}(q) \propto \exp(-\Phi(q; u^\eta) - R(q))$$
   where $\Phi$ is the data misfit and $R(q)$ is a regularization term.

2. Persistent Homology (PH) Construction:

   • Topological Features: The authors utilize persistent homology to quantify the topological features of candidate functions. By discretizing the domain into a simplicial complex and defining sublevel sets, they track the "birth" and "death" of homology classes (connected components in 1D, loops in 2D) as the function value threshold increases.
   • Persistence Pairs: These features are encoded as persistence pairs $(\sigma, \tau)$, where $\sigma$ is the birth simplex and $\tau$ is the death simplex. The persistence is defined as $|q(\tau) - q(\sigma)|$.
   • Persistence Distance: A metric $\|q\|_{per}$ is defined as the sum of persistences of all pairs. In 1D, this is shown to be closely related to the Total Variation (TV), specifically $\|q\|_{per} + \text{global range} = TV(q)$.

3. The PH-Gaussian Prior:
   The authors introduce a regularization term $R_{PH}(q)$ based on the persistence distance. The prior density is defined as:
   $$\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)$$

   • Weighting Strategy: The term $\alpha_j(q)$ is a weight function dependent on the persistence magnitude. Small persistences (likely noise) are heavily penalized (large $\alpha_j$), while large persistences (significant structural features) are penalized less (small $\alpha_j$).
   • Generalization: If $\alpha_j \equiv 1$, the prior reduces to the standard TV prior. Thus, the PH prior is presented as a generalization of TV that incorporates topological weighting.
   • 2D Extension: For 2D square domains, the prior is constructed by applying the 1D PH prior to each row and column of the discretized function matrix and summing the results.

4. Sampling Algorithm:
   To sample from the posterior, the authors employ a Delayed Acceptance Metropolis-Hastings (MH) algorithm. The potential functional $V(q)$ is decomposed into two parts: a computationally expensive data misfit term and a cheaper regularization term. The algorithm uses a two-stage acceptance mechanism:

   • Stage 1: A pre-screening step using the regularization term $R(q)$ to filter proposals.
   • Stage 2: A correction step using the full posterior potential if the proposal passes the first stage.
     This approach improves computational efficiency, particularly when evaluating the forward operator is expensive.

Key Contributions

• Novel Prior Formulation: The paper introduces a hybrid prior that combines the probabilistic framework of Gaussian measures with topological constraints derived from persistent homology.
• Topological Regularization: Unlike standard TV priors, the PH prior explicitly encodes the topological structure of the unknown function, distinguishing between noise (short-lived topological features) and significant features (long-lived features) through a weighted persistence distance.
• Generalization of TV: The authors demonstrate that the PH prior is a generalization of the TV prior, capable of operating on more general topological spaces and preserving richer structural information.
• Efficient Sampling: The development of a delayed acceptance algorithm tailored to the decomposition of the PH-Gaussian potential allows for efficient sampling in high-dimensional inverse problems.

Numerical Results

The authors validate the method through numerical experiments in both 1D and 2D settings, comparing the PH-Gaussian prior against:

1. Standard Gaussian priors (with various covariance operators).
2. TV-Gaussian priors.

Findings:

• Non-Smooth Functions: In cases involving piecewise constant functions or functions with sharp jumps (Examples 3, 4, 5, 6, 7), standard Gaussian priors fail to capture discontinuities, producing oscillatory or overly smooth results. Both TV-Gaussian and PH-Gaussian priors significantly outperform the standard Gaussian prior.
• Oscillatory Functions: In Example 1 (a Weierstrass-type function), the PH-Gaussian prior demonstrated superior accuracy compared to the TV-Gaussian prior, with visibly lower error traces. The PH prior effectively suppressed high-frequency noise while preserving the complex oscillatory structure.
• 2D Reconstruction: In 2D examples involving circular discontinuities, the PH-Gaussian prior provided reconstructions comparable to or slightly better than the TV-Gaussian prior, particularly in handling complex topological structures (Example 7).
• Computational Cost: The delayed acceptance algorithm was found to be slightly faster than standard pCN schemes. While the PH prior itself is computationally more expensive to evaluate than the TV prior, the improved accuracy justifies the cost in specific non-smooth scenarios.

Significance and Claims

The paper claims that the proposed PH-based prior offers a modest yet consistent improvement over the classical TV prior, particularly for functions with complex topological features or high-frequency oscillations.

The authors emphasize that the primary advantage of their approach is its applicability: by defining the prior with respect to a Gaussian reference measure and requiring only that unknown functions belong to a suitable topological space, the method avoids the strict regularity constraints of classical approaches. This flexibility makes the PH prior a robust tool for inverse problems where the target function is non-smooth and possesses intricate structural information that standard smoothness priors cannot capture. The work positions the PH prior not as a replacement for all existing methods, but as a specialized, topologically informed extension of the TV prior for challenging inverse potential problems.