{"id":"11d9165f-6aec-412a-b350-22ca534844b1","arxiv_id":"2504.16553","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"A hybrid optimization that solves the network's linear output layer by least squares inside each gradient step accelerates and stabilizes PINN training for the scattered Helmholtz equation.","lead":"This paper embeds a least-squares solver into the gradient-descent training of physics-informed neural networks, so the output layer of the network is solved exactly at each step instead of slowly learned. The result is faster and more stable convergence for Helmholtz wavefield simulations, a practical gain for seismic imaging and acoustic modeling.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Eq. 25's PML constraint block [DC, DC] does not implement Eq. 13; the Marmousi LS-GD experiment may minimize a different objective.","rationale":"Reader and I converge on the same weakest point: the PML soft-constraint block in Eq. 25 is algebraically inconsistent with Eq. 13. I considered the more general concern that gradients through the LS solve are not derived; this is less load-bearing because automatic differentiation through a differentiable solve is standard practice, the code is public, and the non-PML experiments would still support the convergence claim. The Eq. 25 block is an internal algebraic inconsistency, not a disagreement with consensus, and it directly affects the only PML experiment. I therefore keep the CONDITIONAL verdict: the core no-PML idea is coherent and reproducible in principle, but the PML results should not be cited as evidence for the stated loss until the block is corrected and the experiment rerun. I agree with the reader's weakest_assumption, so no verdict change is proposed.","tokens_in":11119,"tokens_out":8498,"duration_ms":84218,"concrete_test":"Run the released TensorFlow implementation's PML branch with the third block of D changed from [DC, DC] to [[DC,0],[0,DC]] (and R padded with a second zero block), keeping all hyperparameters and seeds; compare the Marmousi 30 Hz training-loss/MSE curves with Fig. 5. If the curves change materially or the near-source field grows, Eq. 25 is responsible and the PML claims need revision; if they are unchanged, the issue is a benign presentation error.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The PML half of the central claim rests on the least-squares system in Eq. 25. With W stacked as [W_r; W_i], the third block row contributes ||DC(W_r+W_i)||^2 = (beta/N_C) sum_j gamma_j^2 (u_r,j+u_i,j)^2 to the LS objective, whereas the stated soft constraint Eq. 13 is (beta/N_C) sum_j gamma_j^2 (u_r,j^2+u_i,j^2). The correct matrix is block-diagonal [DC,0;0,DC] with two zero blocks in R; the printed [DC,DC] introduces a spurious 2(beta/N_C) sum gamma^2 u_r u_i cross-term, so the PML loss being minimized is not the loss defined in Eqs. 10/13. If the released TensorFlow code follows Eq. 25, the Marmousi convergence curves and MSE in Fig. 5 are evidence for a different objective, and the claim that the method is correctly formulated with PML is not established. The non-PML derivation and experiments are unaffected, so the flaw is scoped to the PML portion but load-bearing for that portion.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a hybrid least-squares/gradient-descent (LS-GD) training scheme for physics-informed neural networks (PINNs) solving the frequency-domain scattered Helmholtz equation. The key idea is to replace the standard gradient-descent update of the linear output layer with an analytic damped least-squares solve, and to define the GD loss as the residual of that solve (Eq. 22), so that the hidden layers are trained against the optimal output layer at each iteration. Derivations are given for the case without PML (Eqs. 17-20) and with PML (Eqs. 23-25), together with implementation details and complexity estimates. Numerical experiments on a 10 Hz simple model and a 30 Hz Marmousi model compare the proposed LS-GD training with standard GD training of the same network architecture.","tokens_in":11344,"tokens_out":7746,"duration_ms":77153,"significance":"The non-PML least-squares derivation is mathematically sound, and the single-loss embedding is a clean variant of the alternating LS/GD schemes of Cyr et al. and Uriarte et al.; if the PML formulation is corrected, the paper would provide a practically useful and low-overhead acceleration for PINN-based Helmholtz solvers. The reported experiments show large qualitative improvements in convergence on benchmark models, and the promise of released code is a strength. However, the PML block in Eq. 25 does not implement the stated soft constraint, and because the Marmousi experiment uses Eq. 25, the current evidence for the PML half of the central claim is compromised. The contribution is incremental rather than a new optimization paradigm, but it is a meaningful one for wavefield simulation applications.","major_comments":[{"comment":"The third block row [D_C, D_C] of the PML least-squares matrix does not implement the soft constraint in Eq. 13. With the output weights stacked as W = [W_r; W_i], this block contributes ||D_C(W_r + W_i)||^2 = (β/N_C) Σ_j γ_j^2 (u_{r,j} + u_{i,j})^2 to the LS objective, whereas Eq. 13 requires (β/N_C) Σ_j γ_j^2 (u_{r,j}^2 + u_{i,j}^2). The correct representation is block-diagonal [D_C, 0; 0, D_C]. As printed, the LS problem minimizes a different objective that includes a spurious cross-term 2(β/N_C)Σ_j γ_j^2 u_{r,j} u_{i,j}. Since the Marmousi experiment (Section 4.2) uses Eqs. 22 and 25, its reported loss and MSE curves are evidence for the stated loss only if the released code actually uses the corrected block and the text is merely misprinted. Please fix the matrix, rerun the PML experiments, and state explicitly which form the code implements.","section":"Section 2.4.2, Eq. (25)"},{"comment":"The abstract and conclusion claim that the LS solver adds minimal computational overhead and that the method converges where standard GD fails, but no wall-clock times, per-epoch costs, or multiple-seed statistics are reported. The complexity analysis in Section 5 counts only the normal-matrix construction and Cholesky solve; it does not account for the cost of building D by forward-mode differentiation at every epoch or for backpropagating through the LS solve when training the hidden layers. Please report runtimes or measured per-epoch costs for the reported experiments and, for the claim that standard GD 'fails' on Marmousi, repeat the comparison with at least a few random seeds so the reader can assess robustness.","section":"Section 4, Figs. 3-5"}],"minor_comments":[{"comment":"Equation 16 introduces W without specifying its dimensions or relation to the later W*; please clarify that W here denotes the output-layer weight matrix, while W* in Eq. 21 is the damped least-squares solution used in the loss of Eq. 22.","section":"Eq. (16) and Eq. (21)"},{"comment":"The wavelength λ is not defined in the text; if it denotes the acoustic wavelength, please state the corresponding frequency and reference velocity used to evaluate it.","section":"Section 2.3, Eq. (14)"},{"comment":"Please state explicitly whether the gradient of the loss in Eq. 22 flows through the damped LS solve W* = (D^T D + εI)^{-1} D^T R or whether W* is treated as a detached constant; this is important for reproducing the reported training dynamics.","section":"Section 3"},{"comment":"The validation MSE is reported excluding the PML region; this exclusion is mentioned only in the figure caption and should be stated in the main text, together with a justification for omitting the PML zone from the reported error.","section":"Figure 5(f) and Section 4.2"},{"comment":"Please correct the typos: 'itterative' (Section 2.4.1), 'trianing' and 'gradualy' (Section 3), 'expriments' (Section 5), and 'oberved' (Section 4.1).","section":"Throughout"}],"recommendation":"major_revision","confidential_remarks":"The PML matrix error is the main obstacle; it is localized and fixable, but the Marmousi experiments must be revalidated with the corrected Eq. 25 or with a clear statement that the code already uses the correct block. I would also ask the authors to verify the released code against the formulas in the paper before acceptance. The contribution is a moderate incremental advance over the cited LS/GD hybrid methods, so the bar for acceptance should be clarity and correctness of the PML formulation plus reproducible experimental comparisons."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Bottom line: the non-PML part of this paper is a clean, useful application of Cyr et al.'s least-squares output-layer idea to the scattered Helmholtz equation. The PML part contains a real linear-algebra bug in the printed formulation, so the Marmousi results as written are evidence for a different objective than the one stated.\n\nWhat is actually new: the authors do not alternate LS and GD; they embed the LS solution into a modified GD loss (Eq. 22) and derive the matrix system for the non-PML Helmholtz case. That derivation, Eqs. 17-22, is correct. The non-PML experiments show a strong, qualitative improvement in convergence, especially with few collocation points, and pointing to a TensorFlow implementation is good practice. The use of forward-mode differentiation is also sensible.\n\nSoft spots, in proportion:\n\n1. Eq. 25 is wrong. The soft-constraint term in Eq. 13 penalizes (1/N_C) sum gamma^2 (u_r^2 + u_i^2). With W stacked as [W_r; W_i], the block row [DC, DC] contributes ||DC(W_r+W_i)||^2, i.e. (u_r+u_i)^2, which includes a spurious 2 u_r u_i cross-term. The correct matrix is block-diagonal [DC,0;0,DC]. This is not typographical: it changes the loss being minimized, so Fig. 5's convergence curves and MSE cannot be attributed to the loss in Eqs. 10/13 unless the code silently uses the corrected matrix. A referee should ask the authors to run the PML experiment with the corrected block and report whether the curves change.\n\n2. The experiments report 'best predictions' without seeds, error bars, or wall-clock times. That is a fixable discipline problem, but it makes the convergence claims weaker than they look.\n\n3. The gradient of the new loss with respect to hidden-layer parameters is asserted, not derived. Since W* depends on the penultimate features through D^T D and D^T R, backprop through the LS solve is not automatic in the usual sense; you need an implicit-differentiation formula or a careful explanation of what the autodiff graph computes. The paper skips this.\n\nThe citation pattern looks honest. The self-citations are to earlier work used as inputs (background wavefield, PML form), not as support for the new claim.\n\nVerdict: this deserves a serious referee, not a desk rejection. The non-PML half is publishable as a modest contribution on its own; the PML claim needs the matrix corrected and the experiments re-run. If the corrected formulation and code agree, I would take it. If not, a shorter paper on the non-PML case would still be worth having.","headline":"Non-PML LS-embedding for PINN Helmholtz is a clean incremental contribution; the PML half has a load-bearing matrix bug and the Marmousi results need re-running before they can be trusted.","tokens_in":11929,"tokens_out":5339,"would_cite":false,"duration_ms":52355,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"A least-squares solve for the output layer, embedded into the PINN loss, makes scattered-wavefield Helmholtz training converge fast and stably where plain gradient descent stalls.","keywords":["physics-informed neural networks","Helmholtz equation","least-squares optimization","hybrid optimization","seismic wavefield simulation","perfectly matched layer","gradient descent","scattered wavefield"],"falsifier":"Compute, at a fixed epoch of the published code, the Eq. 22 loss using the Eq. 25 matrix and compare it with the sum of the PML PDE residual (Eq. 12) and $\\beta$ times the soft constraint (Eq. 13) on the same wavefield; if the soft-constraint contribution equals $\\|D_C(W_r+W_i)\\|^2$ rather than $\\|D_C W_r\\|^2+\\|D_C W_i\\|^2$, the PML experiments minimize a different objective and the paper's PML claim is not established.","tokens_in":2041,"feed_emoji":"🌊","tokens_out":2611,"duration_ms":93741,"temperature":0.7,"pith_summary":"The paper claims that the slow, unstable training of physics-informed neural networks (PINNs) for the frequency-domain Helmholtz equation can be cured by solving the output layer's weights with a least-squares (LS) solve inside the gradient-descent loss at every step. Instead of backpropagating through all layers, the network first computes the best possible linear output layer for the current hidden features, then trains the hidden layers against that optimal layer. On scattered acoustic wavefield benchmarks, the authors report that this LS-embedded loss converges from the first epochs, reaches lower errors than standard GD, and succeeds where plain GD fails, including a 30 Hz Marmousi test with perfectly matched layers. The added cost is small because the LS solve works on a small normal matrix, where the matrix size is the width of the penultimate layer. If correct, the method gives a drop-in training upgrade for PINN-based wavefield simulation.","feed_headline":"One least-squares solve per step fixes slow wavefield neural nets","feed_subtitle":"Embedding an LS solver for the output layer beats plain gradient descent on the Helmholtz equation, even at 30 Hz.","key_machinery":"The load-bearing object is the least-squares-embedded loss $L = \\sum_j |DW^* - R|^2$ of Eq. 22. Here $D$ is a matrix assembled from the PDE residual operator acting on the penultimate-layer outputs $H$ and from the source soft-constraint term; $R$ contains the background-wavefield forcing; and $W^*$ is the damped least-squares output weight from Eq. 21. The mechanism is that each training step solves the small $P \\times P$ normal system $(D^T D + \\epsilon I)W^* = D^T R$ by Cholesky decomposition, then propagates gradients of this loss through $W^*$ back to the hidden layers using automatic differentiation, with forward-mode differentiation providing the derivatives of $H$ efficiently. This turns the linear output layer into an exactly solved subproblem and leaves gradient descent with only the nonlinear hidden features to learn.","core_discovery":"The central discovery is that the loss function for PINN training can be redefined as $L = \\sum_j |DW^* - R|^2$, where $W^* = (D^T D + \\epsilon I)^{-1} D^T R$ is the damped least-squares solution for the output weights given the current penultimate-layer features $H$. Because $W^*$ is the exact minimizer of the quadratic PDE-plus-constraint loss for the output layer, every gradient-descent step optimizes the hidden layers against the best possible linear readout rather than a readout partway through a gradient trajectory. The paper derives $D$ and $R$ for the scattered Helmholtz equation in two settings: with real and imaginary parts decoupled when no perfectly matched layer is used, and with a block matrix that couples them when PML coordinate stretching is included. The numerical claim is that this change alone, with no architectural modification, makes 10 Hz and 30 Hz scattered wavefield PINNs converge rapidly and stably, even with 500 collocation points or only 8 penultimate neurons, where standard GD training stalls.","pith_inferences":["Editorial inference: The same loss reformulation should transfer to any PDE residual that is linear in a final linear network layer, such as time-harmonic elasticity, diffusion, or Maxwell-type equations, where the LS step would be a drop-in change to the training loop.","Editorial inference: Because the LS solve removes the output-layer coordinates from the GD landscape, the hidden layers become the sole nonlinear feature learners; one testable expectation is that architecture choices such as depth, width, and activation will matter more under LS-GD than under plain GD.","Editorial inference: The method's per-batch exact solves also suggest a natural combination with adaptive sampling: if collocation batches are chosen to emphasize high-residual regions, the LS output layer will reflect that emphasis immediately, potentially making the method more responsive than GD to adaptive collocation strategies."],"forward_implications":["With the LS-embedded loss, the PINN converges from the first epochs even with very few collocation points, so the number of collocation points required per epoch can be drastically reduced compared to plain GD training.","The method consistently achieves lower validation error than standard GD on the simple velocity model across different collocation counts and different penultimate layer widths, including widths as small as 8 neurons.","In the Marmousi 30 Hz test with PML, the LS-GD method reaches in about 15,000 epochs an accuracy that plain GD does not reach after 150,000 epochs, and it does so with markedly smaller loss fluctuations.","The added computational cost of the LS step is modest, scaling as $O((N+N_C)P^2)$ to assemble the normal matrix and $O(P^3)$ for the Cholesky solve, where $P$ is the penultimate layer width.","The same loss reformulation applies both with and without perfectly matched layers, so the training upgrade can be used in realistic absorbing-boundary seismic simulations."],"supporting_citations":[{"why":"Establishes the hybrid least-squares/gradient-descent strategy for linear-output-layer networks that this paper adapts to the Helmholtz PINN loss.","marker":"(Cyr et al., 2020)"},{"why":"Supplies the scattered Helmholtz equation and the real/imaginary separated PDE loss used for the no-PML case.","marker":"(Alkhalifah et al., 2021)"},{"why":"Provides the PML-augmented scattered Helmholtz equation and the analytical background wavefield used in the PML formulation.","marker":"(Abedi et al., 2025)"},{"why":"Introduces the perfectly matched layer technique on which the PML attenuation factors are built.","marker":"(Berenger, 1994)"},{"why":"The source of the soft-constraint term $L_C$ that keeps the scattered solution from blowing up near the source.","marker":"(Huang and Alkhalifah, 2023)"},{"why":"Defines the physics-informed neural network training framework that the paper modifies.","marker":"(Raissi et al., 2019)"},{"why":"Background for forward-mode differentiation and least-squares optimization in variational PINNs, cited for the derivative-computation strategy.","marker":"(Uriarte et al., 2025)"},{"why":"Introduces the Marmousi velocity model used as the challenging benchmark.","marker":"(Versteeg, 1994)"},{"why":"Source of the sinusoidal positional encoding applied to the input coordinates.","marker":"(Huang et al., 2021)"}],"fun_headline_variants":["One least-squares solve per step accelerates wavefield PINNs","Embedded LS solver fixes slow acoustic wavefield PINNs","Exact output-layer solve boosts Helmholtz PINN convergence","LS-in-loss speeds PINNs on scattered wavefields"],"cache_read_input_tokens":13952,"weakest_assumption_plain":"The method assumes the PML soft-constraint term can be written as one shared block applied to the sum of the real and imaginary wavefield parts, but the written penalty requires separate blocks on each part; if the implementation follows the matrix in Eq. 25, the Marmousi training optimizes a different objective than the paper states, and the reported PML convergence would not demonstrate the claimed loss.","fun_headline_variants_meta":{"raw":{"variants":["One least-squares solve per step accelerates wavefield PINNs","Embedded LS solver fixes slow acoustic wavefield PINNs","Exact output-layer solve boosts Helmholtz PINN convergence","LS-in-loss speeds PINNs on scattered wavefields"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.00026,"raw_usage":{"total_tokens":1600,"prompt_tokens":964,"completion_tokens":636,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":580,"completion_tokens_details":{"reasoning_tokens":569}},"tokens_in":580,"tokens_out":636,"duration_ms":6815,"temperature":1.0,"reasoning_tokens":569,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-16T11:00:49.525645+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Compute, at a fixed epoch of the published code, the Eq. 22 loss using the Eq. 25 matrix and compare it with the sum of the PML PDE residual (Eq. 12) and $\\beta$ times the soft constraint (Eq. 13) on the same wavefield; if the soft-constraint contribution equals $\\|D_C(W_r+W_i)\\|^2$ rather than $\\|D_C W_r\\|^2+\\|D_C W_i\\|^2$, the PML experiments minimize a different objective and the paper's PML claim is not established.","supporting_citations":[{"cited_title":"C., Lan, Y., and Xu, J","cited_arxiv_id":null,"evidence_quote":"Establishes the hybrid least-squares/gradient-descent strategy for linear-output-layer networks that this paper adapts to the Helmholtz PINN loss."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the scattered Helmholtz equation and the real/imaginary separated PDE loss used for the no-PML case."},{"cited_title":"Gabor-Enhanced Physics-Informed Neural Networks for Fast Simulations of Acoustic Wavefields","cited_arxiv_id":"2502.17134","evidence_quote":"Provides the PML-augmented scattered Helmholtz equation and the analytical background wavefield used in the PML formulation."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Introduces the perfectly matched layer technique on which the PML attenuation factors are built."},{"cited_title":"and Alkhalifah, T","cited_arxiv_id":null,"evidence_quote":"The source of the soft-constraint term $L_C$ that keeps the scattered solution from blowing up near the source."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Defines the physics-informed neural network training framework that the paper modifies."},{"cited_title":"M., and Rojas, S","cited_arxiv_id":null,"evidence_quote":"Background for forward-mode differentiation and least-squares optimization in variational PINNs, cited for the derivative-computation strategy."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Introduces the Marmousi velocity model used as the challenging benchmark."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Source of the sinusoidal positional encoding applied to the input coordinates."}],"review_version":1}