Seismic Hessian vector products: a complete Lagrangian formulation and error assessment against automatic differentiation

Summary Hessian–vector products (HVPs) provide matrix-free access to second-order information in full-waveform inversion (FWI), enabling the use of truncated-Newton updates with improved convergence. This may result in enhanced imaging of small-scale scatterers in strongly heterogeneous, high-contrast media. In the time domain, standard adjoint-based HVP derivations are often presented either via perturbation theory or via a Lagrangian (Lagrange-multiplier) framework. Although these routes are formally equivalent, their correctness hinges on a precise treatment of transpose/adjoint placement, boundary/endpoint terms, and the discrete inner product used by the numerical solver. Revisiting these derivations, we show that a commonly used Lagrangian expression for the first HVP contribution can be incorrect when the wave equation operator L is not symmetric, and that rewriting the first term of the HVP, $\mathcal {H}_1$, using L† (or equivalently, swapping the contraction order) requires accounting for the associated boundary terms, where the L† indicates the adjoint of L. To resolve this ambiguity, we provide a complete Lagrangian formulation that explicitly incorporates initial and boundary conditions, together with a Green-identity-based treatment that distinguishes purely algebraic transpose identities from true operator adjoints obtained by integration by parts. We validate the resulting gradient and HVP expressions against reference HVPs computed by automatic differentiation (AD) using an independent JAX computational graph, which we term FDCompGraph. FDCompGraph reproduces the same discrete time stepping, source/receiver injection, and convolutional perfectly matched layer (CPML) treatment as the PDE solver, enabling apples-to-apples validation after matching wavefields under identical discretization and CPML settings. Using the first-order velocity–stress acoustic vertically transversely isotropic (VTI) system with CPML (a representative case where L is non-symmetric), we demonstrate that incorrect handling of $\mathcal {H}_1$ can yield errors up to approximately 400 per cent for HVP in ε and approximately 40 per cent for HVP in δ, whereas the HVPs in Vp and ρ are numerically unaffected. We further assess numerical errors unique to HVP evaluation, which requires perturbed wavefields and their time/space derivatives. In acoustic VTI tests using the damping sponge boundary condition, we observe HVP sensitivity losses of approximately 10 per cent amplitude compared with the results obtained using the CPML boundary condition. These errors may be difficult to diagnose because the commonly used finite-difference check, H δm ≈ g(m + δm) − g(m), is itself a linear approximation, particularly when the gradient g is computed using wavefields affected by the same imperfect boundary absorption. To make accurate HVPs practical at scale, we develop Devito–Hdm, a storage- and I/O-efficient Devito-based engine for co-simulating the required forward, adjoint, tangent, and second-adjoint system, with CPML-consistent checkpointing and partially on-the-fly correlations. As an application, we perform three-parameter VP–ε–δ acoustic VTI FWI using Newton–CG across models with increasing scatterer density and compare it against standard NLCG and L-BFGS. This application shows that incorporating HVP information overall yields a more reliable recovery of small-scale scatterer perturbations at approximately the same computational cost.