Neural Eikonal Solvers

• Neural eikonal solvers are deep learning methods that compute first-arrival travel times in heterogeneous media, offering mesh-free and highly parallelizable solutions.
• They employ advanced architectures such as PINNs, DeepONet, and equivariant neural fields to enforce physics-based loss functions and improve solution stability.
• These techniques enable rapid seismic imaging, velocity inversion, and tomography by outperforming traditional methods in accuracy, scalability, and generalization.

Neural eikonal solvers comprise a class of machine learning-based methodologies for computing solutions to the eikonal equation, a non-linear Hamilton–Jacobi PDE that models first-arrival travel times in heterogeneous media. These solvers replace or augment classical numerical schemes (e.g., fast marching/sweeping) with deep neural architectures—typically physics-informed neural networks (PINNs), neural operators, or hybrid schemes—yielding mesh-free, highly parallelizable, and operator-learned surrogates for geodesic traveltimes, seismic imaging, tomography, and geometric computing. Recent directions include equivariant neural fields for grid-free, geometry-adaptive generalization, robust factorizations for caustics, and operator-based approaches for inversion and multi-parameter tomographic mapping.

1. Mathematical Foundations and Problem Formulations

The canonical steady-state eikonal equation in a domain $\Omega\subset\mathbb{R}^d$ (typically $d=2,3$) is

$$|\nabla T(x)| = \frac{1}{v(x)}, \qquad x \in \Omega,$$

with a boundary or initial condition $T(x_0) = 0$ at a source $x_0$. Here, $v(x)$ denotes the velocity field, and $T(x)$ is the first-arrival travel time. In seismic or elasticity contexts, multiparameter systems for P- and S-waves appear: $$|\nabla T_p(x)|^2 = \frac{1}{v_p^2(x)},\qquad |\nabla T_s(x)|^2 = \frac{1}{v_s^2(x)}$$ for P- and S-wave speeds $v_p, v_s$.

Variations in neural approaches include:

• One-point eikonal: $T(x; x_0)$ for fixed $d=2,3$0 (pinned source).
• Two-point eikonal: $d=2,3$1, symmetric in $d=2,3$2, $d=2,3$3 (arbitrary source-receiver).
• Inverse mapping: $d=2,3$4 (velocity inversion).

Mesh-based formulations are replaced by mesh-free surrogates, often leveraging factorization to absorb source-based singularities via

$d=2,3$5

where $d=2,3$6.

2. Neural Architectures for Eikonal Solving

Neural eikonal solvers adopt a range of architectures, with design dictated by the desired mapping, problem geometry, and properties of the target solution.

Physics-Informed Neural Networks (PINNs)

PINNs treat $d=2,3$7 or a related factor (e.g., slowness or normalized time) as the output of a neural network, trained by minimizing a physics-based residual: $d=2,3$8 where $d=2,3$9. For two-point formulations, symmetry-enforced architectures (e.g., $$|\nabla T(x)| = \frac{1}{v(x)}, \qquad x \in \Omega,$$0) improve reciprocity.

Deep Operator Networks (DeepONet)

DeepONets encode operator mappings from function spaces (e.g., $$|\nabla T(x)| = \frac{1}{v(x)}, \qquad x \in \Omega,$$1 or $$|\nabla T(x)| = \frac{1}{v(x)}, \qquad x \in \Omega,$$2) to function-valued outputs using a "branch" network for the input function and a "trunk" network for the coordinate query, merged via an inner product. For example, forward DeepONet receives $$|\nabla T(x)| = \frac{1}{v(x)}, \qquad x \in \Omega,$$3 and $$|\nabla T(x)| = \frac{1}{v(x)}, \qquad x \in \Omega,$$4 as input, predicts $$|\nabla T(x)| = \frac{1}{v(x)}, \qquad x \in \Omega,$$5; inverse DeepONet receives $$|\nabla T(x)| = \frac{1}{v(x)}, \qquad x \in \Omega,$$6 as input, predicts $$|\nabla T(x)| = \frac{1}{v(x)}, \qquad x \in \Omega,$$7 (Mei et al., 2023).

Local Neural Solvers in Fast Marching/Sweeping

The local update in classical fast marching (FMM) may be replaced with a compact feedforward network trained to regress local distances from neighbor values, retaining the O($$|\nabla T(x)| = \frac{1}{v(x)}, \qquad x \in \Omega,$$8) global ordering (Lichtenstein et al., 2019).

Equivariant Neural Fields (ENFs)

ENF-based neural eikonal solvers model $$|\nabla T(x)| = \frac{1}{v(x)}, \qquad x \in \Omega,$$9, encoding families of solutions parameterized by latent context $T(x_0) = 0$0, structured via group actions for equivariance under Euclidean or more general homogeneous manifold transforms. This enables grid-free, geometry-adaptive learning with explicit steerability (García-Castellanos et al., 21 May 2025).

Table: Representative Network Choices

Architecture	Input	Output
PINN	$T(x_0) = 0$1 (and/or $T(x_0) = 0$2, $T(x_0) = 0$3)	$T(x_0) = 0$4, $T(x_0) = 0$5
DeepONet	$T(x_0) = 0$6 (branch), $T(x_0) = 0$7, $T(x_0) = 0$8	$T(x_0) = 0$9
Local solver (FMM)	$x_0$0	$x_0$1
ENF	$x_0$2, $x_0$3, $x_0$4	$x_0$5

Architectural advances involve residual blocks, locally adaptive activations (e.g., arctan with trainable slope), bounded activations (to constrain the physical range), and group-invariant pooling for equivariant representation (García-Castellanos et al., 21 May 2025, Waheed et al., 2020, Grubas et al., 2022).

3. Loss Functions and Training Strategies

Loss design is dictated by the need to enforce the eikonal PDE, boundary conditions, and, for some approaches, data consistency or operator mapping accuracy. Key loss formulations include:

• PDE Residual Minimization: $x_0$6 penalized (L1 or L2).
• Factored Residuals: Using analytically known factors to regularize singularities and improve convergence, e.g., via $x_0$7 as in PINNPStomo's $x_0$8 (Song et al., 2024).
• Hamiltonian Loss: L1 norm of the Hamiltonian $x_0$9 provides robustness near caustics and accommodates non-smooth behavior (Grubas et al., 2022).
• Data-informed Losses: For inverse problems or when data are available, MSE to observed $v(x)$0.
• Regularization: Positivity/physical range penalties, boundary enforcement, and (in DeepONet) potentially normalization—but explicit regularization may be omitted in high-performing models (Mei et al., 2023).

Adaptive optimizers (Adam, L-BFGS-B), dynamic gradient-weight balancing, two-stage optimization (coarse then fine collocation), and data-driven sampling (weighted by error concentration) are common strategies for accelerating or stabilizing training (Waheed et al., 2020, Waheed et al., 2021, Grubas et al., 2022, Song et al., 2024).

4. Advances in Factorization, Robustness, and Generalization

Factorization techniques address the $v(x)$1 singularity at sources and improve learning dynamics:

• Classical Multiplicative Factorization: $v(x)$2, $v(x)$3 reflecting the homogeneous solution (Waheed et al., 2020, Waheed et al., 2021).
• New Factored Forms: PINNPStomo employs $v(x)$4, $v(x)$5, making the factorization source-velocity independent and eliminating the need for a background velocity $v(x)$6 (Song et al., 2024).
• Bounding and Symmetrization: Use of sigmoid mappings to constrain output range, reciprocity averaging for two-point networks (Grubas et al., 2022).

For challenging regimes with caustics (multi-valued solutions, non-differentiable wavefronts), tailored architectures (Gaussian activation, L1 Hamiltonian loss), input scaling, and explicit factorization enable highly accurate and stable learning, with performance superior to standard PINNs or FMMs (e.g., relative MAE $v(x)$7 0.2–0.8% vs. FMM on Marmousi) (Grubas et al., 2022).

Equivariance and steerability, as in E-NES, allow for transferable neural representations across rotated, translated, or otherwise transformed domains, and enable user-controllable solution fields via latent group actions (García-Castellanos et al., 21 May 2025).

5. Operator Learning and Inverse Problems

Operator learning frameworks generalize the solution concept beyond fixed initial/boundary conditions:

• Forward Operator Learning: $v(x)$8; DeepONet models can generate travel times for arbitrary new sources without retraining, provided the velocity model class matches training (Mei et al., 2023).
• Inverse Operator Learning: $v(x)$9; inverse DeepONet infers underlying velocity models from measured first-arrival times, enabling parametric seismic inversion (Mei et al., 2023).
• Multi-parameter Inversion: PINNPStomo uses coupled networks and a joint loss for simultaneous recovery of $T(x)$0, achieving multiparameter inversion in both 2D and 3D settings (Song et al., 2024).
• Transfer and Surrogate Modeling: PINN surrogates trained on multiple source positions enable rapid travel time evaluation for new sources, with transfer learning accelerating retraining under velocity perturbations (Waheed et al., 2020, Waheed et al., 2021).

6. Computational Performance and Benchmarking

Neural eikonal solvers are evaluated versus traditional schemes across accuracy, speed, memory, and scalability:

• Accuracy: Neural approaches consistently achieve lower absolute and relative travel time errors compared to first- and second-order fast marching/sweeping, and handle challenging regimes (sharp contrasts, caustics) where traditional methods degrade (Lichtenstein et al., 2019, Grubas et al., 2022).
• Speed and Scalability: Training can be computationally intensive, but once trained, inference is grid-free and massively parallel: e.g., EikoNet computes $T(x)$1 queries in 0.424s on a V100 GPU, massively undercutting FMM lookup table requirements (Smith et al., 2020).
• Generalization: Approaches such as DeepONet and ENF provide generalization to new sources/receivers, velocity fields, or domain geometries, while ENF generalizes to non-Euclidean manifolds (spheres, hyperbolic spaces) with group equivariance (García-Castellanos et al., 21 May 2025, Mei et al., 2023).
• Robustness: Gaussian activations and robust L1 Hamiltonian losses for NES increase reliability near caustics and improve training speed and error by 10–60$T(x)$2 over previous neural solvers (Grubas et al., 2022).
• Operator Surrogacy: Neural surrogates admit O(1) inference in travel-time tables, advantageous in iterative inversion or real-time localization (Waheed et al., 2021, Waheed et al., 2020).

7. Applications, Extensions, and Open Directions

Neural eikonal solvers are deployed in seismic forward modeling, velocity inversion, tomography, hypocenter localization, and geometric computing. Noteworthy features and future directions include:

• Mesh-free flexibility: All leading approaches are mesh-free, enabling seamless handling of irregular topographies, free surfaces, and unstructured domains (Waheed et al., 2020, Waheed et al., 2021, Song et al., 2024).
• Seismic tomography and multiparameter inversion: PINNPStomo establishes simultaneous P/S inversion from joint travel times with no background speed assumption (Song et al., 2024).
• Scalability to 3D/complex physics: References extend PINN frameworks to anisotropy, attenuation, and sharp heterogeneity with only residual modification (Waheed et al., 2020, Waheed et al., 2021), while E-NES and PINNPStomo demonstrate 3D scalability and data-driven training (García-Castellanos et al., 21 May 2025, Song et al., 2024).
• Operator learning open questions: Stability, uniqueness, and robustness, especially in inverse problems, as well as the efficacy of physics-informed loss terms (vs. data-driven), remain areas of investigation (Mei et al., 2023).
• Meta-learning and adaptive sampling: Unsolved issues involve convergence in highly heterogeneous fields (spectral bias), which may be addressed via active collocation, meta-learning, or hybrid PINN/operator schemes (Waheed et al., 2020, Waheed et al., 2021).

A plausible implication is that future neural eikonal solvers will integrate operator learning, geometric equivariance, and robust physics-informed regularization to achieve real-time, high-fidelity, and physically interpretable solutions across a spectrum of domains and acquisition geometries.