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

Dichiarazione del Problema

Il documento affronta il problema inverso di ricostruire un coefficiente di potenziale $q$ in un'equazione differenziale alle derivate parziali (PDE) ellittica a partire da osservazioni distribuite e rumorose della soluzione $u$. L'equazione governante è:
$$\begin{cases} -\Delta u + qu = f & \text{in } \Omega, \\ u = 0 & \text{su } \partial\Omega, \end{cases}$$
dove $\Omega \subset \mathbb{R}^s$ ($s=1,2$) è un dominio limitato, $f$ è una sorgente nota e il potenziale $q$ appartiene a un insieme ammissibile di funzioni limitate. Le osservazioni sono date da $u^\eta(x) = u(x) + \eta(x)$, dove $\eta$ rappresenta il rumore di misura.

La sfida centrale risiede nella natura della funzione incognita $q$. In molte applicazioni pratiche (ad es. imaging, elastografia), $q$ può presentare brusche discontinuità o caratteristiche non regolari. Gli approcci bayesiani convenzionali basati su prior gaussiani spesso falliscono in questi scenari perché impongono una eccessiva regolarità, portando a ricostruzioni sfocate che non riescono a catturare salti bruschi o discontinuità a gradino. Sebbene i prior di Variazione Totale (TV) affrontino questo problema promuovendo soluzioni a costanza a tratti, essi presentano dei limiti nel definire distribuzioni su spazi topologici generali e nel catturare informazioni strutturali complesse.

Metodologia

Gli autori propongono un nuovo prior bayesiano basato sulla persistenza dell'omologia (PH) all'interno di un framework ibrido. La metodologia è composta dai seguenti elementi:

1. Framework Bayesiano:
   Il problema è formulato utilizzando il teorema di Bayes. La distribuzione a posteriori $\mu_{post}$ è derivata da una funzione di verosimiglianza (basata su rumore gaussiano) e da una misura prior $\mu_{pr}$. Per garantire che il prior sia ben definito in spazi a dimensione infinita, gli autori costruiscono il prior rispetto a una misura di riferimento gaussiana $\mu_0$ (tipicamente definita su $L^2(\Omega)$). La densità a posteriori è espressa tramite la derivata di Radon-Nikodym:
   $$\frac{d\mu_{post}}{d\mu_0}(q) \propto \exp(-\Phi(q; u^\eta) - R(q))$$
   dove $\Phi$ è il mismatch dei dati e $R(q)$ è un termine di regolarizzazione.

2. Costruzione della Persistenza dell'Omologia (PH):

   • Caratteristiche Topologiche: Gli autori utilizzano l'omologia persistente per quantificare le caratteristiche topologiche delle funzioni candidate. Discretizzando il dominio in un complesso simpliciale e definendo insiemi di sottolivello, tracciano la "nascita" e la "morte" delle classi di omologia (componenti connesse in 1D, cicli in 2D) al variare della soglia del valore della funzione.
   • Coppie di Persistenza: Queste caratteristiche sono codificate come coppie di persistenza $(\sigma, \tau)$, dove $\sigma$ è il simplesso di nascita e $\tau$ è il simplesso di morte. La persistenza è definita come $|q(\tau) - q(\sigma)|$.
   • Distanza di Persistenza: Viene definita una metrica $\|q\|_{per}$ come la somma delle persistenze di tutte le coppie. In 1D, si dimostra che questa è strettamente correlata alla Variazione Totale (TV), specificamente $\|q\|_{per} + \text{range globale} = TV(q)$.

3. Il Prior PH-Gaussiano:
   Gli autori introducono un termine di regolarizzazione $R_{PH}(q)$ basato sulla distanza di persistenza. La densità del prior è definita come:
   $$\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)$$

   • Strategia di Ponderazione: Il termine $\alpha_j(q)$ è una funzione di peso dipendente dall'entità della persistenza. Le piccole persistenze (probabilmente rumore) sono pesantemente penalizzate (grande $\alpha_j$), mentre le grandi persistenze (caratteristiche strutturali significative) sono penalizzate meno (piccolo $\alpha_j$).
   • Generalizzazione: Se $\alpha_j \equiv 1$, il prior si riduce allo standard prior TV. Pertanto, il prior PH è presentato come una generalizzazione della TV che incorpora una ponderazione topologica.
   • Estensione 2D: Per domini quadrati 2D, il prior è costruito applicando il prior PH 1D a ciascuna riga e colonna della matrice della funzione discretizzata e sommandone i risultati.

4. Algoritmo di Campionamento:
   Per campionare dalla distribuzione a posteriori, gli autori utilizzano un algoritmo di Metropolis-Hastings (MH) ad accettazione ritardata. Il funzionale del potenziale $V(q)$ è decomposto in due parti: un termine di mismatch dei dati computazionalmente costoso e un termine di regolarizzazione più economico. L'algoritmo utilizza un meccanismo di accettazione a due stadi:

   • Stadio 1: Un passaggio di pre-screening utilizzando il termine di regolarizzazione $R(q)$ per filtrare le proposte.
   • Stadio 2: Un passaggio di correzione utilizzando il potenziale a posteriori completo se la proposta supera il primo stadio.
     Questo approccio migliora l'efficienza computazionale, specialmente quando l'operatore forward è costoso da valutare.

Contributi Chiave

• Nuova Formulazione del Prior: Il documento introduce un prior ibrido che combina il framework probabilistico delle misure gaussiane con vincoli topologici derivati dall'omologia persistente.
• Regolarizzazione Topologica: A differenza dei classici prior TV, il prior PH codifica esplicitamente la struttura topologica della funzione incognita, distinguendo tra rumore (caratteristiche topologiche a breve durata) e caratteristiche significative (caratteristiche a lunga durata) attraverso una distanza di persistenza ponderata.
• Generalizzazione della TV: Gli autori dimostrano che il prior PH è una generalizzazione del prior TV, capace di operare su spazi topologici più generali e di preservare informazioni strutturali più ricche.
• Campionamento Efficiente: Lo sviluppo di un algoritmo ad accettazione ritardata adattato alla decomposizione del potenziale PH-gaussiano permette un campionamento efficiente in problemi inversi ad alta dimensionalità.

Risultati Numerici

Gli autori validano il metodo attraverso esperimenti numerici sia in contesti 1D che 2D, confrontando il prior PH-gaussiano con:

1. Standard Gaussian priors (con vari operatori di covarianza).
2. Prior TV-gaussiani.

Risultati:

• Funzioni Non Regolari: Nei casi che coinvolgono funzioni a costanza a tratti o funzioni con salti bruschi (Esempi 3, 4, 5, 6, 7), i prior gaussiani standard falliscono nel catturare le discontinuità, producendo risultati oscillatori o eccessivamente smussati. Sia i prior TV-gaussiani che i PH-gaussiani superano significativamente il prior gaussiano standard.
• Funzioni Oscillanti: Nell'Esempio 1 (una funzione di tipo Weierstrass), il prior PH-gaussiano ha dimostrato una precisione superiore rispetto al TV-gaussiano, con tracce di errore visibilmente inferiori. Il prior PH ha efficacemente soppresso il rumore ad alta frequenza preservando al contempo la complessa struttura oscillante.
• Ricostruzione 2D: Negli esempi 2D che coinvolgono discontinuità circolari, il prior PH-gaussiano ha fornito ricostruzioni comparabili o leggermente migliori rispetto al TV-gaussiano, in particolare nella gestione di strutture topologiche complesse (Esempio 7).
• Costo Computazionale: L'algoritmo ad accettazione ritardata è risultato leggermente più veloce degli schemi pCN standard. Sebbene il prior PH sia computazionalmente più costoso da valutare rispetto al prior TV, la migliore accuratezza giustifica tale costo in specifici scenari non regolari.

Significatività e Rivendicazioni

Il documento sostiene che il proposto prior basato sulla PH offre un miglioramento modesto ma consistente rispetto al classico prior TV, in particolare per funzioni con caratteristiche topologiche complesse o oscillazioni ad alta frequenza.

Gli autori sottolineano che il vantaggio principale del loro approccio è la sua applicabilità: definendo il prior rispetto a una misura di riferimento gaussiana e richiedendo solo che le funzioni incognite appartengano a uno spazio topologico appropriato, il metodo evita i severi vincoli di regolarità dei metodi classici. Questa flessibilità rende il prior PH uno strumento robusto per problemi inversi in cui la funzione target è non regolare e possiede informazioni strutturali intricate che i normali prior di regolarità non riescono a catturare. Il lavoro posiziona il prior PH non come un sostituto di tutti i metodi esistenti, ma come un'estensione specializzata e topologicamente informata del prior TV per problemi inversi di potenziale particolarmente impegnativi.