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

Declaración del Problema

El artículo aborda el problema inverso de la reconstrucción de un coeficiente de potencial $q$ en una ecuación diferencial parcial (EDP) elíptica a partir de observaciones distribuidas y con ruido de la solución $u$. La ecuación que rige el proceso es:
$$\begin{cases} -\Delta u + qu = f & \text{en } \Omega, \\ u = 0 & \text{en } \partial\Omega, \end{cases}$$
donde $\Omega \subset \mathbb{R}^s$ ($s=1,2$) es un dominio acotado, $f$ es una fuente conocida y el potencial $q$ pertenece a un conjunto admisible de funciones acotadas. Las observaciones se dan por $u^\eta(x) = u(x) + \eta(x)$, donde $\eta$ representa el ruido de la medición.

El desafío central radica en la naturaleza de la función desconocida $q$. En muchas aplicaciones prácticas (por ejemplo, imagenología, elastografía), $q$ puede exhibir discontinuidades bruscas o características no suaves. Los enfoques bayesianos convencionales que dependen de priors gaussianos suelen fallar en estos escenarios porque imponen una suavidad excesiva, lo que conduce a reconstrucciones difusas que no pueden capturar saltos abruptos o discontinuidades de tipo escalón. Aunque los priors de Variación Total (TV) abordan esto promoviendo soluciones de constante a trozos, tienen limitaciones para definir distribuciones en espacios topológicos generales y para capturar información estructural compleja.

Metodología

Los autores proponen un novedoso prior bayesiano basado en la Homología Persistente (PH) dentro de un marco híbrido. La metodología consta de los siguientes componentes:

1. Marco Bayesiano:
   El problema se formula utilizando el teorema de Bayes. La distribución posterior $\mu_{post}$ se deriva de una función de verosimilitud (basada en ruido gaussiano) y una medida prior $\mu_{pr}$. Para asegurar que el prior esté bien definido en espacios de dimensión infinita, los autores construyen el prior con respecto a una medida de referencia gaussiana $\mu_0$ (típicamente definida en $L^2(\Omega)$). La densidad posterior se expresa mediante la derivada de Radon-Nikodym:
   $$\frac{d\mu_{post}}{d\mu_0}(q) \propto \exp(-\Phi(q; u^\eta) - R(q))$$
   donde $\Phi$ es el desajuste de los datos y $R(q)$ es un término de regularización.

2. Construcción de la Homología Persistente (PH):

   • Características Topológicas: Los autores utilizan la homología persistente para cuantificar las características topológicas de las funciones candidatas. Al discretizar el dominio en un complejo simplicial y definir conjuntos de subnivel, rastrean el "nacimiento" y la "muerte" de las clases de homología (componentes conexas en 1D, lazos en 2D) a medida que aumenta el umbral del valor de la función.
   • Pares de Persistencia: Estas características se codifican como pares de persistencia $(\sigma, \tau)$, donde $\sigma$ es el simplex de nacimiento y $\tau$ es el simplex de muerte. La persistencia se define como $|q(\tau) - q(\sigma)|$.
   • Distancia de Persistencia: Se define una métrica $\|q\|_{per}$ como la suma de las persistencias de todos los pares. En 1D, se demuestra que está estrechamente relacionada con la Variación Total (TV), específicamente $\|q\|_{per} + \text{rango global} = TV(q)$.

3. El Prior PH-Gaussiano:
   Los autores introducen un término de regularización $R_{PH}(q)$ basado en la distancia de persistencia. La densidad del prior se define como:
   $$\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)$$

   • Estrategia de Ponderación: El término $\alpha_j(q)$ es una función de peso dependiente de la magnitud de la persistencia. Las persistencias pequeñas (probablemente ruido) se penalizan fuertmente ( $\alpha_j$ grande), mientras que las persistencias grandes (características estructurales significativas) se penalizan menos ($\alpha_j$ pequeño).
   • Generalización: Si $\alpha_j \equiv 1$, el prior se reduce al prior TV estándar. Por lo tanto, el prior PH se presenta como una generalización de la TV que incorpora una ponderación topológica.
   • Extensión a 2D: Para dominios cuadrados en 2D, el prior se construye aplicando el prior PH de 1D a cada fila y columna de la matriz de la función discretizada y sumando los resultados.

4. Algoritmo de Muestreo:
   Para muestrear de la posterior, los autores emplean un algoritmo de Metropolis-Hastings (MH) de Aceptación Diferida. El potencial de la función $V(q)$ se descompone en dos partes: un término de desajuste de datos computacionalmente costoso y un término de regularización más económico. El algoritmo utiliza un mecanismo de aceptación de dos etapas:

   • Etapa 1: Un paso de pre-criba utilizando el término de regularización $R(q)$ para filtrar las propuestas.
   • Etapa 2: Un paso de corrección utilizando el potencial posterior completo si la propuesta pasa la primera etapa.
     Este enfoque mejora la eficiencia computacional, particularmente cuando la evaluación del operador de propagación (forward operator) es costosa.

Contribuciones Clave

• Nueva Formulación de Prior: El artículo introduce un prior híbrido que combina el marco probabilístico de las medidas gaussianas con restricciones topológicas derivadas de la homología persistente.
• Regularización Topológica: A diferencia de los priors de TV estándar, el prior PH codifica explícicamente la estructura topológica de la función desconocida, distinguiendo entre el ruido (características topológicas de corta duración) y las características significativas (características de larga duración) a través de una distancia de persistencia ponderada.
• Generalización de la TV: Los autores demuestran que el prior PH es una generalización del prior de la TV, capaz de operar en espacios topológicos más generales y de preservar información estructural más rica.
• Muestreo Eficiente: El desarrollo de un algoritmo de aceptación diferida adaptado a la descomposición del potencial PH-gaussiano permite un muestreo eficiente en problemas inversos de alta dimensión.

Resultados Numéricos

Los autores validan el método mediante experimentos numéricos tanto en entornos 1D como 2D, comparando el prior PH-gaussiano contra:

1. Priors gaussianos estándar (con varios operadores de covarianza).
2. Priors de TV-gaussiano.

Hallazgos:

• Funciones No Suaves: En casos que involucran funciones constantes a trozos o funciones con saltos bruscos (Ejemplos 3, 4, 5, 6, 7), los priors gaussianos estándar fallan en capturar las discontinuidades, produciendo resultados oscilatorios o excesivamente suaves. Los priors de TV-gaussiano y PH-gaussiano superan significativamente al prior gaussiano estándar.
• Funciones Oscilatorias: En el Ejemplo 1 (una función de tipo Weierstrass), el prior PH-gaussiano demostró una precisión superior en comparación con el prior TV-gaussiano, con trazas de error visiblemente menores. El prior PH suprimió eficazmente el ruido de alta frecuencia mientras preservaba la compleja estructura oscilatoria.
• Reconstrucción en 2D: En ejemplos 2D que involucran discontinuidades circulares, el prior PH-gaussiano proporcionó reconstrucciones comparables o ligeramente mejores que el prior TV-gaussiano, particularmente al manejar estructuras topológicas complejas (Ejemplo 7).
• Costo Computacional: El algoritmo de aceptación diferida resultó ser ligeramente más rápido que los esquemas pCN estándar. Aunque el propio prior PH es más costoso de evaluar que el prior TV, la mejora en la precisión justifica el costo en escenarios no suaves específicos.

Significancia y Reivindicaciones

El artículo afirma que el prior propuesto basado en PH ofrece una mejora modesta pero consistente sobre el prior clásico de la TV, particularmente para funciones con características topológicas complejas u oscilaciones de alta frecuencia.

Los autores enfatizan que la principal ventaja de su enfoque es su aplicabilidad: al definir el prior con respecto a una medida de referencia gaussiana y requerir solo que las funciones desconocidas pertenezcan a un espacio topológico adecuado, el método evita las estrictas restricciones de regularidad de los enfoques clásicos. Esta flexibilidad hace que el prior PH sea una herramienta robusta para problemas inversos donde la función objetivo no es suave y posee información estructural intrincada que los priors de suavidad estándar no pueden capturar. El trabajo posiciona al prior PH no como un reemplazo para todos los métodos existentes, sino como una extensión especializada y topológicamente informada del prior de la TV para problemas de potencial inverso desafiantes.