{"id":"0d5a0787-1798-491b-b649-63d1961d2c97","arxiv_id":"2608.05862","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"NestyNet, a segmented softplus model with analytic derivatives and a second-order optimizer, reports 2,100x better function values and 1,400x better derivatives than Adam-trained MLPs on AI Feynman.","lead":"Physics-style functions are hard for standard neural networks to fit accurately, especially their derivatives. This paper introduces NestyNet, a modular surrogate model whose analytic derivatives and tailored optimizer reportedly beat standard networks by hundreds to thousands of times on a 120-equation physics benchmark.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Headline improvement factors hinge on a single fixed MLP baseline; a stronger standard baseline could shrink the reported gains.","rationale":"The reader's weakest_assumption identifies the single most load-bearing concern: the headline improvement factors are measured against one fixed MLP baseline, and the paper provides no evidence that this baseline is representative of standard neural networks in scientific ML. The central claim of the paper is precisely these factors, so if the baseline is weak or undertrained, the quantitative claim is overstated. This is not an internal inconsistency in the method—the analytic derivatives, the LM optimizer, and the released code are real assets—but a comparator problem that a targeted experiment can settle. I agree with the reader that this is load-bearing. I also note a minor internal discrepancy: the abstract's improvement factors differ from the medians in Figure 2 by roughly 10%, and the abstract does not state which aggregation produced its numbers. That is a secondary issue but reinforces the need to verify the headline quantities. The paper has genuine strengths: the same-model optimizer comparison in §8.1 is a good internal control, the code is released with checkpoints and provenance, and the limitations section is candid about scope. The concern is addressable: a hyperparameter sweep of the MLP baseline on a representative subset would show whether the gap persists. If the gap remains large after a fair sweep, the central claim holds; if not, the factors need revision. Therefore the correct verdict is unchanged from the reader's CONDITIONAL: the paper should be accepted pending this verification and the release of the deferred SPLA theory, but the baseline representativeness must be checked before the headline numbers are taken at face value.","tokens_in":31348,"tokens_out":10978,"duration_ms":105793,"concrete_test":"Take a representative subset of 20 AI Feynman equations spanning the reported improvement-factor distribution. For each, train the same MLP baseline (two hidden layers, 64 tanh units, double precision) with a modest but standard hyperparameter sweep: widths {64, 256, 1024}, depths {2, 4}, activations {tanh, SiLU}, with and without residual connections, AdamW with cosine annealing over 5 learning rates in [1e-4, 1e-2], plus full-batch L-BFGS from both random and Adam-initialized starts. Use the same 1.5e4/5e3 data splits and evaluation protocol as §8.3. For each equation, take the best baseline validation error per metric (value, gradient, second-derivative) and recompute the median improvement factors of NestyNet over this best-of-sweep baseline.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central empirical claim—median improvement factors of 2,100x/1,400x/780x over Adam and 540x/450x/250x over L-BFGS-refined MLPs—rests on the comparison in §8.3 against a single fixed baseline: a two-hidden-layer tanh MLP with 64 units per layer, trained in double precision. The paper reports no hyperparameter search, no learning-rate schedule, no training budget, and no architecture variants for this baseline, nor does it report parameter counts or training compute for either method. The internal same-model comparison (§8.1) shows the optimizer matters on the segmented model, but it does not establish that this MLP is representative of standard scientific-ML practice; a wider, residual, SiLU-activated, or Fourier-feature MLP, or one trained with a proper schedule and from multiple restarts, could plausibly achieve errors closer to NestyNet and shrink the headline factors proportionally. The comparison also conflates architecture and optimizer: the fixed MLP is never paired with a second-order solver, so the reported gap partly measures optimizer choice rather than the surrogate family. The paper's own numbers are also slightly inconsistent: the abstract quotes 2,100x/1,400x/780x and 540x/450x/250x, while Figure 2 reports median factors of 2,016x/1,572x/872x and 593x/475x/271x; the source of the abstract values is not explained. This does not change the qualitative story, but it underscores that the headline quantities need careful verification.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper introduces NestyNet, a deterministic segmented softplus neural-network surrogate whose analytic derivatives (gradients, Hessians, Laplacians) and antiderivatives are available in closed form, coupled with a predictive Levenberg–Marquardt optimizer, a segment-native preconditioning framework (SPLA), and a compositional adaptor ecosystem. The central empirical claims are benchmarked on the 120-equation AI Feynman suite: median improvements over a fixed two-hidden-layer tanh MLP trained with Adam of about 2,100x on function values, 1,400x on first derivatives, and 780x on second derivatives, with smaller but still large factors after L-BFGS refinement; an end-to-end speed advantage over vectorized autograd reaching about 44x; two PINNacle PDE wins plus an expressivity study; and an application to recovering the vertical acceleration of the Galactic disk from the collisionless Boltzmann equation.","tokens_in":31632,"tokens_out":8318,"duration_ms":80053,"significance":"If the headline comparisons survive a more rigorous baseline protocol, NestyNet would be a practically valuable substrate for scientific machine learning: the analytic derivative and antiderivative machinery is correct and cheap, the canonical initialization gives reproducible fits, and the benchmark protocol is mostly disciplined (shared splits and transforms, a single declared suite, no per-equation refits, and row-level provenance with the code release). The paper is also commendably explicit about what it does not claim: SPLA is not needed for the main results, the supervised PDE fits are not PINN solutions, and the interpolation law is an empirical regularity rather than a theorem. The main weakness is that the quantitative magnitude of the claimed gains rests on a single unsearched MLP baseline and on single-run medians, so the improvement factors should be treated as upper-bound estimates until the baseline is strengthened or the claims are narrowed.","major_comments":[{"comment":"The headline improvement factors are ratios against a single fixed baseline: a two-hidden-layer tanh MLP with 64 units per layer, trained in double precision. No hyperparameter search, learning-rate schedule, architecture variants, restarts, or training budgets are reported for this baseline, and its parameter count is not given. Because the abstract's 2,100x/1,400x/780x and 540x/450x/250x claims are measured against this one network, the representativeness of the baseline is load-bearing. Please either add a small, declared baseline-selection protocol on the same 120 equations (for example, a random search over width, depth, activation, and learning rate with the same splits and transforms) and report the resulting distribution of improvement factors, or explicitly restrict all headline claims to 'the fixed 64x64 tanh MLP baseline used here.' Without this, the quantitative magnitude of the claimed gains is not established.","section":"§8.3, Figure 2"},{"comment":"The numbers in the abstract do not match the numbers in the figures. The abstract quotes median improvements of 2,100x/1,400x/780x for Adam and 540x/450x/250x for L-BFGS, while Figure 2 reports 2,016x/1,572x/872x and 593x/475x/271x, and §8.4 reports yet another set (2,016x/960x/1,466x) for its RMS-norm analysis. Similarly, the integral gain is 5,260x in Figure 2 but roughly 4,800x in the conclusions. The source of the abstract values is not explained. Since these medians are the paper's central quantitative claims, please harmonize the values and state exactly which per-axis metrics, panels, and norm definitions produce each quoted factor.","section":"Abstract; §8.3–§8.4, Figures 2–3"},{"comment":"All improvement factors are single-run medians with no repeated-seed variability. The MLP baseline is randomly initialized, and no number of seeds, standard deviations, or per-equation variability is reported; the deterministic canonical initialization of NestyNet does not by itself determine the baseline's noise. Please report multiple seeds for the baseline (and ideally for NestyNet too) or provide a clear justification for why the medians are insensitive to seed choice, so the reader can judge whether the 98–100% win rates and the factor magnitudes are stable.","section":"§8.3, Figures 2–3"},{"comment":"The conclusion states that 'all of these results were obtained with a single, fixed hyperparameter configuration applied across the entire benchmark,' but the numerical sections use different model sizes and recipes: §8.2 crosses two model families and four training recipes, §8.3 uses dual-layer 20-segment models, and §8.7 sweeps segment counts up to 8,960 parameters. Please specify precisely which hyperparameters were fixed across which benchmark suite and which ones varied by design; as written, the reproducibility claim is ambiguous.","section":"§10; §8.2–§8.3"}],"minor_comments":[{"comment":"The caption says color encodes the benchmark equation index, but no color bar is shown; please add a color bar or remove the color encoding.","section":"Figures 2–3"},{"comment":"The column headings appear garbled ('W all (s)', 'V al.'); please typeset them properly and define each column in the caption.","section":"Table 2"},{"comment":"The statement that torch.compile 'silently miscompiles the forward-mode autograd Jacobian to zero' is surprising and affects the choice of eager autograd as the baseline; please cite a PyTorch issue or include a minimal reproducer in the released code.","section":"§8.5, footnote 8"},{"comment":"Section 8.2 retains 88 problems after excluding targets with dynamic range above four decades, while §8.3 reports results for all 120 problems with an asinh conditioning policy; please add one sentence explaining how the two protocols relate.","section":"§8.2 vs §8.3"},{"comment":"The Armijo backtracking along the LM direction is mentioned, but no default constants or safeguarding details are given; please provide the Armijo parameters or a citation.","section":"§5.2"},{"comment":"The curvature-budget inequality is stated without derivation; a one-line derivation from Eq. (8) and g(1−g)≤1/4 would help readers verify the bound.","section":"§8.4, Eq. (41)"}],"recommendation":"major_revision","confidential_remarks":"The paper is suitable in scope for astro-ph.IM and for the broader scientific-ML readership, and the code release with row-level provenance is a genuine strength. The main editorial risk is that the headline improvement factors are quoted in the abstract before the baseline-representativeness issue is resolved; I would encourage the authors to add a small baseline search or to rephrase the claims as relative to the specific fixed MLP. The companion papers are all in preparation and should not be load-bearing for the present claims."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: this paper is genuinely useful and the central engineering is sound, but treat the headline improvement factors as provisional. The integration of a segmented softplus surrogate with a carefully built Levenberg-Marquardt solver is a real contribution. The analytic derivative machinery - polylog ladder, closed-form line and box antiderivatives, O(S) Laplacian - is correct and clearly presented. The benchmark is more honest than most: same data splits, declared suite, no per-equation refits, and the GP comparison on problem #29 is a nice sanity check. The interpolation law in Section 8.4 is a known inequality, but the empirical tightness across 360 fits is interesting, and the curvature-tracking diagnostics are well argued. The soft spots are real but addressable. The headline medians come from a single fixed MLP baseline (two-hidden-layer tanh, 64 units), with no hyperparameter search, no training budget, and no error bars; a stronger baseline could shrink the factors. The abstract quotes 2100x/1400x/780x while Figure 2 reports 2016x/1572x/872x - the source of the abstract values isn't explained. That should be fixed before publication. SPLA theory is deferred to an unreleased companion paper, but the paper explicitly says SPLA wasn't needed for the main results, so the deferral is more a completeness issue than a load-bearing one. And the title's broad claim that physics functions are hard to fit with neural networks overstates what one benchmark shows. None of this changes the qualitative conclusion. The same-model optimizer comparison shows LM beating L-BFGS and Adam on the segmented model, and the analytic derivative cost advantages are concrete and measured. The code and checkpoints accompanying the release are a real asset. I'd send this to peer review, asking for a more rigorous baseline comparison and clarification of the numbers. For anyone working on scientific surrogates or PINNs, this is worth a careful look.","headline":"Solid, well-engineered framework with believable analytic-derivative gains; headline factors need a harder baseline and number reconciliation before I'd quote them.","tokens_in":590,"tokens_out":664,"would_cite":true,"duration_ms":27724,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"On the 120-equation AI Feynman benchmark, NestyNet reports median 2,100x better function values, 1,400x better first derivatives, and 780x better second derivatives than Adam-trained MLPs, with all derivatives analytic.","keywords":["neural networks","surrogate models","Levenberg-Marquardt","analytic derivatives","softplus segments","Laplacian","AI Feynman benchmark","Galactic dynamics"],"falsifier":"Train a substantially stronger neural baseline on the same 120 AI Feynman equations and data splits, for example a residual network with several hidden layers of 256 units each, well-tuned initialization, and full-batch L-BFGS; if its median function and derivative errors come within a small factor of NestyNet's, the central improvement factors shrink proportionally and the paper's premise that geometry rather than capacity is the bottleneck would be called into question.","tokens_in":31158,"feed_emoji":"🧮","tokens_out":12483,"duration_ms":97797,"temperature":0.7,"pith_summary":"Many smooth functions that matter in physics are exactly the ones standard neural networks fit badly, and this paper's claim is why: the loss landscape induced by stiff, multiscale targets is full of strongly correlated, sloppy parameter directions that first-order optimizers like Adam cannot traverse efficiently. The remedy is a coupled model-and-optimizer, NestyNet, built from a segmented softplus surrogate whose gradients, Hessians, Laplacians, and antiderivatives are closed-form, trained by a tailored Levenberg-Marquardt scheme. On the 120-equation AI Feynman benchmark the framework reports median improvements of 2,100x for function values, 1,400x for first derivatives, and 780x for second derivatives over Adam-trained MLPs, with 540x, 450x, and 250x remaining after L-BFGS refinement, and up to roughly 44x faster than vectorized automatic differentiation. The paper also derives an interpolation law: for fitted errors that are spectrally narrow-band, gradient error is the geometric mean of value and curvature errors, which turns value accuracy plus curvature control into a certificate for derivative accuracy. A sympathetic reader would care because trustworthy derivatives are what make surrogates usable for physics-informed learning, equation discovery, and inverse problems.","feed_headline":"NestyNet fits physics functions 2,100x more accurately","feed_subtitle":"A segmented softplus network trained with Levenberg-Marquardt returns exact gradients, Hessians, and Laplacians.","key_machinery":"The central object is the segmented softplus network, $f_o(x)=\\sum_s a_{os}\\sigma(z_{os})$ with $\\sigma(z)=\\log(1+e^z)$ and each preactivation $z_{os}=K_{os}\\cdot x+b_{os}$ affine in the inputs; the softplus is a smooth hinge, so the model is a sum of smooth pieces that can represent multiscale structure while remaining twice differentiable. All differential operators descend from the gate factor $g=\\sigma'(z)=1/(1+e^{-z})$: gradients are $\\sum a\\,g\\,K$, Hessians are $\\sum a\\,g(1-g)\\,KK^T$, and the Laplacian collapses to $\\sum a\\,g(1-g)\\,\\|K\\|^2$, so it costs $O(SWN_x)$ instead of forming a dense Hessian. The same softplus atom sits on the polylogarithm ladder, since $\\sigma(z)=-\\mathrm{Li}_1(-e^z)$ and $d/dz\\,\\mathrm{Li}_\\nu(-e^z)=\\mathrm{Li}_{\\nu-1}(-e^z)$, which supplies exact antiderivatives along lines and over boxes. The optimizer is a predictive Levenberg-Marquardt scheme that orders parameters by output-segment, exploits the exact split of $J^T J$ into within-segment blocks and cross-segment coupling, starts from a deterministic canonical initialization built on the identity $\\sigma(z)-\\sigma(-z)=z$, and adds quadratic path extrapolation, Woodbury fast solves, and an L-BFGS rescue for stagnation.","core_discovery":"The central claim is that the difficulty in fitting many physics functions with neural networks is not representational power but parameter-space geometry: the right parameters exist, but stiff, strongly correlated loss landscapes make them hard to find with first-order training. NestyNet meets that geometry with a deterministic segmented softplus architecture and a predictive second-order optimizer, and the paper reports that on the AI Feynman suite this combination lowers median function-value error by 2,100x, first-derivative error by 1,400x, and second-derivative error by 780x relative to a fixed two-hidden-layer tanh MLP trained by Adam, with gains of 540x, 450x, and 250x after L-BFGS refinement. Because every derivative is analytic, the improved accuracy does not come at the cost of speed: end-to-end Levenberg-Marquardt steps run up to about 44x faster than a vectorized automatic-differentiation baseline. The same machinery makes gradients, Hessians, Laplacians, and line or box antiderivatives available in closed form, and the paper demonstrates the payoff on PDE residuals and on recovering the vertical acceleration of a model Galactic disk from phase-space densities.","pith_inferences":["Beyond the benchmark, the flat cost of the analytic Laplacian in input dimension suggests high-dimensional PDE workloads such as Fokker-Planck, Schrödinger, or Hamilton-Jacobi equations could see the largest relative speedups, since autograd Laplacians cost one backward pass per dimension.","A testable diagnostic follows from the interpolation law: any fitted surrogate can report its spectral tightness ratio on held-out probes, and values well below 1 would signal broadband error where gradient accuracy cannot be inferred from value accuracy alone.","The exact polylogarithm antiderivatives could be used to integrate fitted densities without quadrature in downstream normalization or Monte Carlo tasks, a possibility the paper does not develop.","If the gains are representative of physics functions more broadly, the same coupled surrogate-optimizer recipe could improve equation-discovery pipelines that depend on stable mixed partials, rather than only pointwise fits."],"forward_implications":["Derivative-informed tasks inherit the accuracy gains: physics-informed residuals, Hamiltonian dynamics, and inverse-force estimation get closed-form Hessians, Laplacians, and mixed partials, and operator probes show the gain holds for directional derivatives, curvature contractions, and Laplacians.","Value accuracy cannot be converted into gradient accuracy without simultaneous curvature control; the interpolation law makes that coupling explicit and gives a capacity-selection rule that needs no derivative data.","Surrogates can be integrated exactly: line, axis-parallel, and hyperrectangle integrals of the segmented model are polylogarithm evaluations, so integrating a fitted surrogate adds no numerical error of its own.","The framework extends beyond the segmented family: an automatic-differentiation adaptor feeds general neural models through the same predictive LM optimizer, and the matrix-free SPLA preconditioner keeps LM feasible at nearly $10^4$ parameters with near-linear scaling."],"supporting_citations":[{"why":"Supplies the 120-equation AI Feynman benchmark on which all headline accuracy and derivative improvements are measured.","marker":"S.-M. Udrescu & M. Tegmark 2020"},{"why":"Standardizes the AI Feynman problems as SRBench, fixing domains, splits, and evaluation protocols used in the comparisons.","marker":"W. La Cava et al. 2021"},{"why":"Introduces the Levenberg damping that the predictive second-order optimizer is built around.","marker":"K. Levenberg 1944"},{"why":"Adds the curvature-adapted Marquardt damping and the normal-equations framework the paper tailors to segmented models.","marker":"D. W. Marquardt 1963"},{"why":"Defines the Adam optimizer used as the first-order baseline in all AI Feynman comparisons.","marker":"D. P. Kingma & J. Ba 2015"},{"why":"Introduces the L-BFGS quasi-Newton method used both as the refinement baseline and as the rescue step inside the optimizer.","marker":"J. Nocedal 1980"},{"why":"Provides the limited-memory BFGS implementation details that make L-BFGS refinement practical at full batch.","marker":"D. C. Liu & J. Nocedal 1989"},{"why":"Defines the physics-informed neural network residual formulation whose derivative sensitivity motivates the paper's accuracy targets.","marker":"M. Raissi et al. 2019"},{"why":"Provides the PINNacle benchmark with 22 standardized PDE problems, two of which the paper uses to demonstrate PDE transfers.","marker":"Z. Hao et al. 2024"},{"why":"Supplies the sloppy-models analysis of parameter-sensitivity eigenvalue spectra that the paper uses to explain why first-order training stalls on stiff physics targets.","marker":"M. K. Transtrum et al. 2011"}],"fun_headline_variants":["NestyNet: 2,100x more accurate physics fits with analytic derivatives","NestyNet: physics fits 2,100x better, with exact derivatives","NestyNet: analytic derivatives, 2,100x more accurate physics fits","NestyNet: physics functions 2,100x more accurate, exact derivatives","NestyNet: 2,100x accuracy gain for physics functions"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The headline accuracy gains are measured against one fixed baseline, a two-hidden-layer tanh network with 64 units per layer trained in double precision, so the claim assumes that this baseline fairly represents how standard neural networks perform on physics functions.","fun_headline_variants_meta":{"raw":{"variants":["NestyNet: 2,100x more accurate physics fits with analytic derivatives","NestyNet: physics fits 2,100x better, with exact derivatives","NestyNet: analytic derivatives, 2,100x more accurate physics fits","NestyNet: physics functions 2,100x more accurate, exact derivatives","NestyNet: 2,100x accuracy gain for physics functions"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001417,"raw_usage":{"total_tokens":5809,"prompt_tokens":1123,"completion_tokens":4686,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":739,"completion_tokens_details":{"reasoning_tokens":4592}},"tokens_in":739,"tokens_out":4686,"duration_ms":26935,"temperature":1.0,"reasoning_tokens":4592,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-07T22:08:41.381590+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Train a substantially stronger neural baseline on the same 120 AI Feynman equations and data splits, for example a residual network with several hidden layers of 256 units each, well-tuned initialization, and full-batch L-BFGS; if its median function and derivative errors come within a small factor of NestyNet's, the central improvement factors shrink proportionally and the paper's premise that geometry rather than capacity is the bottleneck would be called into question.","supporting_citations":[],"review_version":1}