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

Technical Summary: Scalable Parallel 3-D TEM Inversion via Rational Approximation of the Matrix Exponential

Problem Statement
Transient Electromagnetics (TEM) is a critical technique in applied geophysics for mineral exploration, hydrogeology, and environmental studies. However, inverting 3-D conductivity models in non-trivial environments is computationally demanding. Traditional approaches rely on implicit time-stepping schemes (e.g., backward Euler) to solve Maxwell's equations in the time domain. These methods organize the computational workflow around temporal marching, where the cost of forward modeling and inversion is tightly coupled to the number of time steps required to resolve the transient response. As single-core performance gains plateau due to Moore's Law limitations, these inherently sequential algorithms struggle to scale efficiently on modern parallel hardware.

Methodology
The authors propose a novel inversion framework that decouples computational cost from the number of observation times by reformulating the forward problem using the action of a matrix exponential.

1. Rational Approximation of the Matrix Exponential:
   Instead of time-stepping, the transient solution $u(t)$ is expressed as the action of the matrix exponential $\exp(-t M^{-1}K)$ on an initial vector. The authors employ a rational near-best approximation (RBA) of this exponential, constructed using the RKFIT algorithm. This approximation represents the solution as a linear combination of shifted linear systems:
   $$u(t_j) \approx 2 \text{Re} \left( \sum_{i=1}^m \alpha_i^{(j)} (K - \xi_i M)^{-1} f \right)$$
   where $\{ \xi_i \}$ are poles shared across all time points, and $\{ \alpha_i^{(j)} \}$ are precomputed coefficients. Crucially, the shifted system matrices $(K - \xi_i M)$ are independent of time.

2. Parallel Forward and Sensitivity Computation:
   The reformulation allows the $m$ shifted linear systems to be solved independently and in parallel. The time dependence is entirely contained in the residuals (coefficients $\alpha$), meaning the computation of forward responses and sensitivities is effectively independent of the number of desired observation times $K$.

   • Forward Solver: Uses parallel direct solvers (MUMPS) to factorize and solve the shifted systems.
   • Sensitivities (Jacobian): The Jacobian is computed implicitly using the same set of shifted systems. The derivative of the mass matrix is applied to the precomputed solution vectors $g_i = (K - \xi_i M)^{-1}f$. This avoids the sequential forward/adjoint recursions typical of time-stepping methods.

3. Inversion Framework:

   • Optimization: The inverse problem is solved using a Gauss–Newton framework. The linearized least-squares problem is solved via the LSQR algorithm, which handles the Jacobian as an implicit operator (matrix-free).
   • Regularization: Smoothness constraints are formulated using Raviart–Thomas (RT0) elements to handle piecewise constant model parameters on unstructured meshes, ensuring a consistent discretization of the gradient.
   • Implementation: The code is written in Julia, utilizing the Gridap.jl framework for finite-element assembly and a custom MPI interface for MUMPS to manage parallel factorizations on shared-memory architectures.

Key Contributions

• Paradigm Shift in Algorithmic Design: The paper demonstrates a shift from adapting sequential time-marching methods to parallel hardware, to designing mathematical formulations (RBA) that inherently expose concurrency.
• Decoupling of Time Channels: The computational cost is determined by the number of poles $m$ (typically small, e.g., 16–21) rather than the number of time channels $K$ (often hundreds). This allows for high-resolution time sampling without increasing the linear algebra workload.
• Scalable Inversion Workflow: By embedding parallelism at the level of the mathematical model, the authors achieve scalable inversion where the dominant cost (solving shifted systems) can be distributed across MPI ranks.
• Open-Source Release: The implementation is released as open-source software to support reproducible research in the geophysical community.

Results
The method was validated on a synthetic 3-D model with approximately 700,000 degrees of freedom, featuring conductive and resistive anomalies in a half-space.

• Reconstruction Accuracy: The inversion successfully recovered the locations and general shapes of the anomalies. Conductive blocks were well-resolved, while resistive blocks showed reduced contrast, consistent with the known limitations of TEM sensitivity and smoothness regularization.
• Performance:
  • Scaling: Strong scaling experiments demonstrated near-ideal performance up to 8 MPI ranks for both factorization and back-substitution. Efficiency degraded slightly at 16 ranks due to shared-memory NUMA (Non-Uniform Memory Access) effects and memory bandwidth contention.
  • Memory: The peak memory consumption was approximately 620 GB, required to store the factorizations of 21 shifted systems simultaneously.
  • Comparison: Compared to implicit Euler time-stepping (which requires 31 real-valued factorizations for the same time window), the RBA approach requires 21 complex-valued factorizations. While complex arithmetic is more expensive, the RBA method eliminates the sequential time-marching critical path, offering superior potential for wall-clock time reduction on distributed architectures.

Significance and Claims
The authors claim that this work provides a scalable alternative to standard 3-D TEM inversion, specifically addressing the bottleneck of temporal recursion. The primary advantage is not necessarily a reduction in serial arithmetic work, but the reduction of the parallel critical path and the independence of cost from the number of time channels.

The paper acknowledges current limitations, noting that the simultaneous storage of multiple sparse factorizations creates a significant memory footprint and that shared-memory architectures face bandwidth bottlenecks. The authors suggest that future work should focus on transitioning to distributed high-performance computing (HPC) clusters and potentially replacing direct solvers with preconditioned iterative solvers to further improve scalability and reduce memory usage. The study concludes that the proposed framework successfully recovers 3-D structures and establishes a foundation for more efficient, large-scale electromagnetic inversions.