{"id":"aadc0ea8-176a-47db-9d6f-bf58c5381c22","arxiv_id":"2507.10804","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":6,"one_line_summary":"A combined PSF and pseudo-differential-operator probing method accelerates seismic full-waveform inversion and enables more trustworthy MCMC uncertainty quantification for high-rank Hessian problems.","lead":"This paper introduces a faster way to approximate the Hessian operator of seismic full-waveform inversion, and shows it speeds up both deterministic imaging and Bayesian uncertainty quantification. The authors combine two existing approximation methods and validate on a synthetic quadratic model and the Marmousi benchmark, reporting more reliable posterior sampling than low-rank or prior-only baselines.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Marmousi UQ validation uses a low-rank 'benchmark' that the paper elsewhere argues underestimates variance; without a rank-converged or independent reference, the central claim that high-rank Hessians yield accurate posterior variances remains unverified.","rationale":"The reader's weakest_assumption is the spatial and frequency locality of the misfit Hessian symbol, which is indeed a methodological assumption and is acknowledged in Section 5. I considered that concern, but the paper provides empirical support for locality in Figs 4, 5, 12, and 13, and the idealized quadratic model is intentionally constructed to satisfy the low-rank-symbol assumption, so locality alone is not the most decisive gap. The more decisive gap is the absence of any true posterior reference for the realistic Marmousi test. The low-rank 'benchmark' is not a ground truth: low-rank truncation is precisely the failure mode the paper documents for pCN and gpCN-LR, and the Table 4 note explicitly says the compared chains have not converged. Without either a rank-convergence study or exact marginal variances, the central claim rests on comparing one approximation to another approximation of the same type. The ideal quadratic model mitigates this, but it is a synthetic operator built to satisfy the method's assumptions; it cannot carry the realistic claim alone. The proposed test, exact CG solves for the five marginal variances, would settle whether PSF+ sample variances are accurate. This concern does not require changing the reader's CONDITIONAL verdict; it sharpens and justifies it.","tokens_in":19325,"tokens_out":6218,"duration_ms":81749,"concrete_test":"Compute exact posterior marginal variances at the five evaluation points of Fig. 16 by solving H_map y_i = e_i with preconditioned CG using the exact matrix-free Hessian-vector product, and compare these values with the PSF+ MCMC sample variances shown in Fig. 18. Also compute the same comparison for a rank-converged low-rank reference, e.g., r = 500 or larger until the generalized eigenvalues of (H_d, R) decay, to test whether the 'brute-force benchmark' is itself converged in rank. If the PSF+ variances deviate from the exact CG values beyond the Monte Carlo error of the chains, the central claim that high-rank Hessian approximations accurately reflect the posterior variance is not supported for Marmousi.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim (Section 6) is that MCMC-gpCN samples using high-rank Hessian approximations 'more accurately reflect the variance of the posterior distribution.' The ideal quadratic model supports this because the true posterior is known exactly, but the Marmousi validation in Section 4.2.2 has no exact reference. There, the authors establish a benchmark as a 'brute-force low-rank approximation of the Hessian at the MAP model, followed by MCMC-gpCN sampling.' This is the load-bearing weakness. Section 3.1 and the Section 4.1.3 results show that low-rank approximations with insufficient rank underestimate posterior variance and can produce misleading ESS values, and Section 3.7 states that low-frequency eigenvalues are large and significantly influence samples. If the benchmark uses too small a rank, or if its MCMC chains are not converged, then its marginal variances and histograms are themselves biased. The Table 4 note explicitly concedes that the compared chains have not converged, and no convergence diagnostics are reported for the benchmark chain. 'Closely matches the benchmark' would then only indicate that PSF+ shares the benchmark's low-rank truncation error, not that it reproduces the posterior. Section 5 also acknowledges that the methods rely on spatial and frequency locality, which may fail with sharp parameter variations, and Appendix A3 concedes that PDO overestimates high-frequency components; these are real approximation limits, but the benchmark issue is more central because it undermines the evidence for the strongest claim in the realistic setting.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper adapts two established high-rank Hessian approximation methods (the point spread function method and the pseudo-differential operator probing method) and proposes a new PSF+ method that combines row and column symbol information to approximate the misfit Hessian of a seismic inverse problem. These approximations are used as preconditioners for L-BFGS in deterministic full-waveform inversion and as the proposal covariance in a generalized preconditioned Crank-Nicolson (gpCN) MCMC sampler for Bayesian uncertainty quantification. The methods are validated on an ideal quadratic model with analytically known posterior covariance and on the Marmousi model, where a low-rank Hessian-based gpCN run is used as the reference. The paper concludes that high-rank Hessian approximations are essential for accurate posterior variance estimation and can substantially improve MCMC efficiency.","tokens_in":19634,"tokens_out":5193,"duration_ms":63851,"significance":"If the central claim holds, the paper would make a meaningful contribution to large-scale Bayesian seismic UQ by showing that high-rank Hessian approximations can capture the directional scalings of the posterior in a regime where low-rank approximations are known to fail. The ideal quadratic model is a strong, falsifiable benchmark because the true posterior covariance is known analytically and the PSF+ method reproduces it, and the PSF/PDO duality derivation in Eqs. (39)-(42) is algebraically consistent. The paper also demonstrates concrete L-BFGS preconditioning benefits. However, the Marmousi validation currently lacks an independent reference, and the paper itself reports non-converged chains and acknowledges approximation limitations (Section 5, Appendix A3), so the strongest claims about posterior variance accuracy are not yet fully supported.","major_comments":[{"comment":"The Marmousi benchmark is constructed as a 'brute-force low-rank approximation of the Hessian at the MAP model' followed by MCMC-gpCN sampling, but Sections 3.1 and 3.7 argue that low-rank approximations with insufficient rank underestimate posterior variance, and Table 4's own note states that the compared chains 'have not yet converged.' Therefore the statement that gpCN-PSF+ 'closely matches the benchmark' could simply reflect a shared low-rank truncation or non-convergence bias. The authors should provide a higher-rank or independently verified reference (for example, a Hessian approximation with explicit rank-convergence checks, or an independent sampling method with convergence diagnostics) before the Marmousi experiment can support the central claim about posterior variance accuracy.","section":"Section 4.2.2, Table 4"},{"comment":"The ESS values for gpCN-PSF+ on the Marmousi model (7.05, 6.14, 6.57, 4.88, 16.00) are not systematically higher than those of pCN or gpCN-LR, yet the text dismisses ESS and autocorrelation as 'heuristic' while relying on histograms that are compared only with the low-rank benchmark. Since no quantitative distance between the PSF+ marginal histograms and the benchmark is reported, and since non-convergence can affect all chains including the benchmark, the evidence for improved posterior variance accuracy in Marmousi is currently visual and indirect. Quantitative comparisons (for example, variance ratios against a trusted covariance, coverage of a reference interval, or a suitable distance between marginals) with convergence diagnostics are needed.","section":"Section 4.2.2, Table 4"},{"comment":"The construction of UQ proposal samples via pointwise square roots of symbol rows and columns requires discarding negative symbol values, but the effect of this truncation on the resulting covariance is not analyzed. In addition, Appendix A3 states that the PDO method 'tends to overestimate high-frequency components' and that this 'poses challenges for UQ' because the proposal random vector is not band-limited; nevertheless, PDO is used as a UQ method in Sections 4.1.3 and 4.2.2. The paper should quantify the resulting bias in posterior marginal variances, for example by comparing PDO-based posterior variance against the analytic reference in the ideal quadratic model, and explain why the acknowledged high-frequency overestimate does not invalidate the UQ conclusions.","section":"Section 3.7 and Appendix A3"},{"comment":"The methods rely on the assumption that the misfit Hessian has strong locality in both spatial and frequency domains, and Section 5 acknowledges that this 'may be violated in models with sharp parameter variations.' The Marmousi model contains strong velocity contrasts including the 400 m water layer, but no quantitative approximation error (for example, the relative operator norm error ||Hd - H_approx||_F / ||Hd||_F) is reported for the PSF, PDO, or PSF+ approximations on this model. A numerical assessment of the approximation error and its spatial distribution would be needed to rule out that the locality assumption is silently violated in the very setting used for the UQ validation.","section":"Section 5 and Figures 12-13"}],"minor_comments":[{"comment":"There are several typographical errors: 'self-adjont' should be 'self-adjoint', 'Too see' should be 'To see', and the Figure 1 caption 'An illustration of the duality ... this shown' should be 'this is shown.' These should be corrected.","section":"Section 3.5, Eq. (39) and Figure 1"},{"comment":"The sentence 'we introduce low-rank approximation as both a baseline and an auxiliary technique...' contains a duplicated word ('are are') in the preceding paragraph: 'two established Hessian approximation methods that are are widely used.' Please fix.","section":"Section 3, introductory paragraph"},{"comment":"The MCMC settings are not fully specified: the step-size parameter beta in Eq. (14), the burn-in length, the number of chains, and the total number of samples after burn-in are not reported. Providing these details in the appendix would make the experiments reproducible and help interpret the reported ESS values.","section":"Sections 4.1.3 and 4.2.2"},{"comment":"The definition of the symbol class is imprecise. The standard definition of S^1 requires derivative bounds of the form |∂_ξ^α ∂_x^β s(x,ξ)| ≤ C_{αβ}(1+|ξ|)^{1-|α|}; the text's wording 'bounded by polynomials of ξ of corresponding orders' is too vague. Similarly, Eq. (35) is an asymptotic expansion, not an exact polynomial expansion, and should be described as such.","section":"Section 3.3, Eq. (27)"},{"comment":"The wave simulation code is proprietary, and the paper states it can be substituted with any standard code that supports gradient and Hessian-vector products. It would be helpful to specify at least the grid spacing, time-stepping scheme, source frequencies, and number of wave equation solves per Hessian-vector product, so that the reported numerical results can be meaningfully reproduced with an open-source solver.","section":"Data Availability"}],"recommendation":"major_revision","confidential_remarks":"The strongest part of the paper is the ideal quadratic model experiment, where the analytic reference makes the PSF+ variance improvement unambiguous. The Marmousi UQ section is the weakest support for the central claim because the benchmark is itself a low-rank approximation and the paper concedes the chains have not converged. I would suggest requesting a rank-convergence study or an independent reference for the Marmousi case, and a quantitative treatment of the acknowledged PDO high-frequency overestimate and the negative-symbol truncation in Section 3.7. The proprietary wave simulation code is a reproducibility limitation that the editor may wish to flag."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Here's my read. The genuinely new thing is PSF+, which merges the PSF and PDO methods by exploiting the duality between symbol rows and columns, and it's well motivated. The ideal quadratic model gives a clean ground truth and the results there are convincing: PSF+ tracks the true posterior standard deviation much better than pCN or gpCN-LR, and ESS is clearly higher. That's real evidence. The L-BFGS preconditioning results are also solid and show the Hessian approximations help deep recovery.\n\nThe soft spots are concentrated in the Marmousi UQ section. The benchmark there is a brute-force low-rank approximation of the Hessian, with no rank-convergence check. Since the paper itself shows that low-rank approximations with insufficient rank underestimate posterior variance, \"closely matches the benchmark\" doesn't tell you much. The authors even note in Table 4 that the chains have not converged, and they present no diagnostics for the benchmark chain. So the strongest claim—that high-rank approximations more accurately reflect posterior variance in a realistic setting—remains unverified. This is a genuine weakness, but it's addressable: a rank-converged reference or an independent ground-truth estimate would fix it. The PDO high-frequency overestimation (Appendix A3) is real but secondary, and the locality assumption is openly acknowledged in Section 5. The proprietary wave solver is an annoyance for reproducibility, though the approximation code is on GitHub.\n\nOverall, the paper is worth a serious referee. The algorithmic contribution is novel and the ideal quadratic model is a legitimate benchmark. The Marmousi validation needs work before publication, but the core idea stands.","headline":"A genuinely new PSF+ Hessian approximation with a convincing ideal-model test, but the Marmousi UQ validation leans on a low-rank benchmark that undercuts the paper's strongest claim.","tokens_in":20219,"tokens_out":1564,"would_cite":true,"duration_ms":17794,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65N21","86A15","65C05","62F15"],"pacs":[],"model":"deepseek-v4-flash","headline":"High-rank Hessian approximations built from point-spread functions and pseudo-differential probing make seismic inversion and Bayesian uncertainty quantification practical, with MCMC variances that track the true posterior.","keywords":["seismic full-waveform inversion","Bayesian uncertainty quantification","high-rank Hessian approximation","point spread function method","pseudo-differential operator probing","MCMC-gpCN","Marmousi model","Hessian preconditioning"],"falsifier":"Run PSF+ on a synthetic model with a sharp high-contrast interface, compare the approximate Hessian's action on many random vectors against exact Hessian-vector products, and check whether gpCN samples still reproduce the variance of a converged reference sampler; growing relative error with contrast strength, or sample variance falling below the reference, would show the locality assumption has failed.","tokens_in":19113,"feed_emoji":"🌊","tokens_out":12697,"duration_ms":134974,"temperature":0.7,"pith_summary":"The paper sets out to make seismic full-waveform inversion and Bayesian uncertainty quantification tractable by approximating the misfit Hessian—the operator that encodes local curvature of the inversion objective as well as the covariance of the Laplace-approximated posterior—without forming it explicitly. Its central claim is that high-rank approximations built from point-spread functions and from pseudo-differential-operator probing, plus a new combined method called PSF+, reproduce the directional scalings of the posterior well enough that a generalized preconditioned Crank-Nicolson (gpCN) Markov chain Monte Carlo sampler explores the high-dimensional model space efficiently. On a synthetic quadratic problem with known uncertainty and on the Marmousi model, the authors find that samples from unapproximated or low-rank proposals concentrate in narrow regions, underestimate posterior variance, and overstate effective sample size, while the high-rank approximations yield variances that track the true posterior. If correct, this gives a practical route to trustworthy uncertainty estimates in seismic imaging at a fraction of the cost of exact Hessian construction.","feed_headline":"Seismic uncertainty sampling gets a high-rank Hessian speedup","feed_subtitle":"MCMC samples that used to underestimate variance now track the true posterior at far lower wave-equation cost.","key_machinery":"The central object is the Hessian $H = H_d + \\Gamma_{pr}^{-1}$ of the negative log-posterior, accessed only through matrix-vector products that each cost two wave-equation solves per source. The argument treats $H_d$ as a pseudo-differential operator with symbol $s(x,\\xi)$, a phase-space scaling function; PSF and PDO are dual ways of sampling that symbol, respectively along spatial rows via delta-function responses and along frequency columns via sinusoidal probing vectors. The load-bearing identities are $H_d\\phi_\\xi(x)=\\phi_\\xi(x)s(x,\\xi)$ and $s(x_k,\\xi)=(2\\pi)^2\\widehat{p_k}(\\xi)$, which let a few operator applications recover enough of the symbol to reconstruct a low-rank separated approximation $s(x,\\xi)\\approx\\sum_{k=1}^r a_k(x)b_k(\\xi)$ that is applied in $O(rN\\log N)$ time. For UQ, the symbol is square-rooted pointwise to keep the operator symmetric positive definite, and a high-pass filter plus low-rank correction separates the well-localized mid-frequency symbol from the smooth low-frequency part.","core_discovery":"On the paper's own terms, the central claim is that the seismic misfit Hessian $H_d$, although too large to store and too high-rank for low-rank compression, can be approximated as a pseudo-differential operator with a symbol $s(x,\\xi)$ that is low-rank in a separated sense. The PSF method samples rows of this symbol by applying $H_d$ to delta functions; the PDO method samples columns by applying it to sums of sinusoids; the new PSF+ method combines both sets of samples by extracting a row basis with a singular value decomposition (SVD) and refining column coefficients through a small Tikhonov-regularized least-squares problem. After a high-pass filter removes nonlocal low-frequency content and a low-rank correction restores that content, the approximation is symmetric positive definite, so its inverse square root can seed Gaussian proposals for the gpCN sampler. The numerical evidence is that these proposals mix far faster than pCN or low-rank gpCN and, more importantly, that the spread of the samples matches the true posterior variance, whereas the baselines underestimate it and give misleadingly high effective sample sizes.","pith_inferences":["An implication the paper leaves implicit is that the row-column symbol completion strategy should transfer to other inverse problems governed by partial differential equations whose Hessians are pseudo-differential operators, such as elastic or electromagnetic inversion, though no such experiment appears here.","The paper's diagnostic lesson suggests a practical safeguard it does not itself implement: before trusting effective sample size for a Hessian-informed sampler, compare sample marginals against the Laplace approximation's marginals.","A natural testable extension, if sharp parameter contrasts violate symbol locality, is adaptive symbol sampling that places point-spread functions or probing frequencies where the symbol changes fastest, or a multiscale split that treats the nonlocal part separately."],"forward_implications":["Full-waveform inversion can use the inverse Hessian approximation as an L-BFGS initial Hessian or Newton-Krylov preconditioner, cutting iteration counts and improving recovery in deeper, less-illuminated regions.","Bayesian UQ with gpCN becomes practical for high-rank seismic posteriors, since the proposal covariance captures directional scalings that low-rank and uninformative proposals miss.","Common MCMC diagnostics such as autocorrelation and effective sample size can report good mixing even when a chain is confined to a narrow region; the paper's histograms and variance maps expose that false confidence.","The PSF+ method inherits both the spatial locality of PSF sampling and the frequency locality of PDO probing, and in the benchmark experiments it is more accurate than either method alone.","The pseudo-differential row-column structure remains valid in three dimensions, so the approach is extensible to 3D seismic problems."],"supporting_citations":[{"why":"Establishes that the seismic Hessian is a pseudo-differential operator and gives the symbol expansion used for interpolation.","marker":"Bao & Symes 1996"},{"why":"Develops the pseudo-differential scaling and probing approach that the PDO method adapts.","marker":"Nammour & Symes 2011"},{"why":"Provides the product-convolution point-spread-function approximation underlying the PSF method.","marker":"Alger et al. 2019"},{"why":"Demonstrates band-limited high-frequency symbol decay, motivating the PDO interpolation caveat and the square-root handling for UQ.","marker":"Demanet et al. 2012"},{"why":"Introduces Hessian-informed MCMC sampling for seismic posterior inference, the lineage gpCN extends.","marker":"Martin et al. 2012"},{"why":"Supplies the generalized preconditioned Crank-Nicolson proposal used for the MCMC sampler.","marker":"Pinski et al. 2015"},{"why":"Provides the comparative assessment of Hessian-informed MCMC methods and the low-rank Hessian baseline used here.","marker":"Kim et al. 2023"},{"why":"Uses pointwise symbol operations in PSF deconvolution, supporting the square-root symbol construction for symmetric positive definite approximations.","marker":"Yang et al. 2022"},{"why":"Supplies the randomized SVD low-rank Hessian machinery used as baseline and as low-frequency correction.","marker":"Villa et al. 2021"}],"fun_headline_variants":["High-rank Hessian approximations boost seismic inversion and UQ","Faster seismic inversion and uncertainty via high-rank Hessian","New Hessian approximation speeds up seismic inversion and MCMC","High-rank Hessian shortcuts accelerate seismic inversion and UQ","Seismic inversion and uncertainty get high-rank Hessian boost"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the misfit Hessian behaves like a pseudo-differential operator with a symbol smooth and separable in space and frequency, so sampled point-spread functions stay localized and probed symbol columns can be interpolated; Section 5 states this can be violated in models with sharp parameter variations, and if it fails the approximations and the uncertainty gains break down.","fun_headline_variants_meta":{"raw":{"variants":["High-rank Hessian approximations boost seismic inversion and UQ","Faster seismic inversion and uncertainty via high-rank Hessian","New Hessian approximation speeds up seismic inversion and MCMC","High-rank Hessian shortcuts accelerate seismic inversion and UQ","Seismic inversion and uncertainty get high-rank Hessian boost"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000438,"raw_usage":{"total_tokens":2294,"prompt_tokens":1081,"completion_tokens":1213,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":697,"completion_tokens_details":{"reasoning_tokens":1132}},"tokens_in":697,"tokens_out":1213,"duration_ms":9617,"temperature":1.0,"reasoning_tokens":1132,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T17:26:16.697666+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run PSF+ on a synthetic model with a sharp high-contrast interface, compare the approximate Hessian's action on many random vectors against exact Hessian-vector products, and check whether gpCN samples still reproduce the variance of a converged reference sampler; growing relative error with contrast strength, or sample variance falling below the reference, would show the locality assumption has failed.","supporting_citations":[{"cited_title":"& Symes, W","cited_arxiv_id":null,"evidence_quote":"Establishes that the seismic Hessian is a pseudo-differential operator and gives the symbol expansion used for interpolation."},{"cited_title":"Matrix probing: a randomized preconditioner for the wave-equation H essian, Applied and Computational Harmonic Analysis\\/ , 32 (2), 155--168","cited_arxiv_id":null,"evidence_quote":"Demonstrates band-limited high-frequency symbol decay, motivating the PDO interpolation caveat and the square-root handling for UQ."},{"cited_title":"C., Burstedde, C., & Ghattas, O., 2012","cited_arxiv_id":null,"evidence_quote":"Introduces Hessian-informed MCMC sampling for seismic posterior inference, the lineage gpCN extends."},{"cited_title":"J., Simpson, G., Stuart, A","cited_arxiv_id":null,"evidence_quote":"Supplies the generalized preconditioned Crank-Nicolson proposal used for the MCMC sampler."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the comparative assessment of Hessian-informed MCMC methods and the low-rank Hessian baseline used here."},{"cited_title":"An efficient and stable high-resolution seismic imaging method: point-spread function deconvolution, Journal of Geophysical Research: Solid Earth\\/ , 127 (7), e2021JB023281","cited_arxiv_id":null,"evidence_quote":"Uses pointwise symbol operations in PSF deconvolution, supporting the square-root symbol construction for symmetric positive definite approximations."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the randomized SVD low-rank Hessian machinery used as baseline and as low-frequency correction."}],"review_version":1}