pith. sign in

arxiv: 2605.19951 · v1 · pith:LUAK2ZKAnew · submitted 2026-05-19 · 🧮 math.NA · cs.NA· physics.geo-ph

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

Ralph-Uwe B\"orner , Stefan G\"uttel , Thomas G\"unther This is my paper

Pith reviewed 2026-05-20 04:09 UTC · model grok-4.3

classification 🧮 math.NA cs.NAphysics.geo-ph
keywords transient electromagnetic inversion3D inversionmatrix exponentialrational approximationGauss-Newton methodparallel direct solverselectromagnetic modeling
0
0 comments X

The pith

Rational near-best approximation of the matrix exponential isolates time dependence in the residuals so that forward responses and sensitivities in 3D transient electromagnetic inversion become independent of the number of observation times

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper presents a parallel Gauss-Newton method for large-scale three-dimensional transient electromagnetic inversion that replaces repeated time-stepping with a rational near-best approximation to the matrix exponential. The approximation is expressed through partial fractions, converting the problem into a collection of shifted linear systems whose solutions are combined using time-dependent residuals. Because the time dependence now lives only in those residuals, the cost of computing forward responses and sensitivities no longer grows with the number of observation times. The method is implemented with parallel direct solvers on shared-memory hardware, smoothness regularization on Raviart-Thomas elements, and an implicit Jacobian solved by LSQR. Experiments recover a synthetic conductivity model with roughly 700,000 degrees of freedom, demonstrating that the approach makes dense time sampling computationally affordable.

Core claim

By representing the matrix exponential via a rational near-best approximation and its partial-fraction expansion, the time dependence of the electromagnetic fields is moved entirely into the residuals of the rational functions; consequently the forward responses and sensitivities required by the Gauss-Newton iteration become independent of the number of desired observation times.

What carries the argument

Rational near-best approximation of the matrix exponential expressed as a partial-fraction sum, reducing the transient problem to a set of shifted linear systems solved once and then combined with time-dependent residual factors.

If this is right

  • Forward modeling and sensitivity calculations scale independently of the number of observation times.
  • Parallel direct solvers for the shifted systems enable efficient shared-memory execution via MPI.
  • Smoothness regularization formulated with Raviart-Thomas elements integrates directly into the Gauss-Newton framework.
  • Synthetic conductivity structures with approximately 700,000 degrees of freedom are recoverable without time-stepping cost.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same residual-based decoupling could be applied to other time-domain geophysical inverse problems that involve matrix exponentials.
  • Replacing the direct solvers with preconditioned iterative methods would reduce the memory footprint caused by storing multiple sparse factorizations simultaneously.
  • The method would permit inversion at finer time discretizations to resolve fast transients without a proportional rise in runtime.

Load-bearing premise

The rational approximation remains accurate enough across the conductivity models and time windows that arise during the nonlinear inversion iterations and does not introduce systematic bias into the recovered conductivity structure.

What would settle it

Run the same synthetic 3-D data set through both the rational-approximation inversion and a conventional time-stepping forward solver; any systematic difference in the recovered conductivity features that grows with the number of time channels would falsify the claim of unbiased independence.

Figures

Figures reproduced from arXiv: 2605.19951 by Ralph-Uwe B\"orner, Stefan G\"uttel, Thomas G\"unther.

Figure 1
Figure 1. Figure 1: Schematic illustration of the parallel distribution of shifted linear systems arising from the rational approximation. The system matrices K, M, the poles ξi , and the right-hand side vector f are distributed across computational resources. Each worker is assigned one or more poles ξi , forming, factorizing, and solving the corresponding linear systems (K − ξiM)gi = f independently. The resulting solutions… view at source ↗
Figure 2
Figure 2. Figure 2: Taylor remainder test for the forward operator and its linearization. Red markers show the norm of the first-order difference e0(h) (25), while the blue markers show the norm of the second-order remainder e1(h) (22). Solid lines indicate the expected asymptotic behaviour of order O(h) and O(h 2 ), respectively. the relative mismatch [PITH_FULL_IMAGE:figures/full_fig_p010_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Three-dimensional view of the synthetic model and survey configuration. The square transmitter loop is located at the surface, together with the receiver grid, while two conductive (red) and two resistive (blue) blocks are embedded in a homogeneous half-space at shallow depth. The uniform background is omitted for clarity. receiver array consists of 49 surface stations arranged in a regular 7 × 7 grid with… view at source ↗
Figure 4
Figure 4. Figure 4: Three-dimensional view of the reconstructed conductivity model obtained from the synthetic data after 25 iterations. The inversion recovers the locations and general shapes of the conductive blocks embedded in the resistive half-space. The anomalies appear slightly smoothed and reduced in amplitude due to regularization effects. behaviour is expected due to the applied smoothness regularization and the lim… view at source ↗
Figure 5
Figure 5. Figure 5: Slice through the reconstructed conductivity model obtained after 25 Gauss–Newton iterations. The slice illustrates the lateral extent and depth positioning of the recovered anomalies. While the main structures are clearly identified, the anomalies are moderately smoothed and exhibit reduced contrast, reflecting the influence of regularization and the limited resolution of the data. ceivers and time channe… view at source ↗
Figure 6
Figure 6. Figure 6 [PITH_FULL_IMAGE:figures/full_fig_p017_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Heat map of the weighted data residuals Wd(d(m) − d obs) after 25 Gauss–Newton iterations. Rows correspond to time channels, and columns correspond to receiver locations arranged according to the acquisition grid. The colour scale is symmetric about zero, highlighting positive and negative weighted deviations of the model response from the observed data. strates that the computational cost is dominated by … view at source ↗
Figure 8
Figure 8. Figure 8: Inversion convergence and performance. Bars show the computational time required for each Gauss–Newton (GN) iteration (left axis), while the red curve indicates the corresponding evolution of the data misfit χ 2 (right axis). Numerical annotations within the bars denote the number of LSQR iterations performed at each Gauss–Newton step. The inset displays the evolution of the regularization parameter λ thro… view at source ↗
Figure 9
Figure 9. Figure 9: Strong-scaling behaviour of the shifted linear systems for an RBA approximation with 16 poles. The left panels show the measured runtimes for the sparse factorization and subsequent back-substitution stages, while the right panel shows the corresponding parallel efficiency relative to the single-rank run [PITH_FULL_IMAGE:figures/full_fig_p019_9.png] view at source ↗
Figure 10
Figure 10. Figure 10: Measured and predicted runtimes per Gauss–Newton iteration for an RBA approximation with 21 poles. The predicted runtimes are obtained from a linear performance model consisting of a fixed factorization cost and a contribution proportional to the number of LSQR iterations. Numerical labels indicate the corresponding number of LSQR iterations for each Gauss–Newton step. implementation, the MUMPS factors ar… view at source ↗
read the original abstract

We present a novel parallel implementation for large-scale three-dimensional electromagnetic inversion based on a Gauss-Newton framework combined with a rational near-best approximation of the matrix exponential for transient simulations. The method employs parallel direct solvers for the shifted linear systems arising from the partial fraction representation of the rational approximation and demonstrates efficient parallel execution on a shared-memory architecture using MPI. A key property of the approach is that the time dependence is entirely contained in the residuals of the employed rational functions, such that the computation of forward responses and sensitivities becomes effectively independent of the number of desired observation times. Model regularization is done with smoothness constraints, formulated with Raviart-Thomas elements. The linearized inverse problems are solved using LSQR, using an implicit parallel Jacobian operator. Numerical experiments demonstrate the successful recovery of a synthetic 3-D conductivity structure with approximately 700,000 degrees of freedom. The study further discusses computational bottlenecks related to memory consumption and shared-memory scalability arising from the simultaneous storage of multiple sparse matrix factorizations. Possible improvements based on preconditioned iterative solvers and distributed high-performance computing architectures are outlined. The implementation in the Julia programming language is released as open-source software to support reproducible research and further development by the geophysical inversion community.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 2 minor

Summary. The manuscript presents a parallel implementation of 3-D transient electromagnetic (TEM) inversion within a Gauss-Newton framework. Forward responses and sensitivities are computed via a rational near-best approximation to the matrix exponential, with time dependence isolated in the rational-function residuals so that the cost becomes independent of the number of observation times. Shifted linear systems are solved with parallel direct factorizations; the inverse problem is solved with LSQR using an implicit Jacobian operator; regularization uses smoothness constraints discretized with Raviart-Thomas elements. Numerical experiments recover a synthetic conductivity model with roughly 700 000 degrees of freedom, and the Julia code is released as open source.

Significance. If the rational approximation remains uniformly accurate across the sequence of conductivity models visited during nonlinear iterations, the method would provide a practical route to large-scale 3-D TEM inversion whose computational cost scales with the number of spatial degrees of freedom rather than the number of time samples. The open-source release and emphasis on shared-memory parallelism are positive contributions to reproducible geophysical inversion research.

major comments (2)
  1. [Numerical experiments] Numerical experiments section: the synthetic recovery is shown but no quantitative approximation-error norms (e.g., ||exp(A) - r(A)|| for the conductivity models at each Gauss-Newton iteration) or direct comparison against a reference time-stepping scheme are reported. Without these metrics it is impossible to verify that approximation errors do not systematically bias the recovered conductivity structure.
  2. [Rational approximation] Rational-approximation construction (likely §3): the poles and residues are chosen once (near-best approximation) and then held fixed while the system matrix A(σ) changes with each conductivity update. No a-priori error bound or numerical study is supplied that quantifies how the approximation error grows when σ varies over the range encountered in the inversion loop.
minor comments (2)
  1. The description of the implicit parallel Jacobian operator would benefit from an explicit algorithmic outline or pseudocode showing how the action on a vector is assembled from the already-factorized shifted systems.
  2. Memory-consumption discussion in the conclusions could be strengthened by reporting the actual peak RAM usage and the number of simultaneous factorizations stored for the largest run.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the constructive report and the positive assessment of the parallel implementation and open-source release. We address the two major comments below, agreeing that additional quantitative validation will strengthen the manuscript.

read point-by-point responses

  1. Referee: Numerical experiments section: the synthetic recovery is shown but no quantitative approximation-error norms (e.g., ||exp(A) - r(A)|| for the conductivity models at each Gauss-Newton iteration) or direct comparison against a reference time-stepping scheme are reported. Without these metrics it is impossible to verify that approximation errors do not systematically bias the recovered conductivity structure.

    Authors: We agree that the current numerical experiments would benefit from explicit quantitative error metrics. In the revised manuscript we will add a dedicated subsection reporting the matrix-exponential approximation error ||exp(A(σ)) - r(A(σ))|| evaluated at the conductivity models obtained after each Gauss-Newton iteration. We will also include a direct comparison of forward responses computed with the rational approximation against a reference time-stepping scheme on a smaller-scale test problem, thereby confirming that the approximation errors remain small enough not to introduce systematic bias in the recovered conductivity. revision: yes

  2. Referee: Rational-approximation construction (likely §3): the poles and residues are chosen once (near-best approximation) and then held fixed while the system matrix A(σ) changes with each conductivity update. No a-priori error bound or numerical study is supplied that quantifies how the approximation error grows when σ varies over the range encountered in the inversion loop.

    Authors: The rational approximant was constructed once as a near-best approximation over a spectral interval chosen to encompass the eigenvalues expected for the conductivity range explored by the inversion. While the manuscript does not contain an a-priori error bound, we will add a numerical study in the revised version that evaluates the approximation error for a sequence of conductivity models sampled from the inversion trajectory. This study will quantify the variation of the error with σ and demonstrate that it remains negligible relative to the discretization and regularization errors. revision: yes

Circularity Check

0 steps flagged

No circularity: time-independence follows from standard rational approximation theory

full rationale

The paper derives the claimed time-independence of forward responses and sensitivities directly from the partial-fraction decomposition of a rational approximant to the matrix exponential. This property is a standard consequence of the approximation (time factors appear only in the scalar residuals) and does not reduce to any fitted parameter, self-defined quantity, or load-bearing self-citation within the present work. The implementation uses established parallel direct solvers, LSQR, and Raviart-Thomas regularization; the synthetic experiment validates feasibility but supplies no definitional input to the core claim. The derivation chain therefore remains self-contained against external mathematical results.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

The central claim rests on the accuracy of a rational approximation drawn from existing approximation theory and on the correctness of parallel direct solvers for the shifted systems; no new free parameters are introduced in the abstract description.

axioms (1)
  • domain assumption Rational near-best approximations can represent the action of the matrix exponential for the discretized Maxwell operator to sufficient accuracy for inversion purposes.
    Invoked when stating that time dependence factors entirely into the rational residuals.

pith-pipeline@v0.9.0 · 5755 in / 1261 out tokens · 47639 ms · 2026-05-20T04:09:46.881011+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

14 extracted references · 14 canonical work pages

  1. [1]

    B¨orner, R.-U

    Julia: A fresh approach to numerical computing,SIAM Review,59(1), 65–98, doi: 10.1137/141000671. B¨orner, R.-U. & G ¨uttel, S.,

    work page doi:10.1137/141000671

  2. [2]

    B¨orner, R.-U., Ernst, O

    Fast parallel transient electromagnetic modelling using a uniform- in-time approximation to the exponential,Geophysical Journal International,243(1), ggaf319, doi: 10.1093/gji/ggaf319. B¨orner, R.-U., Ernst, O. G., & G ¨uttel, S.,

    work page doi:10.1093/gji/ggaf319

  3. [3]

    Cai, H., Liu, M., Zhou, J., Li, J., & Hu, X.,

    Three-dimensional transient electromagnetic mod- elling using Rational Krylov methods,Geophysical Journal International,202(3), 2025–2043, doi: 10.1093/gji/ggv224. Cai, H., Liu, M., Zhou, J., Li, J., & Hu, X.,

    work page doi:10.1093/gji/ggv224 2025

  4. [4]

    Danisch, S

    SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications,Computers and Geosciences,85, 142–154, doi: 10.1016/j.cageo.2015.09.015. Danisch, S. & Krumbiegel, J.,

    work page doi:10.1016/j.cageo.2015.09.015 2015

  5. [6]

    Han, B., Li, Y ., & Li, G.,

    Inversion of time domain three-dimensional electromagnetic data,Geophysical Journal International,171(2), 550–564, doi: 10.1111/j.1365- 246x.2007.03365.x. Han, B., Li, Y ., & Li, G.,

    work page doi:10.1111/j.1365- 2007

  6. [7]

    Heagy, L

    3D forward modeling of magnetotelluric fields in general anisotropic media and its numerical implementation in Julia,Geophysics,83(4), F29–F40, doi: 10.1190/geo2017-0515.1. Heagy, L. J., Cockett, R., Kang, S., Rosenkjaer, G. K., & Oldenburg, D. W.,

    work page doi:10.1190/geo2017-0515.1

  7. [8]

    Lee, E.-J., Huang, H., Dennis, J

    A frame- work for simulation and inversion in electromagnetics,Computers and Geosciences,107, 1–19, doi: 10.1016/j.cageo.2017.06.018. Lee, E.-J., Huang, H., Dennis, J. M., Chen, P., & Wang, L.,

    work page doi:10.1016/j.cageo.2017.06.018 2017

  8. [9]

    Liu, S., Liu, R., Zhao, H., Sun, H., & Liu, D.,

    An optimized parallel LSQR algorithm for seismic tomography,Computers and Geosciences,61, 184–197, doi: 10.1016/j.cageo.2013.08.013. Liu, S., Liu, R., Zhao, H., Sun, H., & Liu, D.,

    work page doi:10.1016/j.cageo.2013.08.013 2013

  9. [10]

    Nocedal, J

    3D parallel inversion of time-domain airborne EM data,Applied Geophysics,13(4), 701–711. Nocedal, J. & Wright, S. J., 2006.Numerical Optimization, Springer, New York, 2nd edition., doi: 10.1007/978-0-387-40065-5. Paige, C. C. & Saunders, M. A.,

    work page doi:10.1007/978-0-387-40065-5 2006

  10. [11]

    Rochlitz, R., Seidel, M., & B ¨orner, R.-U.,

    custEM: Customizable finite-element simulation of complex controlled-source electromagnetic data,GEOPHYSICS,84(2), F17–F33, doi: 10.1190/geo2018-0208.1. Rochlitz, R., Seidel, M., & B ¨orner, R.-U.,

    work page doi:10.1190/geo2018-0208.1

  11. [12]

    Ruthotto, L., Treister, E., & Haber, E.,

    Evaluation of three approaches for simulating 3-D time-domain electromagnetic data,Geophysical Journal International,227(3), 1980–1995. Ruthotto, L., Treister, E., & Haber, E.,

    work page 1980

  12. [13]

    Weiss, M., Rochlitz, R., & G ¨unther, T.,

    Electromagnetic Modeling Using Adaptive Grids – Error Estimation and Geometry Representation,Surveys in Geophysics,45(1), 277–314, doi: 10.1007/s10712-023-09794-9. Weiss, M., Rochlitz, R., & G ¨unther, T.,

    work page doi:10.1007/s10712-023-09794-9

  13. [14]

    Werthm¨uller, D., Rochlitz, R., Castillo-Reyes, O., & Heagy, L.,

    Evaluation of an iterative framework for geophysical electromagnetic forward and inverse modelling problems,Geophysical Journal International,243(1), doi: 10.1093/gji/ggaf290. Werthm¨uller, D., Rochlitz, R., Castillo-Reyes, O., & Heagy, L.,

    work page doi:10.1093/gji/ggaf290

  14. [15]

    Towards an open-source landscape for 3-D CSEM modelling,Geophysical Journal International,227(1), 644–659, doi: 10.1093/gji/ggab238

    work page doi:10.1093/gji/ggab238