Pith. sign in

REVIEW 4 major objections 5 minor 14 references

Introduction to optimization methods for training SciML models

T0 review · 4 major / 5 minor · reviewed 2026-08-03 · deepseek-v4-flash

Pith's one-line read Training physics-constrained neural networks is governed by the PDE operator's spectrum, not by data statistics, so stochastic first-order methods often fail and curvature-aware or deterministic optimizers are needed.

desk verdict A useful tutorial on SciML optimization, but its central 'operator spectrum governs loss landscape' claim is a heuristic presented as a derivation, and a few formal slips need fixing. read the letter →

arxiv 2601.10222 v2 pith:ITNKW42J submitted 2026-01-15 math.NA cs.AIcs.NAmath.OC

classification math.NAcs.AIcs.NAmath.OC MSC 65K1068T0790C5365N22
keywords scientificmachinelearningphysics-informedneuralnetworksoptimizationtangentkernelspectralbiaspreconditioningstochasticgradientdescentcollocationsampling
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper argues that optimization problems in scientific machine learning are structurally different from those in classical machine learning. In classical ML, losses are sums of independent sample terms, so stochastic gradient methods make sense. In physics-informed settings, the loss includes a differential operator that couples every point in the domain, making the landscape stiff, anisotropic, and dominated by the operator's spectral properties rather than by the dataset. Because of this, the paper claims, effective training of PINNs generally requires large-batch or deterministic algorithms, curvature-aware methods, and careful choice of collocation points, while data-driven SciML models can retain standard stochastic training.

What carries the argument

The neural tangent kernel (NTK) is the central object. For a network h_θ, the NTK is K(x,x') = ∇_θ h_θ(x)^T ∇_θ h_θ(x'), and its empirical version Θ = (1/m) J J^T shares nonzero eigenvalues with the Gauss-Newton Hessian. The paper uses the NTK to show that GD acts as a kernel gradient flow, that error modes decay as (1 - α λ_i)^k, and that preconditioned methods effectively replace the NTK by Θ^{(M)} = (1/m) J M^{-1} J^T. A Fourier scaling argument — applying p derivatives multiplies a Fourier mode by ω^p, so H ~ ω^{2p} — is the load-bearing mechanism that ties operator order to ill-conditioning. The paper also uses Weyl's inequality on the sum of rank-one NTK contributions to show that adap

What would settle it

Measure the Hessian or NTK eigenvalue spectrum of a practical finite-width PINN on a nonlinear PDE (e.g., Burgers' equation or Allen–Cahn) and check whether the condition number actually scales as (ω_max/ω_min)^{2p} and whether NTK eigenvalues separate cleanly by frequency. If not, the coupling is more complex than the Fourier argument claims.

Watch

Extended reading notes

Core claim

The paper's central claim is that the optimization landscape of physics-constrained SciML is shaped by the spectrum of the underlying differential operator, not by data statistics. Using the neural tangent kernel framework, it shows that gradient descent in parameter space induces a kernel flow in function space, and that the eigenvalues of the NTK govern how quickly different error modes decay. The paper derives a Fourier-mode scaling argument: for a differential operator of order p, the Jacobian scales as ω^p and the Hessian as ω^{2p}, so the condition number of the objective grows as (ω_max/ω_min)^{2p}. This explains why higher-order PDEs are harder to train, why spectral bias is stronger

Load-bearing premise

The central claim rests on treating the network output as a sum of independent Fourier modes that the PDE operator multiplies by ω^p, and on assuming the Hessian is dominated by the Gauss-Newton term; in a finite-width network with nonlinear activations and a nonlinear operator, this is an idealization.

Editorial extensions

If this is right

  • In physics-constrained PINN training, mini-batch SGD can stagnate because stochastic noise projects onto flat NTK directions while informative directions are frozen; deterministic or very-large-batch regimes are often preferable.
  • Collocation-point sampling is an optimization parameter: adding points in high-residual regions raises the smallest NTK eigenvalue, improving the condition number without changing the network.
  • The conditioning penalty worsens with differential order: for a p-th order PDE, the condition number scales as (ω_max/ω_min)^{2p}, so higher-order operators need more aggressive preconditioning.
  • Hybrid schedules that switch from Adam to L-BFGS after the gradient norm stabilizes can combine robust early progress with fast asymptotic convergence.
  • Adaptive gradient methods and second-order methods act as preconditioners that flatten the effective NTK spectrum, narrowing the gap between fast and slow error modes.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The Fourier scaling argument makes a testable prediction: for a fixed architecture, a PINN's Hessian or NTK condition number should grow roughly with the operator order p; measuring this on a nonlinear PDE could show where the decoupled-mode idealization breaks down.
  • If the operator-spectrum view is correct, then architecture choices such as Fourier features and multiscale representations are not merely accuracy aids — they are optimization preconditioners that reshape the NTK spectrum.
  • For data-driven SciML (operator learning, surrogate models, latent dynamics), the paper's own logic implies that standard stochastic methods remain appropriate because the loss decomposes over samples; the challenge is only in physics-constrained terms.
  • The same spectral lens could guide the design of spatially stratified mini-batches, explicitly keeping the effective kernel well-conditioned during stochastic training rather than relying on uniform random sampling.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 5 minor

Summary. This paper is an introductory survey/tutorial on optimization methods for training scientific machine learning (SciML) models. It contrasts classical ML, where finite-sum, sample-separable objectives are typically trained with stochastic first-order methods, with physics-constrained SciML, where PDE-induced losses are claimed to be non-separable, globally coupled, anisotropic, and stiff. The paper argues that the conditioning of these problems is governed by the spectral properties of the underlying differential operator, motivating deterministic and curvature-aware optimizers such as L-BFGS, inexact Newton, and hybrid Adam-to-L-BFGS schedules, while data-driven SciML retains the characteristics of standard stochastic ML. It reviews the neural tangent kernel, spectral bias, SGD convergence theory, adaptive sampling, momentum/adative methods, quasi-Newton and Hessian-free methods, and recent preconditioning approaches organized by data, parameter, and function space. Worked numerical examples and accompanying notebooks are provided.

Significance. If the central narrative is sound, this survey offers a useful synthesis for practitioners and researchers, making the case that problem structure should drive optimizer choice in SciML. The paper is valuable as a teaching resource: it derives NTK-based function-space dynamics, gives explicit mode-wise decay rates, provides a helpful taxonomy of preconditioning strategies, and ships Jupyter notebooks that make the numerical experiments reproducible. Many parts are honest about open problems, repeatedly noting the absence of principled guidelines and theoretical guarantees. However, the manuscript's core dichotomy between physics-constrained and data-driven SciML rests on a heuristic Fourier-scaling argument that is presented as established fact, and the paper contains several formal misstatements in its convergence theorems and definitions. The survey is likely to be useful to the community if these load-bearing issues are corrected and the central claims are appropriately qualified.

major comments (4)
  1. [Section 3.4 (Impact of physical constraints on loss landscape and the NTK of PINNs)] The derivation from the Fourier decomposition h_θ = Σ ω ĥ(ω)e^{iωx}, through ∂^p_x h_θ ~ ω^p, J_k ~ ω^p, H(θ)=(1/m)J_k^T J_k ~ ω^{2p}, to the condition-number scaling κ(H)≈(ω_max/ω_min)^{2p} is load-bearing for the paper's central claim. As written, it assumes (i) Fourier modes diagonalize the network Jacobian, (ii) the PDE operator acts as a scalar multiplier ω^p (true only for linear constant-coefficient operators), and (iii) the Hessian is dominated by the Gauss-Newton term. For finite-width nonlinear networks and nonlinear operators such as Burgers, modes couple, the NTK is not constant, and higher-order Hessian terms are not negligible. The text does not flag these as assumptions. Since this scaling motivates the practical recommendations in Sections 5 and 7 (deterministic/large-batch training, L-BFGS, operator-spectrum preconditioning), the authors should either re-frame the argume
  2. [Section 4.1.3 (Theorem 4)] Theorem 4 states a bound on min_{k=0,...,K−1} E∥∇f(θ̄_K)∥^2, but θ̄_K is not defined in this theorem; it is the averaged iterate from Theorem 2. The standard nonconvex SGD guarantee is min_k E∥∇f(θ_k)∥^2, not a statement about the averaged iterate. In addition, the step size α = sqrt(2/(L_f L_max K)) uses L_f, which is undefined, and L_max is only introduced in Assumption 1 with the index n in place of m. Because this theorem is used to conclude that large curvature (large L_max) slows SGD, the statement should be corrected and its proof or source checked.
  3. [Definition 1 and Assumption 1 (Section 4.1)] Definition 1 defines a smooth function by ∥∇f(θ1)−∇f(θ2)∥ ≤ ∥θ1−θ2∥, which arbitrarily fixes the Lipschitz constant to 1 and leaves the stated constant L unused. The standard definition is ∥∇f(θ1)−∇f(θ2)∥ ≤ L∥θ1−θ2∥. In Assumption 1, L_max is written as max_{i=1,...,n} L_i, but the summands f_i are indexed by i=1,...,m. These errors propagate into the constants in Theorems 1–4 and weaken the paper's claims about how the condition number affects convergence.
  4. [Abstract and Section 1 vs Sections 5.3, 8.1–8.3] The abstract and Section 1 assert that SciML losses are 'often non-separable and globally coupled' and that optimization behavior is 'governed by the spectral properties of the underlying physical models.' Yet the manuscript itself repeatedly acknowledges that the supporting theory is incomplete: Section 5.3.1 says 'the SciML literature currently provides no principled guidelines for constructing mini-batches tailored to PINNs'; Section 8.2 states 'principled guidelines for constructing stable parameter decompositions and convergence guarantees in stochastic regimes remain largely absent'; Section 8.3 says multilevel strategies' theoretical guidance ... remains only partially understood.' The strong introductory claims should be softened to reflect this reality, or the paper should explicitly label the operator-spectrum viewpoint as a working hypothesis rather than an established dichoto
minor comments (5)
  1. [Section 3.1.1] 'where ˙θ(t) denotes the derivatives of ˙θ with respect to t' should read 'the derivative of θ with respect to t.'
  2. [Proof of Theorem 1 (Section 4.1.1)] The line 'Ek∥θk+1−θ⋆∥2 ≤ (1−αµ)−2α(f(θk)−f(θ⋆)) + α^2_k σ^2' is missing the norm term and a parenthesis; it should read '(1−αµ)∥θk−θ⋆∥² − 2α(f(θk)−f(θ⋆)) + α^2_k ...'.
  3. [Assumption 1 (Section 4.1)] The index in 'L_max := max_{i=1,...,n} L_i' should be i=1,...,m, consistent with the number of summands in (2).
  4. [Algorithm 1 (Section 4)] The input line specifies 'an initial iterate x_0', but the iterates are denoted θ throughout. Use θ_0 for consistency.
  5. [Section 4.2 (Figure 7)] In the caption of Figure 7, 'SGD (blue) utilizes |I_k|=16 and α_k=10−3' might be confusing because the blue curve is the baseline; consider rewording to 'baseline SGD'.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: the spectral-stiffness narrative is a heuristic review claim, not a derivation from fitted outputs or author-imported uniqueness results.

full rationale

This is a tutorial/review rather than a derivation whose conclusion is built into its inputs. The central claim—that physics-constrained SciML losses are stiff, anisotropic, and governed by the differential operator spectrum—is supported in Section 3.4 by a Fourier-mode scaling argument (h_theta ~ sum of Fourier modes, ∂_x^p ~ ω^p, H ~ ω^{2p}). That argument is an idealization: it assumes decoupled Fourier modes, scalar ω^p action of the PDE operator, and dominance of the Gauss-Newton Hessian term. These are substantive heuristic assumptions and potential correctness risks, but they are not circular: no parameter is fitted to a target quantity and then renamed as a prediction, and no claimed result is defined in terms of the very quantity it is supposed to explain. The NTK material in Sections 3.1–3.3 is standard external theory (Jacot et al. [32], Lee et al. [47]). The self-citations appearing in Sections 8.2–8.3 (e.g., [43], [44], [48], [57], [72]) are survey pointers to multilevel and domain-decomposition training methods; they are not used as load-bearing justification for the central spectral-stiffness thesis, and no author-imported uniqueness theorem is invoked. The manuscript also repeatedly labels its practical recommendations as heuristic or open, which further indicates the absence of an overclaimed derivation. Accordingly, no circular step meeting the evidence threshold is present.

Assumptions & free parameters 0 free parameters · 4 assumptions · 0 invented entities

The central claim of this paper is expository, so the ledger contains no fitted parameters and no invented entities. The arguments rely on two heuristic domain assumptions — the frozen-NTK/Gauss-Newton approximation of the Hessian and a Fourier-mode scaling for PDE operators — plus standard matrix facts such as Weyl's inequality. These are appropriate for a survey, but they are idealizations for nonlinear finite-width networks.

assumptions (4)
  • domain assumption Empirical NTK/Gauss-Newton approximation: the Hessian is dominated by (1/m)J^T J and the NTK is nearly constant during training.
    Used in Sections 3.1-3.2 to identify Hessian and NTK spectra and to derive the mode-decay recursion e_{k+1} ≈ (I−αΘ)e_k. Strictly valid in the wide-network/lazy-training regime, which is an idealization for finite-width PINNs.
  • domain assumption Fourier-mode scaling: the PDE operator multiplies a Fourier component of frequency ω by ω^p, giving J ~ ω^p and H ~ ω^{2p}.
    Used in Section 3.4 to conclude κ(H) ≈ (ω_max/ω_min)^{2p}. This ignores mode coupling introduced by nonlinear network maps and general nonlinear differential operators.
  • standard math Weyl monotonicity for PSD matrices: adding a positive semidefinite term A+B with B ⪰ 0 does not decrease any eigenvalue.
    Invoked in Section 5.2 to argue that adding collocation points can only increase λ_min of the empirical NTK. This is standard and correct.
  • standard math SGD convergence results assume smoothness, (strong) convexity, and bounded gradient variance.
    Section 4.1 reviews standard theorems from the literature; the review does not prove them.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Introduction to optimization methods for training SciML models." pith.science (2026). https://pith.science/paper/ITNKW42J

@misc{pith2026260110222,
  author       = {Pith},
  title        = {Pith review of: Introduction to optimization methods for training SciML models},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/ITNKW42J}},
  note         = {Machine review of arXiv:2601.10222}
}
read the original abstract

Optimization is central to both modern machine learning (ML) and scientific machine learning (SciML), yet the structure of the underlying optimization problems differs substantially across these domains. Classical ML typically relies on stochastic, sample-separable objectives that favor first-order and adaptive gradient methods. In contrast, SciML often involves physics-informed or operator-constrained formulations in which differential operators induce global coupling, stiffness, and strong anisotropy in the loss landscape. As a result, optimization behavior in SciML is governed by the spectral properties of the underlying physical models rather than by data statistics, frequently limiting the effectiveness of standard stochastic methods and motivating deterministic or curvature-aware approaches. This document provides a unified introduction to optimization methods in ML and SciML, emphasizing how problem structure shapes algorithmic choices. We review first- and second-order optimization techniques in both deterministic and stochastic settings, discuss their adaptation to physics-constrained and data-driven SciML models, and illustrate practical strategies through tutorial examples, while highlighting open research directions at the interface of scientific computing and scientific machine learning.

Figures

Figures reproduced from arXiv: 2601.10222 by the authors.

Figure 1
Figure 1. Left: Example of binary classification in the [PITH_FULL_IMAGE:figures/full_fig_p004_1.png] view at source ↗
Figure 2
Figure 2. Example of spectral bias to learn target function [PITH_FULL_IMAGE:figures/full_fig_p014_2.png] view at source ↗
Figure 3
Figure 3. Comparison of the convergence behaviors of GD, the adaptive gradi [PITH_FULL_IMAGE:figures/full_fig_p015_3.png] view at source ↗
Figures from the paper (8 more)
Figure 4
Figure 4. Figure 4: Two-dimensional projection [49] of the objective function landscape [PITH_FULL_IMAGE:figures/full_fig_p017_4.png]
Figure 5
Figure 5. Figure 5: Trajectories of GD (blue) and SGD on a smooth quadratic function [PITH_FULL_IMAGE:figures/full_fig_p021_5.png]
Figure 6
Figure 6. Figure 6: Effect of batch size on convergence of the mini-batch GD method [PITH_FULL_IMAGE:figures/full_fig_p024_6.png]
Figure 7
Figure 7. Figure 7: Example of reducing SGD noise by using learning rate scheduling [PITH_FULL_IMAGE:figures/full_fig_p025_7.png]
Figure 8
Figure 8. Figure 8: (Left) Domain Ω = (0, 1) with localized forcing near x ≈ 0.9. Uni￾form sampling is depicted in orange, while adaptive refinement is depicted in blue. (Right) Level sets of θ ⊤Hmθ = 1 become more isotropic after refinement, reflecting improved conditioning. 5.3 Mini-bat…
Figure 9
Figure 9. Figure 9: Sketch of heavy-ball momentum (left) and NAG (right) iteration. [PITH_FULL_IMAGE:figures/full_fig_p031_9.png]
Figure 10
Figure 10. Figure 10: Example of performance of SGD and adaptive first-order methods [PITH_FULL_IMAGE:figures/full_fig_p034_10.png]
Figure 11
Figure 11. Figure 11: PINN training for Poisson problem in 1D using Adam, L-BFGS, and [PITH_FULL_IMAGE:figures/full_fig_p038_11.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

14 extracted references · 5 linked inside Pith

  1. [1]

    Ahamed, N

    [1]S. Ahamed, N. Zakariaei, E. Haber, and M. Eliasof,Multiscale training of convolutional neural networks, arXive-prints, (2025), pp.arXiv–

  2. [2]

    [42]D. A. Knoll and D. E. Keyes,Jacobian-free Newton–Krylov methods: a survey of approaches and applications, Journal of Computational Physics, 193 (2004), pp. 357–397. [43]A. Kopaničáková, H. Kothari, G. E. Karniadakis, and R. Krause, Enhancing training of physics-informed neural networks using domain decomposition–based preconditioning strategies, SIAM ...

  3. [9]

    [52]L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis,Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nature Machine Intelligence, 3 (2021), pp. 218–229. [53]L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis,DeepXDE: A deep learning library for solving differential equations, SIAM review, 63 (2021), p...

  4. [34]

    ,Extended physics-informed neural networks (XPINNs): A general- ized space-time domain decomposition based deep learning framework for nonlinear partial differential equations, Communications in Computational Physics, 28 (2020). [35]A. D. Jagtap, E. Kharazmi, and G. E. Karniadakis,Conserva- tive physics-informed neural networks on discrete domains for con...

  5. [62]

    [66]M. A. Nabian, R. J. Gladstone, and H. Meidani,Efficient training of physics-informed neural networks via importance sampling, Computer- Aided Civil and Infrastructure Engineering, 36 (2021), pp. 962–977. [67]J. Nocedal and S. Wright,Numerical optimization, Springer Science & Business Media,

  6. [1998]

    Verfürth,A posteriori error estimation and adaptive mesh-refinement techniques, Journal of Computational and Applied Mathematics, 50 (1994), pp

    [88]R. Verfürth,A posteriori error estimation and adaptive mesh-refinement techniques, Journal of Computational and Applied Mathematics, 50 (1994), pp. 67–83. [89]N. Vyas, D. Morwani, R. Zhao, M. Kwun, I. Shapira, D. Brand- fonbrener, L. Janson, and S. Kakade,Soap: Improving and stabiliz- ing Shampoo using Adam, arXiv preprint arXiv:2409.11321, (2024). [9...

  7. [2006]

    Ponce, R

    [68]C. Ponce, R. Li, C. Mao, and P. Vassilevski,Multilevel-in-width training for deep neural network regression, Numerical Linear Algebra with Applications, 30 (2023), p. e2501. [69]A. Quarteroni, P. Gervasio, and F. Regazzoni,Combining physics- based and data-driven models: advancing the frontiers of research with sci- entific machine learning, arXiv pre...

  8. [2012]

    [30]Z. Hu, A. D. Jagtap, G. E. Karniadakis, and K. Kawaguchi,Aug- mented Physics-Informed Neural Networks (APINNs): A gating network- based soft domain decomposition methodology, Engineering Applications of Artificial Intelligence, 126 (2023), p. 107183. [31]B. Irwin and E. Haber,Secant penalized BFGS: a noise robust quasi- Newton method via penalizing th...

Show all 14 references
  1. [2014]

    [78]C. J. Shallue, J. Lee, J. Antognini, et al.,Measuring the effects of data parallelism on neural network training, JMLR, (2019). [79]H.-J. M. Shi, Y. Xie, R. Byrd, and J. Nocedal,A noise-tolerant quasi-Newton algorithm for unconstrained optimization, SIAM Journal on Optimiz...

  2. [2017]

    Kiyani, K

    [40]E. Kiyani, K. Shukla, J. F. Urbán, J. Darbon, and G. E. Kar- niadakis,Optimizing the optimizer for physics-informed neural networks and Kolmogorov-Arnold networks, Computer Methods in Applied Mechan- ics and Engineering, 446 (2025), p. 118308. [41]A. Klawonn, M. Lanser, an...

  3. [2018]

    [82]H. Son, J. W. Jang, W. J. Han, and H. J. Hwang,Sobolev training for physics informed neural networks, arXiv preprint arXiv:2101.08932, (2021). [83]S. Subramanian, R. M. Kirby, M. W. Mahoney, and A. Gholami, Adaptive self-supervision algorithms for physics-informed neural n...

  4. [2019]

    46 [48]Y. Lee, A. Kopaničáková, and G. E. Karniadakis,Two-level overlap- ping additive Schwarz preconditioner for training scientific machine learn- ing applications, Computer Methods in Applied Mechanics and Engineer- ing, 448 (2026), p. 118400. [49]H. Li, Z. Xu, G. Taylor, C...

  5. [2021]

    [20]Z. Gao, L. Yan, and T. Zhou,Failure-informed adaptive sampling for pinns, SIAMJournalonScientificComputing, 45(2023), pp.A1971–A1994. [21]G. Garrigos and R. M. Gower,Handbook of convergence theorems for (stochastic) gradient methods, ArXiv preprint 2301.11235, (2023). [22]...

  6. [2501]

    Amari,Natural gradient works efficiently in learning, Neural Com- putation, 10 (1998), pp

    [2]S.-i. Amari,Natural gradient works efficiently in learning, Neural Com- putation, 10 (1998), pp. 251–276. [3]E. Amid, R. Anil, and M. Warmuth,Locoprop: Enhancing backprop via local loss optimization, in International conference on artificial intelligence and statistics, PML...

Pith tools

Reviewed August 3, 2026 · model on record in the stance chip above.