{"id":"bae122a4-10bd-467e-b2b9-037f0074680f","arxiv_id":"2506.03617","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"MatExPre is a matrix-exponential-based preconditioner that solves high-frequency Helmholtz equations with linear growth in matrix-vector products as frequency increases.","lead":"The paper introduces MatExPre, a new preconditioner that speeds up solving the high-frequency Helmholtz equation, which models waves in materials, using only sparse matrix-vector multiplications. It may make large seismic imaging and other wave simulations faster and more parallel.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The unproved λ*>0 is not the bottleneck because the shift in Algorithm 1 guarantees a positive gap; the real load-bearing gap is the missing control of GMRES/FGMRES on non-normal systems behind the O(ω) SpMV claim.","rationale":"The reader's CONDITIONAL verdict is appropriate, and the reader correctly flags the unproved spectral gap as a weakness of Theorems 2.2 and 3.1. My stress-test shifts the emphasis: for the shifted operator actually used in the main experiments, the inequality Im λ(A_s) ≥ ω²s λ_min(D) makes the gap positive unconditionally, so the original conjecture is less load-bearing than it appears. The sharper gap is that Theorem 3.1 concerns a fixed truncated polynomial P_N, whereas Algorithm 1 uses inner GMRES and outer FGMRES; the only quantitative link from spectra to the claimed O(ω) cost is inequality (12), which depends on the uncontrolled eigenvector condition number cond(X). The paper explicitly acknowledges the limits of spectral analysis for non-normal GMRES, yet the O(ω) complexity claim is supported only by measured iteration tables in a limited frequency range. The released code and extensive benchmarks are genuine positive evidence, and the method appears to work in the tested cases, so this is not a refutation. It identifies the precise condition that would need to be checked or proved before the O(ω) claim can be considered robust, and it suggests a concrete measurement that would settle the question.","tokens_in":21152,"tokens_out":18197,"duration_ms":223408,"concrete_test":"Run the released PETSc/SLEPc code with -log_view and record the exact total SpMV counts, including all SLEPc Krylov iterations inside every exp/ψ1 evaluation, for the 2D homogeneous model at freq=40, 80, 160, 320 and, if resources allow, freq=640, with s=1/freq, t=0.4/freq^2, and fprtol=0.08. If total SpMVs per solve grows faster than a constant multiple of freq, the central O(ω) claim fails. As a theory probe, compute full eigendecompositions of the 1D and small 2D PML matrices at the same frequencies and report cond(X); unbounded growth of cond(X) would invalidate the bound (12) as a justification of the observed inner GMRES scaling.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The paper's convergence theory rests on the conjectured strict spectral gap λ* = min_k Im λ_k > 0 (Section 2.1, before Thm 2.2; Section 3.3). The reader is right that this is unproved for the PML/discretized operator. However, for the algorithm actually used in the complexity claims (Alg. 1, Eq. (13)), the preconditioner is built from A_s = A + iω²sD with D a positive-definite diagonal matrix. If A has nonnegative imaginary part, the shift alone forces Im λ(A_s) ≥ ω²s λ_min(D) > 0; with s=1/ω and t=O(ω^{-2}), t λ*_s = O(ω^{-1}), which is exactly the mechanism behind the observed O(ω) inner iterations. Thus the conjectured λ*>0 for the original A is not the most load-bearing condition for the shifted solver. The load-bearing gap is the step from this spectral information to convergence and cost of the nested algorithm. Theorem 3.1 analyzes the fixed polynomial P_N, but Algorithm 1 solves (I−exp(itA_s))w = r by inner GMRES and wraps it in outer FGMRES; the only quantitative bridge is (12), which contains the eigenvector condition number cond(X). The paper itself notes GMRES on non-normal matrices cannot be inferred from eigenvalues [30], and nothing controls cond(X) or the outer FGMRES iteration count for A A_s^{-1} in the regime s=1/ω (where [26] predicts possible deterioration for small shifts). The 7–14 outer iterations and linear inner counts are empirical; if cond(X) or outer counts grow with frequency/model complexity, the O(ω) complexity claim fails. The Marmousi counts (152 at freq=340, 269.9 at freq=680) already show a large non-normality inflation, so the margin is not large.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript proposes MatExPre, a preconditioner for the discretized Helmholtz equation built from the identity A^{-1} = -it(Σ_{n=0}^∞ ψ_0(itA)^n)ψ_1(itA). The authors derive a Schrödinger-like time-domain fixed-point solver, connect it to exponential integrators, analyze the spectral clustering of the truncated polynomial P_N A, and introduce a complex shift A_s = A + iω²sD. Numerical experiments on 2D and 3D homogeneous models and on the SEG/EAGE salt, overthrust, and Marmousi-II benchmarks report 7–14 outer FGMRES iterations and inner GMRES counts that grow roughly linearly with frequency, leading to the claim that the number of SpMV operations grows linearly with frequency.","tokens_in":21575,"tokens_out":10382,"duration_ms":112545,"significance":"If the O(ω) SpMV claim holds robustly, MatExPre would be a useful contribution to high-frequency Helmholtz preconditioning because it avoids matrix factorizations and heavy sequential routines, uses only sparse matrix-vector products, and has a clean algebraic foundation. The identity (8) is parameter-free and elegant, Theorem 3.1 gives a transparent spectral clustering result under Im λ ≥ λ*>0, and the paper is honest about its conjectured strict spectral gap and about known GMRES limitations on non-normal matrices. The source code is publicly available, which strengthens reproducibility. However, the central complexity claim currently rests on numerical extrapolation rather than on a theorem controlling GMRES/FGMRES on non-normal systems, and the reported evidence for linear growth is not uniform across all tables.","major_comments":[{"comment":"The claim that the total number of SpMV operations grows linearly with frequency is not established by the presented theory. Theorem 3.1 analyzes the fixed polynomial P_N, but Algorithm 1 solves (I−exp(itA_s))w = r with inner GMRES and wraps this in outer FGMRES; the total SpMV count also includes the SLEPc evaluations of ψ1(itA_s) and exp(itA_s) inside each inner iteration. The lower-bound argument in Section 3.4 is only a lower bound, so it cannot justify an upper bound. The tables report iteration counts, not total SpMV counts, and the outer count and the per-iteration matrix-function cost are never combined into a single complexity measure. Please either provide a theorem or a careful accounting for the total SpMV count, or state the linear-growth claim explicitly as empirical.","section":"Section 3.4, Algorithm 1, and the Conclusion"},{"comment":"The theory does not control non-normality. After noting, with reference [30], that GMRES convergence cannot be inferred from eigenvalues alone, the only quantitative bridge is inequality (12), which contains cond(X). No control of cond(X) is given for the PML-discretized Helmholtz matrices, and no bound is provided for the outer FGMRES iteration count on A A_s^{-1} in the regime s = 1/ω, where [26] predicts possible deterioration for small shifts. The numerical linear-growth results in Tables 3, 5, 7, 8, and 9 could therefore be specific to the test problems. Please add a discussion or analysis of non-normality, or at least a systematic numerical study of cond(X) and the total SpMV counts.","section":"Section 3.1 and inequality (12)"},{"comment":"The convergence theorems rest on a finite spectral decomposition of the source g and on the conjectured strict inequality λ*>0; Proposition 2.1 only proves Im λ ≥ 0, and the manuscript explicitly labels λ*>0 as a conjecture. As written, this theory does not directly cover Algorithm 1, which uses the shifted operator A_s. Since Im λ(A_s) ≥ ω²s λ_min(D) > 0, the shift removes the need for the conjecture for the implemented method. Please state this explicitly and give the corresponding version of Theorem 3.1 for A_s, including the resulting scaling of t λ*_s, so that the link between the theory and the observed O(ω) inner iteration counts is explicit and not left to the reader.","section":"Section 2.1 and Theorems 2.2/2.3"},{"comment":"The numerical evidence for linear growth is mixed. In Table 2 (s = 0), the first-FGMRES inner iterations are 15, 25, 43, 108 for freq = 40, 80, 160, 320; the ratio over an 8-fold frequency increase is 7.2, and the increase from 160 to 320 is 108/43 ≈ 2.5 for a 2-fold frequency increase. In Table 8, the overthrust model shows inner iterations 8.3 at freq = 40 and 24.8 at freq = 80, a threefold increase for a twofold frequency increase, which is superlinear. Please report total SpMV counts, including matrix-function evaluations, and either fit the scaling or qualify the linear-growth conclusion accordingly.","section":"Tables 2 and 8; Section 4.1.2"}],"minor_comments":[{"comment":"The displayed identity appears to miss a complex conjugation: taking \\bar{Φ}/c² as the test function gives the stated result, whereas the written ∫ Φ∆Φ dx is not equal to −∫|∇Φ|² dx for complex Φ. The conclusion is correct after conjugation, but the proof should be corrected.","section":"Proposition 2.1 proof"},{"comment":"The word 'convergences' should be 'converges' in both theorem statements.","section":"Theorems 2.2 and 2.3"},{"comment":"The estimate ||(ω²I+H)^{-1}|| ≈ O(ω^{-1}) is described as a postulate; this is a wavenumber-explicit regularity estimate and should be cited or proven rather than presented as an ad-hoc assumption.","section":"Section 2.2"},{"comment":"The column alignment in Tables 6 and 7 is difficult to read because frequency, DoF count, outer iterations, and inner iterations are packed together without clear spacing; please reformat these tables.","section":"Tables 6 and 7"},{"comment":"The phrase 'increases almost linearly in freq' for Table 2 should be quantified, since some rows (for example, the 320-row values 108, 76, 95, 110) do not follow a clean linear trend.","section":"Section 4.1.2"},{"comment":"The left panel plots two time series with different scales but a single 'time (s)' axis label; consider separate panels or a legend that clarifies the two scales.","section":"Figure 6"},{"comment":"The abstract states 'rigorous connection', while the body repeatedly relies on conjectures and formal spectral decompositions; consider softening 'rigorous' to 'algebraic' or 'exact at the discrete level'.","section":"Abstract"}],"recommendation":"major_revision","confidential_remarks":"The paper is a reasonable candidate after revision. The main concern is not the conjectured λ*>0—the shift in Algorithm 1 makes a positive spectral gap available—but the absence of a total SpMV accounting and the lack of control of GMRES/FGMRES on non-normal systems. I would not reject on this basis, because the algebraic framework is clean and the numerical experiments are substantial, but the authors should add total SpMV counts, verify solution accuracy against reference solutions, and either prove or explicitly qualify the linear-growth complexity claim."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Colleague,\n\nThis is a useful paper to know about. It constructs a Helmholtz preconditioner from the identity A^{-1} = -it (sum_{n>=0} exp(intA)) psi_1(itA), evaluates the matrix functions by Krylov methods, and wraps the solve in FGMRES with a complex shift. The algebra is clean, and the experimental section covers serious 3D benchmarks (SEG/EAGE salt and overthrust, Marmousi-II) with the code public. Relative to Luo-Liu [37], the new pieces are the PML formulation and the complex shift; the identity itself is standard, but using it this way as a nested Krylov preconditioner is a legitimate contribution.\n\nWhat is done well: Theorem 3.1 is correct under diagonalizability and a positive spectral gap; the numerical-strategy section is honest about non-normality, citing Greenbaum-Ptak-Strakos and exposing cond(X) in (12); the experiments are extensive and the code is available. The authors also state their limitations openly, for example that spectral radius alone does not predict GMRES convergence on inhomogeneous models.\n\nThe main gap is not the conjectured lambda* > 0. The shifted operator A_s = A + i omega^2 s D gives a positive imaginary spectral gap without that conjecture, so that part is repairable and not load-bearing in the algorithm actually tested. The load-bearing gap is the jump from spectral clustering to a cost theorem for the nested solvers. There is no wavenumber-explicit control of cond(X), no outer FGMRES theory, and the O(omega) claim is supported by experiments plus a support-propagation heuristic that works only as a lower bound. The Marmousi numbers, 152 inner iterations at freq 340 and 269.9 at 680, show non-normality inflation, so the margin is not huge. A referee should ask for either a sharper theorem or a clear statement that O(omega) is empirical.\n\nThe other soft spot is the absence of comparisons with existing preconditioners (shifted Laplace with multigrid, sweeping, source transfer, WaveHoltz, Luo-Liu). Without those, it is hard to calibrate how much this advances the state of the art.\n\nI would send this to peer review. It deserves referee time, with major revision expected: head-to-head comparisons, a precise proven-versus-observed statement on complexity, and ideally a condition-number study of the fixed-point operator. If you work on Helmholtz preconditioning, cite it; I would also put it on a reading-group list.","headline":"A clean and reproducible algebraic Helmholtz preconditioner with honest numerics; the O(omega) complexity is empirical, not proven, but the paper deserves serious review.","tokens_in":22124,"tokens_out":5055,"would_cite":true,"duration_ms":53876,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65F08","65F10","65F60","65N22","35J05"],"pacs":[],"model":"deepseek-v4-flash","headline":"MatExPre approximates the inverse Helmholtz operator by a truncated matrix-exponential series and, with a 1/ω complex shift, solves high-frequency Helmholtz systems with total sparse matrix-vector work that grows only linearly in frequency.","keywords":["Helmholtz equation","matrix exponential preconditioner","high-frequency waves","Schrödinger-like time-domain solver","GMRES and FGMRES","complex shift","spectral gap","acoustic wave modeling"],"falsifier":"Compute all eigenvalues of a 2D or 3D PML Helmholtz matrix at several increasing frequencies and check whether min_k Im(λ_k) stays bounded away from zero; if it decays like 1/ω or touches zero, the radius exp(-(N+1)tλ*) stops shrinking and the inner solve stalls. Alternatively, run the reported fixed-ppw experiments at freq = 640 or 1280 and test whether total inner GMRES iterations still grow linearly in ω rather than quadratically.","tokens_in":20961,"feed_emoji":"⚡","tokens_out":7822,"duration_ms":86702,"temperature":0.7,"pith_summary":"The paper introduces MatExPre, a preconditioner for the discretized high-frequency Helmholtz equation built entirely from sparse matrix-vector products. It starts from a Schrödinger-like time-domain reformulation and turns the resulting fixed-point iteration into an algebraic matrix-exponential identity for the inverse operator. The key estimate places the eigenvalues of the preconditioned matrix inside a small disk around 1, so nested Krylov solvers converge quickly. Numerical experiments on 2D and 3D homogeneous, lens, and benchmark seismic models report outer FGMRES iterations of 7–14 and inner iterations growing linearly with frequency, implying total work O(ω) at fixed points per wavelength.","feed_headline":"Matrix-exponential Helmholtz solver hits linear-in-frequency cost","feed_subtitle":"Outer iterations hold near 7–14; total sparse matrix-vector work grows only linearly with frequency.","key_machinery":"The load-bearing object is the matrix-exponential identity $A^{{-1}}$ = -it(Σ_{n=0}^∞ ψ_0(itA)^n)ψ_1(itA), truncated to N+1 terms as P_N. The ψ_l functions are the standard exponential-integrator entire functions defined by ψ_0(z)=e^z and ψ_{l+1}(z)=(ψ_l(z)-1/l!)/z; they turn the time integral of exp(iτA) into a single matrix-function application. P_N shares eigenvectors with A, and the eigenvalues of P_N A are 1 - exp(i(N+1)tλ_k), so convergence is controlled by the imaginary spectral gap λ* = min_k Im(λ_k). The inner GMRES solve of (I-exp(itA))z = r accelerates the geometric series, the outer FGMRES handles the preconditioner's variability, and the complex shift A_s = A + iω²sD is the tuning mechanism that keeps the gap effective at high frequency.","core_discovery":"The central claim is that the inverse of a PML-damped Helmholtz matrix A can be approximated by the truncated ψ-function series P_N = -it(Σ_{n=0}^N ψ_0(itA)^n)ψ_1(itA), and that this approximation is good enough to precondition GMRES at high frequency. The identity comes from writing the time integral of a matrix exponential as tψ_1(itA) and inverting I - exp(itA) by its geometric series. For any eigenvalue λ of A with Im(λ) ≥ λ* > 0, the eigenvalues of the preconditioned matrix P_N A lie in a disk centered at 1 with radius exp(-(N+1)tλ*), which explains why the inner fixed-point equation converges. The paper conjectures λ* > 0 for the damped Helmholtz operator and supports it with a non-negativity proof plus 1D eigenvalue scans, then demonstrates that a complex shift A + iω²sD enlarges the effective gap and stabilizes the outer solver. The resulting algorithm requires only sparse matrix-vector products and achieves, in the reported tests, total work linear in ω.","pith_inferences":["The paper's O(ω) lower-bound argument via wave-support spreading is heuristic; a rigorous complexity statement would need to bound inner GMRES convergence in terms of eigenvector conditioning, which the paper leaves open. If that conditioning degrades with model complexity, the method's advantage on strongly heterogeneous media would shrink.","The radius estimate exp(-(N+1)tλ*) with t ≈ O(ω^{-2}) would naively predict poor approximation unless N grows with ω; the observed linear scaling therefore likely owes as much to GMRES's clustering behavior on non-normal matrices as to the disk radius, a mechanism the paper does not isolate.","A natural test of the spectral-gap conjecture is to use a few Krylov iterations on the PML matrix to estimate λ* on the fly and choose t adaptively; this could remove the hand-tuned t ∝ 1/ω² rule.","If λ* can be shown bounded below for PML-truncated acoustic operators, the same Schrödinger-exponential construction may extend to time-harmonic Maxwell and elasticity, where the operator is also non-self-adjoint after damping; the paper does not claim this."],"forward_implications":["With s = 1/ω and t ∝ 1/ω², outer FGMRES iterations stay nearly constant across frequency (7–14 in the reported tests), while inner GMRES iterations grow roughly linearly, so total SpMV work is O(ω).","The preconditioner needs no matrix factorization and no sequential sweeping pass, so its parallel scalability is set by SpMV and Krylov inner products; reported strong-scaling efficiency is near 1 at freq = 40 and around 80% at freq = 80, with weak scaling near 80%.","The same construction handles PML, first-order absorbing, and sound-hard boundary conditions without modification, as shown on 2D and 3D variable-velocity models including benchmark seismic data.","Eigenvalue clustering gives a direct design rule: choose t and N so (N+1)tλ* is large enough that exp(-(N+1)tλ*) is small, and set the complex shift s to balance outer versus inner iteration counts.","Relaxing the inner tolerance to about 0.08 leaves outer iteration counts unchanged while cutting inner work, which is useful at large scale."],"supporting_citations":[{"why":"The WaveHoltz paper that supplies the time-domain fixed-point principle MatExPre recasts through a Schrödinger propagator.","marker":"[3]"},{"why":"The fixed-point matrix-exponential Helmholtz solver this work extends by adding PML and a modular algebraic formulation.","marker":"[37]"},{"why":"The restarted Krylov method used to apply ψ0(itA) and ψ1(itA) with only sparse matrix-vector products.","marker":"[17]"},{"why":"The action-of-matrix-exponential algorithm whose norm-based iteration estimate justifies choosing t ∝ 1/ω².","marker":"[1]"},{"why":"The exponential-integrator theory defining the ψ-functions and the variation-of-constants identity behind the inverse formula.","marker":"[34]"},{"why":"The shifted-Laplacian convergence theory that informs the complex-shift scaling s = O(1/ω).","marker":"[26]"},{"why":"The variational Sommerfeld radiation-condition analysis used to argue that imaginary parts of eigenvalues are positive.","marker":"[41]"},{"why":"The rigorous limiting-amplitude-principle result cited to justify the time-domain-to-frequency-domain connection.","marker":"[4]"}],"fun_headline_variants":["Matrix-exponential preconditioner cracks high-frequency Helmholtz","Helmholtz at high frequency, solved in linear time","Truncated matrix exponentials precondition Helmholtz fast","Linear-in-frequency Helmholtz via matrix exponential trick","High-frequency Helmholtz gets a matrix-exponential boost"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"Everything rests on the spectral gap λ* = min_k Im(λ_k) being strictly positive for the damped or PML-truncated Helmholtz operator; the paper proves only non-negativity, conjectures λ* > 0 from the Sommerfeld radiation condition and 1D scans, and every convergence rate is exponential in -λ*.","fun_headline_variants_meta":{"raw":{"variants":["Matrix-exponential preconditioner cracks high-frequency Helmholtz","Helmholtz at high frequency, solved in linear time","Truncated matrix exponentials precondition Helmholtz fast","Linear-in-frequency Helmholtz via matrix exponential trick","High-frequency Helmholtz gets a matrix-exponential boost"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000591,"raw_usage":{"total_tokens":2755,"prompt_tokens":911,"completion_tokens":1844,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":527,"completion_tokens_details":{"reasoning_tokens":1768}},"tokens_in":527,"tokens_out":1844,"duration_ms":14848,"temperature":1.0,"reasoning_tokens":1768,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-07T10:58:22.339172+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Compute all eigenvalues of a 2D or 3D PML Helmholtz matrix at several increasing frequencies and check whether min_k Im(λ_k) stays bounded away from zero; if it decays like 1/ω or touches zero, the radius exp(-(N+1)tλ*) stops shrinking and the inner solve stalls. Alternatively, run the reported fixed-ppw experiments at freq = 640 or 1280 and test whether total inner GMRES iterations still grow linearly in ω rather than quadratically.","supporting_citations":[{"cited_title":"Appel\\\" o , F","cited_arxiv_id":null,"evidence_quote":"The WaveHoltz paper that supplies the time-domain fixed-point principle MatExPre recasts through a Schrödinger propagator."},{"cited_title":"Luo and Q","cited_arxiv_id":null,"evidence_quote":"The fixed-point matrix-exponential Helmholtz solver this work extends by adding PML and a modular algebraic formulation."}],"review_version":1}