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

Resumo Técnico: Inversão Paralela Escalável 3-D de TEM via Aproximação Racional do Exponencial de Matriz

Declaração do Problema
A Eletromagnetismo Transiente (TEM) é uma técnica crítica em geofísica aplicada para exploração mineral, hidrogeologia e estudos ambientais. No entanto, inverter modelos de condutividade 3-D em ambientes não triviais é computacionalmente exigente. Abordagens tradicionais dependem de esquemas de passo de tempo implícitos (por exemplo, Euler implícito) para resolver as equações de Maxwell no domínio do tempo. Esses métodos organizam o fluxo de trabalho computacional em torno da marcha temporal, onde o custo da modelagem direta e da inversão está intimamente acoplado ao número de passos de tempo necessários para resolver a resposta transiente. À medida que os ganhos de desempenho de núcleo único atingem o platô devido às limitações da Lei de Moore, esses algoritmos inerentemente sequenciais lutam para escalar eficientemente em hardware paralelo moderno.

Metodologia
Os autores propõem um novo framework de inversão que desacopla o custo computacional do número de tempos de observação, reformulando o problema direto usando a ação de um exponencial de matriz.

1. Aproximação Racional do Exponencial de Matriz:
   Em vez de passo de tempo, a solução transiente $u(t)$ é expressa como a ação do exponencial de matriz $\exp(-t M^{-1}K)$ sobre um vetor inicial. Os autores empregam uma aproximação racional quase-ótima (RBA) deste exponencial, construída usando o algoritmo RKFIT. Esta aproximação representa a solução como uma combinação linear de sistemas lineares deslocados:
   $$u(t_j) \approx 2 \text{Re} \left( \sum_{i=1}^m \alpha_i^{(j)} (K - \xi_i M)^{-1} f \right)$$
   onde $\{ \xi_i \}$ são polos compartilhados em todos os pontos de tempo, e $\{ \alpha_i^{(j)} \}$ são coeficientes pré-calculados. Crucialmente, as matrizes do sistema deslocado $(K - \xi_i M)$ são independentes do tempo.

2. Cálculo Paralelo Direto e de Sensibilidade:
   A reformulação permite que os $m$ sistemas lineares deslocados sejam resolvidos independentemente e em paralelo. A dependência temporal está inteiramente contida nos resíduos (coeficientes $\alpha$), significando que o cálculo de respostas diretas e sensibilidades é efetivamente independente do número de tempos de observação desejados $K$.

   • Solver Direto: Utiliza solvers diretos paralelos (MUMPS) para fatorizar e resolver os sistemas deslocados.
   • Sensibilidades (Jacobian): O Jacobiano é calculado implicitamente usando o mesmo conjunto de sistemas deslocados. A derivada da matriz de massa é aplicada aos vetores de solução pré-calculados $g_i = (K - \xi_i M)^{-1}f$. Isso evita as recursões diretas/adjuntas sequenciais típicas de métodos de passo de tempo.

3. Framework de Inversão:

   • Otimização: O problema inverso é resolvido usando um framework Gauss–Newton. O problema de mínimos quadrados linearizado é resolvido via o algoritmo LSQR, que trata o Jacobiano como um operador implícito (livre de matriz).
   • Regularização: Restrições de suavidade são formuladas usando elementos Raviart–Thomas (RT0) para lidar com parâmetros de modelo constantes por partes em malhas não estruturadas, garantindo uma discretização consistente do gradiente.
   • Implementação: O código foi escrito em Julia, utilizando o framework Gridap.jl para montagem de elementos finitos e uma interface MPI personalizada para o MUMPS para gerenciar fatorações paralelas em arquiteturas de memória compartilhada.

Principais Contribuições

• Mudança de Paradigma no Design Algorítmico: O artigo demonstra uma mudança de adaptar métodos de marcha temporal sequenciais ao hardware paralelo, para formular formulações matemáticas (RBA) que expõem inerentemente a concorrência.
• Desacoplamento de Canais de Tempo: O custo computacional é determinado pelo número de polos $m$ (tipicamente pequeno, por exemplo, 16–21) e não pelo número de canais de tempo $K$ (frequentemente centenas). Isso permite amostragem temporal de alta resolução sem aumentar a carga de trabalho de álgebra linear.
• Fluxo de Trabalho de Inversão Escalável: Ao incorporar paralelismo no nível do modelo matemático, os autores alcançam uma inversão escalável onde o custo dominante (resolver sistemas deslocados) pode ser distribuído entre ranks MPI.
• Lançamento de Código Aberto: A implementação é lançada como software de código aberto para apoiar a pesquisa reprodutível na comunidade geofísica.

Resultados
O método foi validado em um modelo sintético 3-D com aproximadamente 700.000 graus de liberdade, apresentando anomalias condutivas e resistivas em um semi-espaço.

• Precisão de Reconstrução: A inversão recuperou com sucesso as localizações e formas gerais das anomalias. Blocos condutivos foram bem resolvidos, enquanto blocos resistivos mostraram contraste reduzido, consistente com as limitações conhecidas da sensibilidade do TEM e da regularização de suavidade.
• Desempenho:
  • Escala: Experimentos de escala forte demonstraram desempenho quase ideal até 8 ranks MPI tanto para fatoração quanto para substituição para trás. A eficiência degradou ligeiramente em 16 ranks devido a efeitos NUMA (Acesso à Memória Não Uniforme) de memória compartilhada e contenção de largura de banda de memória.
  • Memória: O consumo máximo de memória foi de aproximadamente 620 GB, necessário para armazenar as fatorações de 21 sistemas deslocados simultaneamente.
  • Comparação: Comparado ao passo de tempo Euler implícito (que requer 31 fatorações de valor real para a mesma janela de tempo), a abordagem RBA requer 21 fatorações de valor complexo. Embora a aritmética complexa seja mais cara, o método RBA elimina o caminho crítico de marcha temporal sequencial, oferecendo potencial superior para redução do tempo de relógio em arquiteturas distribuídas.

Significado e Alegações
Os autores afirmam que este trabalho fornece uma alternativa escalável à inversão padrão 3-D de TEM, abordando especificamente o gargalo da recursão temporal. A principal vantagem não é necessariamente uma redução no trabalho aritmético serial, mas a redução do caminho crítico paralelo e a independência do custo em relação ao número de canais de tempo.

O artigo reconhece limitações atuais, observando que o armazenamento simultâneo de múltiplas fatorações esparsas cria uma pegada de memória significativa e que arquiteturas de memória compartilhada enfrentam gargalos de largura de banda. Os autores sugerem que trabalhos futuros devem focar na transição para clusters de computação de alto desempenho (HPC) distribuídos e potencialmente substituir solvers diretos por solvers iterativos pré-condicionados para melhorar ainda mais a escalabilidade e reduzir o uso de memória. O estudo conclui que o framework proposto recupera com sucesso estruturas 3-D e estabelece uma base para inversões eletromagnéticas de grande escala mais eficientes.