Factored eikonal: removing the source singularity

Part 7 — Travel-time, surface-wave, and joint inversion

Learning objectives

Derive the factored form T(x) = T₀(x)·τ(x) and explain why it cures the source singularity
Re-derive the eikonal residual analytically as a polynomial in (τ, ∇τ)
Train a factored EikoNet end-to-end and quantify the order-of-magnitude improvement vs §7.2
Recognise why the factored form also cures the trivial-solution trap structurally
See where factored EikoNet still struggles (caustics, multi-pathing) — pointing to multi-branch neural eikonal solvers (Grubas 2023)

§7.2 hit a hard plateau at ~5% relative travel-time error because it asked a smooth Tanh MLP to represent a function with diverging second derivatives near the source. The factored eikonal of Treister-Haber 2016 / Fomel-Luo-Zhao 2009 / bin Waheed et al 2021 PINNeik dissolves this difficulty in one move. This section derives the formulation, races a factored EikoNet against the §7.1 FSM reference on the SAME setup as §7.2, and quantifies the improvement.

The factored form

Write the travel-time field as the product of an analytic homogeneous-medium part and a smooth multiplicative correction:

T(x) = T_0(x) \cdot \tau(x), \quad T_0(x) = \|x - x_{\mathrm{src}}\| / c_{\mathrm{src}} ,

where $c_{\mathrm{src}} = c(x_{\mathrm{src}})$ is the velocity at the source. $T_0$ is the EXACT eikonal solution if $c$ were uniform everywhere at $c_{\mathrm{src}}$ . The correction $\tau$ is the factor by which the true bent-ray travel time differs from the straight-ray homogeneous answer. For a smooth velocity field, $\tau$ is a smooth, bounded function near 1 — exactly the kind of function a small Tanh MLP excels at.

Three structural properties make this parameterisation magical:

Source BC is automatic. $T(x_{\mathrm{src}}) = T_0(x_{\mathrm{src}}) \cdot \tau(x_{\mathrm{src}}) = 0 \cdot \tau(x_{\mathrm{src}}) = 0$ for ANY $\tau$ . We do not need a BC penalty term in the loss.
1/r source singularity is absorbed. The diverging second derivatives of $T$ near the source live entirely in $T_0$ , which is computed analytically. The network never sees them.
Trivial-solution trap is structurally cured. $|\nabla T| \equiv 0$ requires $\tau \equiv 0$ AND $\nabla \tau \equiv 0$ — a far stronger collapse than $T \equiv 0$ . With $\tau$ initialised near 1, the optimiser starts a long way from any trivial-zero stationary point.

The eikonal residual in factored form

Apply the product rule to compute the gradient of $T = T_0 \tau$ :

\nabla T = \tau \, \nabla T_0 + T_0 \, \nabla \tau .

For the homogeneous $T_0 = \rho / c_{\mathrm{src}}$ with $\rho = |x - x_{\mathrm{src}}|$ :

\nabla T_0 = \frac{x - x_{\mathrm{src}}}{\rho \, c_{\mathrm{src}}}, \qquad |\nabla T_0|^2 = \frac{1}{c_{\mathrm{src}}^2} .

Substituting and simplifying:

|\nabla T|^2 = \frac{1}{c_{\mathrm{src}}^2} \Bigl[ \tau^2 + 2\tau\,\bigl( (x - x_{\mathrm{src}}) \cdot \nabla \tau \bigr) + \rho^2 \, |\nabla \tau|^2 \Bigr] .

The factor $1/c_{\mathrm{src}}^2$ pulls cleanly out front because the cross term $2 \tau T_0 (\nabla T_0 \cdot \nabla \tau) = 2\tau,(\rho/c_{\mathrm{src}}),((x - x_{\mathrm{src}}) \cdot \nabla \tau)/(\rho, c_{\mathrm{src}})$ also has a $1/c_{\mathrm{src}}^2$ factor.

The full eikonal residual is then

r(x) = c(x)^2 \, |\nabla T(x)|^2 - 1 = \frac{c(x)^2}{c_{\mathrm{src}}^2} \Bigl[ \tau^2 + 2\tau ((x - x_{\mathrm{src}}) \cdot \nabla \tau) + \rho^2 |\nabla \tau|^2 \Bigr] - 1 .

This is a POLYNOMIAL in $(\tau, \nabla \tau)$ — no special handling needed, no division by $\rho$ . The widget computes its squared mean over collocation points as the entire training loss.

Network parameterisation

The 2-32-32-1 Tanh MLP outputs a scalar $\delta\tau(x; \theta)$ and we define

\tau(x) = \exp(\delta\tau(x; \theta)) ,

so the random small-weight init gives $\tau \approx 1$ everywhere AND $\tau > 0$ ALWAYS. The exponential is a critical detail: the eikonal residual is invariant under $T \to -T$ (since $|\nabla T|^2 = |-\nabla T|^2$ ), and the factored form inherits this as $\tau \to -\tau$ . With residual values up to $\sim 7$ at init in our gradient model (the deep layers have $c^2/c_{\mathrm{src}}^2 \approx 3.6$ ), an Adam step that crosses $\tau = 0$ would lock the optimiser into the wrong basin where $T < 0$ . Parameterising $\tau = e^{\delta\tau}$ builds positivity into the function space rather than relying on the optimiser to stay on one side. The chain rule for backprop becomes $\nabla \tau = \tau \cdot \nabla \delta\tau$ and $dL/d\delta\tau = \tau \cdot dL/d\tau + \nabla \delta\tau \cdot dL/d(\nabla \tau) \cdot \tau$ — handled analytically in the widget.

This is the same trick used in production PINN-eikonal codes (PINNeik bin Waheed et al 2021 Section 3.1) and in normalising flows / variational autoencoders to enforce positivity of variance parameters.

Try it

Same domain, velocity model, and source-position controls as §7.2. Click ▶ Train. The widget runs:

FSM reference (~2-5 ms).
Factored EikoNet training: 1500 epochs × 120 collocation points, Adam lr=5e-3 with cosine decay over the last third. NO warmstart, NO BC term. Browser wall-clock: ~10-15 s on a modern laptop — already roughly 2-3× faster than §7.2 because there is no BC penalty term and no warmstart anchor evaluation.

Five panels: T_FSM reference, T_factored = T₀·τ, |T_factored − T_FSM| absolute error map, the LEARNED τ field (plasma colormap to distinguish it from travel-time fields), and the loss trace.

Expected behaviour:

Mean RELATIVE travel-time error $\sim 0.3\text{--}0.8%$ of $T_{\max}$ (vs §7.2's ~5%) — an order-of-magnitude improvement.
τ field range $\sim [0.85, 1.20]$ on this gradient model — gently rising with depth (slower true rays at depth means higher slowness, larger correction factor).
Loss trace: smooth monotone descent from $\sim 1$ to $\sim 10^{-5}$ without the staircase plateaus seen in vanilla EikoNet.
Error map is uniform — no source-neighborhood hot spot. The 1/r singularity is gone.

Why this works (operator-theoretic view)

The factored eikonal is a CHANGE OF VARIABLES that rescales the function space the network has to search. Vanilla EikoNet asks the MLP to span a function space that includes radial cones with diverging Hessians; this is impossible for smooth-activation networks of finite width. Factored EikoNet asks the MLP to span a function space of smooth, bounded multiplicative corrections to a fixed analytic basis function $T_0$ . That smaller, friendlier space IS spanned by a small Tanh MLP — even at low width.

This is the same trick used in CHEBYSHEV POLYNOMIAL approximation for problems with known endpoint behaviour: extract the asymptotic structure analytically, then approximate only the smooth residue. In the eikonal case, the singular asymptotic is the radial $\rho/c_{\mathrm{src}}$ cone; the smooth residue is $\tau$ .

Where factored EikoNet still struggles

The factored form ASSUMES a single-arrival regime: T(x) is the FIRST-arrival travel time and the wavefront from the source spreads monotonically. This breaks down at:

Caustics. Wavefronts fold over themselves; multiple ray paths reach the same point. The viscosity solution still selects first-arrival, but $\tau$ becomes non-smooth at the caustic boundary, and the small MLP smears it.
Velocity inversions / shadow zones. Strong negative velocity gradients can create regions reachable only via reflected rays. First-arrival travel time has discontinuous derivatives at the shadow boundary.
Multi-source regimes. If you want $T$ from multiple sources simultaneously, factor by ONE source biases the network. The cure is to condition the network on $x_{\mathrm{src}}$ as an additional input (Smith-Azizzadenesheli-Ross 2021 do this with $x_{\mathrm{src}}$ concatenated to $x$ ), and use a per-source $T_0$ in the factored form. This source-conditioned EikoNet powers the modern microseismic location codes (HypoSVI, Smith et al 2022) referenced in §7.5.

Grubas-Duchkov-Loginov 2023 introduce a "neural eikonal solver" that handles caustics by training multiple branches; we do not implement it here, but the technique is to split $\tau$ into multi-branch outputs and track which is dominant per region.

Operational note: differentiability through source position

The factored form makes one thing slightly TRICKIER than vanilla EikoNet: differentiating $T(x; x_{\mathrm{src}})$ with respect to $x_{\mathrm{src}}$ for source location problems (§7.5) requires keeping track of the $T_0$ piece. But $T_0 = |x - x_{\mathrm{src}}|/c_{\mathrm{src}}$ is analytic so the source-position derivative is closed-form. Hybrid implementation: compute $\partial T/\partial x_{\mathrm{src}}$ as $\tau \cdot \partial T_0/\partial x_{\mathrm{src}} + T_0 \cdot \partial \tau/\partial x_{\mathrm{src}}$ where the second term comes from the network's dependence on $x_{\mathrm{src}}$ (if conditioned).

References

Treister, E., Haber, E. (2016). A fast marching algorithm for the factored eikonal equation. J. Comput. Phys. 324, 210–225. The factored fast-marching method that established the technique for classical solvers.
Fomel, S., Luo, S., Zhao, H. (2009). Fast sweeping method for the factored eikonal equation. J. Comput. Phys. 228(17), 6440–6455. Earlier factored fast-sweeping; same parameterisation idea.
bin Waheed, U., Haghighat, E., Alkhalifah, T., Song, C., Hao, Q. (2021). PINNeik: Eikonal solution using physics-informed neural networks. Computers & Geosciences 155, 104833. The PINN implementation this widget descends from.
Grubas, S., Duchkov, A., Loginov, G. (2023). Neural eikonal solver: Improving accuracy of physics-informed neural networks for solving eikonal equation in case of caustics. J. Comput. Phys. 474, 111789. Multi-arrival extensions.