Scalable parallel 3-D TEM inversion via rational approximation of the matrix exponential

Resumen Técnico: Inversión Paralela Escalable 3-D TEM mediante Aproximación Racional del Exponencial de Matriz

Enunciado del Problema
La Electromagnetismo Transitorio (TEM) es una técnica crítica en geofísica aplicada para la exploración mineral, la hidrogeología y los estudios ambientales. Sin embargo, la inversión de modelos de conductividad 3-D en entornos no triviales es computacionalmente exigente. Los enfoques tradicionales dependen de esquemas de paso de tiempo implícitos (por ejemplo, Euler hacia atrás) para resolver las ecuaciones de Maxwell en el dominio del tiempo. Estos métodos organizan el flujo de trabajo computacional en torno a la marcha temporal, donde el costo de la modelación directa y la inversión está estrechamente acoplado al número de pasos de tiempo necesarios para resolver la respuesta transitoria. A medida que las ganancias de rendimiento de un solo núcleo se estancan debido a las limitaciones de la Ley de Moore, estos algoritmos inherentemente secuenciales luchan por escalar de manera eficiente en hardware paralelo moderno.

Metodología
Los autores proponen un nuevo marco de inversión que desacopla el costo computacional del número de tiempos de observación mediante la reformulación del problema directo utilizando la acción de un exponencial de matriz.

1. Aproximación Racional del Exponencial de Matriz:
   En lugar de usar pasos de tiempo, la solución transitoria $u(t)$ se expresa como la acción del exponencial de matriz $\exp(-t M^{-1}K)$ sobre un vector inicial. Los autores emplean una aproximación racional casi óptima (RBA) de este exponencial, construida utilizando el algoritmo RKFIT. Esta aproximación representa la solución como una combinación lineal de sistemas lineales desplazados:
   $$u(t_j) \approx 2 \text{Re} \left( \sum_{i=1}^m \alpha_i^{(j)} (K - \xi_i M)^{-1} f \right)$$
   donde $\{ \xi_i \}$ son polos compartidos en todos los puntos de tiempo, y $\{ \alpha_i^{(j)} \}$ son coeficientes precalculados. Crucialmente, las matrices del sistema desplazado $(K - \xi_i M)$ son independientes del tiempo.

2. Cálculo Paralelo Directo y de Sensibilidades:
   La reformulación permite que los $m$ sistemas lineales desplazados se resuelvan de forma independiente y en paralelo. La dependencia temporal está contenida enteramente en los residuos (coeficientes $\alpha$), lo que significa que el cálculo de las respuestas directas y las sensibilidades es efectivamente independiente del número de tiempos de observación deseados $K$.

   • Solver Directo: Utiliza solvers directos paralelos (MUMPS) para factorizar y resolver los sistemas desplazados.
   • Sensibilidades (Jacobian): El Jacobiano se calcula implícitamente utilizando el mismo conjunto de sistemas desplazados. La derivada de la matriz de masa se aplica a los vectores de solución precalculados $g_i = (K - \xi_i M)^{-1}f$. Esto evita las recursiones secuenciales directas/adjuntas típicas de los métodos de paso de tiempo.

3. Marco de Inversión:

   • Optimización: El problema inverso se resuelve utilizando un marco de Gauss–Newton. El problema de mínimos cuadrados linealizado se resuelve mediante el algoritmo LSQR, que maneja el Jacobiano como un operador implícito (libre de matrices).
   • Regularización: Las restricciones de suavidad se formulan utilizando elementos Raviart–Thomas (RT0) para manejar parámetros de modelo constantes por tramos en mallas no estructuradas, asegurando una discretización consistente del gradiente.
   • Implementación: El código está escrito en Julia, utilizando el marco Gridap.jl para el ensamblaje de elementos finitos y una interfaz MPI personalizada para MUMPS para gestionar factorizaciones paralelas en arquitecturas de memoria compartida.

Contribuciones Clave

• Cambio de Paradigma en el Diseño Algorítmico: El artículo demuestra un cambio de adaptar métodos de marcha temporal secuenciales al hardware paralelo, a diseñar formulaciones matemáticas (RBA) que exponen inherentemente la concurrencia.
• Desacoplamiento de Canales de Tiempo: El costo computacional está determinado por el número de polos $m$ (típicamente pequeño, por ejemplo, 16–21) en lugar del número de canales de tiempo $K$ (a menudo cientos). Esto permite un muestreo temporal de alta resolución sin aumentar la carga de trabajo de álgebra lineal.
• Flujo de Trabajo de Inversión Escalable: Al incrustar el paralelismo a nivel del modelo matemático, los autores logran una inversión escalable donde el costo dominante (resolver sistemas desplazados) puede distribuirse entre rangos MPI.
• Lanzamiento de Código Abierto: La implementación se lanza como software de código abierto para apoyar la investigación reproducible en la comunidad geofísica.

Resultados
El método se validó en un modelo sintético 3-D con aproximadamente 700.000 grados de libertad, que presentaba anomalías conductoras y resistivas en un semiespacio.

• Precisión de Reconstrucción: La inversión recuperó con éxito las ubicaciones y formas generales de las anomalías. Los bloques conductores se resolvieron bien, mientras que los bloques resistivos mostraron un contraste reducido, consistente con las limitaciones conocidas de la sensibilidad TEM y la regularización de suavidad.
• Rendimiento:
  • Escalado: Los experimentos de escalado fuerte demostraron un rendimiento casi ideal hasta 8 rangos MPI tanto para la factorización como para la sustitución hacia atrás. La eficiencia disminuyó ligeramente a 16 rangos debido a efectos NUMA (Acceso a Memoria No Uniforme) de memoria compartida y contención de ancho de banda de memoria.
  • Memoria: El consumo máximo de memoria fue de aproximadamente 620 GB, necesario para almacenar las factorizaciones de 21 sistemas desplazados simultáneamente.
  • Comparación: En comparación con el paso de tiempo Euler implícito (que requiere 31 factorizaciones de valor real para la misma ventana de tiempo), el enfoque RBA requiere 21 factorizaciones de valor complejo. Si bien la aritmética compleja es más costosa, el método RBA elimina la ruta crítica de marcha temporal secuencial, ofreciendo un potencial superior para la reducción del tiempo de reloj en arquitecturas distribuidas.

Significado y Afirmaciones
Los autores afirman que este trabajo proporciona una alternativa escalable a la inversión estándar 3-D TEM, abordando específicamente el cuello de botella de la recursión temporal. La ventaja principal no es necesariamente una reducción en el trabajo aritmético serial, sino la reducción de la ruta crítica paralela y la independencia del costo respecto al número de canales de tiempo.

El artículo reconoce las limitaciones actuales, señalando que el almacenamiento simultáneo de múltiples factorizaciones dispersas crea una huella de memoria significativa y que las arquitecturas de memoria compartida enfrentan cuellos de botella de ancho de banda. Los autores sugieren que el trabajo futuro debería centrarse en la transición a clústeres de computación de alto rendimiento (HPC) distribuidos y, potencialmente, en reemplazar los solvers directos con solvers iterativos preconicionados para mejorar aún más la escalabilidad y reducir el uso de memoria. El estudio concluye que el marco propuesto recupera con éxito estructuras 3-D y establece una base para inversiones electromagnéticas más eficientes y a gran escala.