Goal-Oriented Real-Time Bayesian Inference for Linear Autonomous Dynamical Systems With Application to Digital Twins for Tsunami Early Warning

Resumen Técnico: Inferencia Bayesiana en Tiempo Real Orientada a Objetivos para Sistemas Dinámicos Autónomos Lineales

Planteamiento del Problema
El artículo aborda el desafío de resolver problemas inversos lineales de gran escala gobernados por modelos de ecuaciones diferenciales parciales (EDP) de alta fidelidad dentro de un marco bayesiano, específicamente para aplicaciones en tiempo real. La aplicación objetivo es un gemelo digital para la alerta temprana de tsunamis en zonas de subducción, como la Zona de Subducción de Cascadia (CSZ). La tarea específica consiste en inferir el movimiento del lecho marino espaciotemporal inducido por terremotos (el parámetro $m$) a partir de observaciones dispersas de presión en el fondo marino (los datos $d$) y, posteriormente, predecir las alturas de las olas del tsunami (cantidades de interés, o QoIs, denotadas como $q$) en regiones costeras vulnerables.

El problema se caracteriza por una dimensionalidad extrema: un modelo de EDP de ondas acústicas-gravitatorias discretizado para la CSZ genera aproximadamente $O(10^9)$ parámetros de inversión. Los métodos tradicionales de inferencia bayesiana, que dependen del muestreo o de aproximaciones de Laplace de la posterior, son computacionalmente intratables para tales escalas. Además, las técnicas de modelado de sustitutos (surrogates) estándar (por ejemplo, procesos gaussianos, operadores neuronales) son ineficaces porque el mapa de parámetro-a-observable para problemas de propagación de ondas no admite una representación de subespacio de baja dimensión; los valores singulares del mapa decaen lentamente debido a la naturaleza de alta frecuencia de las ecuaciones diferenciales parciales hiperbólicas que las gobiernan. Adicionalmente, la necesidad de ejecución en tiempo real (dentá segundos) imposibilita el uso de resolvedores iterativos convencionales que requieren múltiples resoluciones de la EDP directa y adjunta durante la fase de inferencia.

Metodología
Los autores proponen un marco que explota la estructura autónoma del sistema dinámico directo (específicamente, un sistema lineal e invariante en el tiempo o LTI) para descomponer la solución en una fase de precomputación fuera de línea (offline) y una fase de ejecución en tiempo real (online). La idea matemática central es que los mapas de parámetro-a-observable (p2o) y parámetro-a-QoI (p2q) para tales sistemas exhiben una estructura de bloque Toeplitz debido a la invarianza de desplazamiento temporal.

La metodología se desarrolla a través de cuatro fases distintas:

1. Fase 1 (Offline): Precomputación de Mapas.
   En lugar de resolver el problema directo para cada parámetro, el marco computa el mapa p2o ($F$) y el mapa p2q ($F_q$) resolviendo las ecuaciones de la EDP adjunta. Debido a la estructura de bloque Toeplitz, solo se requiere una resolución de la adjunta por cada ubicación de sensor (para $F$) y por cada ubicación de QoI (para $F_q$). Estas soluciones se desplazan en el tiempo para construir los mapas completos. Este paso es computacionalmente costoso pero se realiza una sola vez.

2. Fase 2 (Offline): Representación Compacta de la Posterior.
   La covarianza posterior de los parámetros, $\Gamma_{post}(m) = (F^* \Gamma_{noise}^{-1} F + \Gamma_{prior}^{-1})^{-1}$, es demasiado grande para almacenarse directamente ($O(10^9 \times 10^9)$). Los autores emplean la fórmula de Sherman–Morrison–Woodbury para desplazar la inversión desde el espacio de parámetros de alta dimensión hacia el espacio de datos de menor dimensión ($O(N_d N_t \times N_d N_t)$). Esto requiere precomputar una matriz $K = \Gamma_{noise} + F \Gamma_{prior} F^*$ y una matriz $G^* = \Gamma_{prior} F^*$. Crucialmente, si la covarianza de la previa es invariante en el tiempo, $G^*$ también hereda la estructura de bloque Toeplitz, permitiendo un almacenamiento y cómputo eficientes.

3. Fase 3 (Offline): Precomputación Orientada a Objetivos.
   Para permitir la predicción en tiempo real de las QoIs sin inferir explícitamente todo el campo de parámetros, el marco precomputa un mapa de datos-a-QoI (d2q) ($Q$). Este mapa relaciona directamente los datos observados con la media y la covarianza posterior de la QoI, evitando el paso intermedio de la inferencia de parámetros durante la fase online.

4. Fase 4 (Online): Inferencia y Predicción en Tiempo Real.
   Dados nuevos datos observacionales, el sistema realiza:

   • Inferencia de Parámetros: Computa la estimación de Máxima Verosimilitud a Posteriori (MAP) del movimiento del lecho marino utilizando las matrices precomputadas y multiplicaciones matriz-vector (matvecs) rápidas.
   • Predicción de QoI: Computa directamente las alturas de las olas del tsunami y sus incertidumbres utilizando el mapa d2q precomputado.

   El cuello de botella computacional en la fase online se elimina mediante un algoritmo de Transformada Rápida de Fourier (FFT) acelerado por múltiples GPUs. Este algoritmo explota la estructura de bloque Toeplitz para realizar matvecs con $F$, $F^*$ y sus adjuntas en tiempo $O(N \log N)$, evitando cualquier resolución de la EDP durante la fase online.

Contribuciones Clave

• Inferencia Exacta en Tiempo Real: El artículo presenta un marco que resuelve el problema inverso discretizado de forma exacta (hasta el error de redondeo) sin aproximaciones ni modelos sustitutos, logrando un rendimiento en tiempo real para problemas con $O(10^8)$ a $O(10^9)$ parámetros.
• Descomposición Offline-Online: Al aprovechar la naturaleza autónoma del sistema, las resoluciones de la EDP computacionalmente intensivas se trasladan enteramente a una fase offline. La fase online depende únicamente de operaciones de álgebra lineal sobre matrices precomputadas.
• Aceleración Basada en FFT: Los autores introducen un método para computar matvecs de la Hessiana utilizando algoritmos basados en FFT en clústeres de múltiples GPUs, logrando aceleraciones de aproximadamente $44,000\times$ a $160,000\times$ en comparación con los métodos convencionales libres de matriz (matrix-free) que involucran resoluciones de la EDP directa/adjunta.
• Eficiencia Orientada a Objetivos: El marco apunta explícitamente a la predicción de las QoIs, permitiendo el mapeo directo de los datos a las QoIs (vía el mapa d2q), lo que reduce drásticamente los requisitos de almacenamiento y cómputo en el paso de predicción online.
• Manejo de Hessianas de Alto Rango: El método aborda con éxito problemas inversos donde la Hessiana tiene un decaimiento lento de sus valores singulares (alto rango), un escenario donde las técnicas estándar de aproximación de bajo rango fallan.

Resultos
Los autores demuestran el marco en una configuración de prueba 3D de las ecuaciones de ondas acústicas-gravitatorias que representa un escenario de zona de subducción:

• Escala: El problema de prueba involucra aproximadamente 132 millones de parámetros ($1.3 \times 10^8$), 24,500 puntos de datos y 800 QoIs.
• Rendimiento:
  • Offline: La precomputación de los mapas y matrices necesarios requirió resolver 65 EDPs adjuntas, tomando aproximadamente 8.6 horas en un clúster de supercomputación.
  • Online: El proceso completo de extremo a extremo de inferir el movimiento del lecho marino y predecir las alturas de las olas del tsunami (incluyendo la cuantificación de la incertidumbre) se completó en fracciones de segundo (específicamente, ~57 ms para la inferencia y ~1 ms para la predicción de la QoI en una computadora portátil para el paso de la QoI).
• Precisión: Para datos sintéticos con un 6% de ruido aditivo, el error relativo en el movimiento del lecho marino inferido fue de aproximadamente 9.5%, mientras que el error relativo en las alturas de las olas del tsunami predichas (QoIs) fue de solo ~1.6%. Los resultados indican que, si bien los componentes de alta frecuencia del movimiento del lecho marino están pobremente resueltos (debido a los datos dispersos), los componentes de baja frecuencia que impulsan la formación de tsunamis se capturan con precisión, lo que conduce a predicciones de QoI altamente exactas.
• Comparación: En comparación con un enfoque convencional de gradiente conjugado, que requeriría días para resolver el mismo problema, el método propuesto ofrece una aceleración de aproximadamente 400 millones de veces para la etapa online.

Significancia y Reivindicaciones
El artículo afirma proporcionar un avance fundamental en la construcción de gemelos digitales para la infraestructura crítica y la gestión de desastres. Al permitir la inferencia bayesiana orientada a objetivos que es tanto exacta (fiel al modelo de la EDP de alta fidelidad) como en tiempo real, el marco hace que sea factible construir gemelos digitales para sistemas previamente considerados intratables debido a su escala y restricciones computacionales.

Los autores enfatizan que su enfoque no depende de sustitutos basados en datos o aproximaciones de bajo rango, que suelen ser insuficientes para los problemas de propagación de ondas. En su lugar, explotan la estructura matemática específica de los sistemas dinámicos autónomos para lograr la eficiencia. El objetivo final es aplicar este marco a la Zona de Subducción de Cascadia para proporcionar alertas tempranas de tsunamis generados por terremotos de megacabalgamiento, potencialmente salvando vidas al reducir el tiempo requerido para la inferencia y la predicción de horas a segundos. El artículo señala que, aunque el caso de prueba actual es más pequeño que la escala completa de la CSZ, la metodología está diseñada para escalar al régimen de $O(10^9)$ parámetros requerido para el gemelo digital completo.