{"id":"7cd306ec-f532-41f4-b060-a786bbbd24a6","arxiv_id":"2505.12646","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":1,"one_line_summary":"An implicit Hessian-vector product algorithm for FEM-based differentiable physics, built from JAX JVP/VJP primitives and validated with finite differences and Taylor tests, accelerates Newton-CG on nonlinear inverse problems.","lead":"This paper derives and implements a way to compute second derivatives (Hessians) for finite-element simulations in automatic differentiation frameworks like JAX, and tests it on four inverse-design problems. A generalist might read it to learn when exact second-order optimization beats cheaper first-order methods in physics-based design.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Finite-difference comparison in benchmarks is uncontrolled: no FD step size or tuning is reported, and SciPy's default step may be far from optimal, so the claim that exact Hessians are needed for robust Newton-CG is not yet established.","rationale":"The paper's algorithmic core—derivation of implicit Hessian-vector products and the Taylor remainder test in Fig. 5—is credible and independently checked. The Taylor test shows third-order convergence of the second-order expansion residual, which is strong evidence that the Hessian-vector product is implemented correctly. The finite-difference comparison in Section 3.3.1 also behaves as expected, with relative differences decreasing as h shrinks. The load-bearing weakness lies in the benchmark interpretation: the paper uses the failure of Newton-CG (FD) in two of four examples to argue that exact AD Hessians are necessary, but it never reports the finite-difference step size or other tuning parameters for the optimizer. Since SciPy's built-in default FD step is likely far from optimal for Hessian-vector products, the observed failure may be an artifact of an untuned configuration rather than evidence against finite-difference Hessians in general. This concern aligns directly with the reader's weakest assumption, and the concrete test proposed here would settle it: if a tuned centered-difference Hessian-vector product makes Newton-CG (FD) converge on the failing benchmarks, the central conclusion about FD instability would need to be weakened. The reader's CONDITIONAL verdict already accounts for this uncertainty, so no verdict change is needed. I also note a secondary inconsistency in the benchmarks: Section 4.4 recommends Newton-CG (FD) over Newton-CG (AD) for shape optimization, yet Section 4.5 and the abstract imply exact Hessians accelerate shape optimization; this should be corrected in a revision regardless of the FD tuning question.","tokens_in":29467,"tokens_out":9533,"duration_ms":97110,"concrete_test":"Re-run the source-field identification benchmark (Section 4.1) and the traction-force identification benchmark (Section 4.2) with Newton-CG, supplying hessp that evaluates the centered-difference Hessian-vector product of Eq. (28) with a user-controlled step size, sweeping h over 1e-3, 1e-4, 1e-5, and 1e-6, and with the same stopping tolerances as the Newton-CG (AD) runs. If either problem converges to the same objective value as Newton-CG (AD) for at least one h, the claim that finite-difference Hessians are inherently unreliable would be falsified. Report the chosen h, the gradient tolerance, and the convergence history for each run, matching the reporting level used for Newton-CG (AD).","verdict_should_be":"UNCHANGED","load_bearing_attack":"The benchmark-derived conclusion that finite-difference Hessian-vector products are unreliable and that exact AD Hessians are needed for robust Newton-CG convergence (Sections 4.1, 4.2, and 4.5) rests on an uncontrolled comparison. In Section 3.3.1, the authors validate their AD Hessian-vector product against a centered finite-difference approximation of the implicit gradient (Eq. 28) using hand-picked step sizes h = 10^-4 to 10^-1, and they report agreement near the 10^-7 level at the smallest step. This shows that finite-difference Hessian-vector products are not inherently inaccurate when the step size is chosen appropriately. Yet in Section 4.1 and 4.2, Newton-CG (FD) is reported as failing to converge without any specification of the finite-difference step size, the gradient tolerance, or whether the FD Hessian-vector product uses Eq. (28) or SciPy's internal approximation. SciPy's Newton-CG, when no hessp is supplied, defaults to a relative step of about sqrt(eps) ~ 1.5e-8 for float64, which is two orders of magnitude smaller than the smallest h validated in Section 3.3.1 and well below the optimal h ~ eps^(1/3) ~ 6e-6 for a second-derivative estimate. At such small step sizes, round-off error and solver tolerance in the gradient evaluation can easily dominate the Hessian-vector product, which would explain the observed non-convergence. Without a tuned-FD comparison, the central claim that exact Hessians avoid 'the instability and inaccuracy of finite difference Hessian approximations' is not supported. Furthermore, Section 4.4 reports Newton-CG (FD) as the fastest method in shape optimization, and Section 4.5 claims Newton-CG methods converge faster for shape optimization 'because of the precise information of second-order derivatives'—a statement that does not distinguish AD from FD and is in tension with the FD-instability narrative.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript derives an implicit Hessian-vector product for PDE-constrained optimization in a discretize-then-optimize framework, implements it in JAX-FEM using JVP/VJP primitives, and validates the implementation against finite-difference approximations and a Taylor remainder test. Four benchmark inverse problems are used to compare L-BFGS-B, Newton-CG with finite-difference Hessian-vector products (Newton-CG (FD)), and Newton-CG with the proposed implicit-differentiation Hessian-vector products (Newton-CG (AD)). The paper concludes that Newton-CG (AD) is robust and accelerates convergence for nonlinear problems, while L-BFGS-B is sufficient for linear problems, and that finite-difference Hessian approximations are unstable and inaccurate.","tokens_in":29773,"tokens_out":8503,"duration_ms":86535,"significance":"The core methodological contribution is sound and useful: the derivation in Section 2.3 is algebraically consistent, the algorithm in Section 3.1 is a natural and efficient composition of first-order implicit differentiation with JVP/VJP primitives, and the Taylor remainder test in Figure 5 provides independent, quantitative evidence that the Hessian-vector product is computed correctly to the expected third-order convergence. The public code is a strength, and the paper fills a real gap in differentiable FEM by making curvature information available as a primitive. The practical benchmark claims are only partially supported, however, because the comparison with finite-difference Hessians is not controlled and several benchmark settings are unreported. With those gaps fixed, the paper would make a solid contribution to the differentiable-physics community.","major_comments":[{"comment":"The claim that exact Hessians avoid \"the instability and inaccuracy of finite difference Hessian approximations\" is not established by the reported experiments, because Newton-CG (FD) is run with SciPy's default finite-difference Hessian-vector product and no step size, gradient tolerance, or other tuning is reported. The authors' own validation in Section 3.3.1 (Figures 3 and 4) shows that the finite-difference Hessian-vector product is strongly step-size dependent: the relative difference ranges from about 1e-3 at h=0.1 to about 1e-7 at h=1e-4. SciPy's default step for Newton-CG without a supplied Hessian is typically of order sqrt(machine epsilon) ~ 1.5e-8, which is two orders of magnitude smaller than the smallest h validated here and far below the usual optimum for a second-derivative estimate. The non-convergence of Newton-CG (FD) in Sections 4.1 and 4.2 may therefore be an artifact of the default step rather than of finite-difference Hessians in general. Please add a tuned-FD comparison, for example with h in the range 1e-4 to 1e-3 using Equation (28), and report all optimizer parameters; if tuned FD still fails to converge, the robustness claim is supported.","section":"Sections 4.1, 4.2, and 4.5"},{"comment":"The benchmark protocol is under-specified and not fully reproducible: the regularization weights α introduced in Equations (39), (43), and (48) are never assigned numerical values, the optimizer tolerances, initial guesses, and linear solver settings are not reported, all convergence histories come from single runs, and the shape-optimization design mapping from the rotation angles to the density field in Section 4.4 is not defined. Consequently, the quantitative timing comparisons and the Section 4.5 recommendations (\"Newton-CG faster for nonlinear, L-BFGS-B for linear\") are not yet supported by the evidence as presented. Please report the complete settings for each benchmark and, where feasible, include multiple restarts or perturbation studies to show that the qualitative conclusions are robust.","section":"Section 4 (Eqs. (39), (43), (48), (49))"}],"minor_comments":[{"comment":"The text repeatedly says \"PDF-constrained optimization problem\" where it should say \"PDE-constrained optimization problem.\"","section":"Section 4 (multiple places)"},{"comment":"The x-axis tick labels are corrupted in the typeset version (for example, \"0⊿000960\" instead of a numeric value with a power-of-ten notation); the figures should be regenerated with proper mathematical notation.","section":"Figures 3 and 4"},{"comment":"The step-size sweep stops at h=10^-4; adding values near the expected optimum for a second-derivative estimate, roughly h ~ eps^(1/3), would make the accuracy validation more convincing.","section":"Section 3.3.1"},{"comment":"The derivation implicitly assumes that the bordered KKT matrix, in particular ∂r/∂y, is nonsingular at every evaluation point; this implicit-function-theorem hypothesis should be stated explicitly as an assumption on the admissible parameter set.","section":"Section 2.3, Eq. (15)"},{"comment":"The mapping from the rotation angles θ to the density field ρ(x) through the sigmoid function is not defined; for reproducibility, the manuscript should give the exact formula and any smoothing or length-scale parameters.","section":"Section 4.4"},{"comment":"Steps 1 and 2 require solving a nonlinear forward problem and a linear adjoint problem, respectively; the intended numerical solver and convergence criterion for these solves should be specified.","section":"Algorithm 1"}],"recommendation":"major_revision","confidential_remarks":"The derivation and the Taylor remainder verification are solid, and the code availability is a clear strength. The main gate for acceptance is the uncontrolled finite-difference comparison in the benchmarks: if a tuned-FD Newton-CG comparison is added and the benchmark settings are fully reported, the paper's central practical claims can be properly evaluated. I would also encourage the editor to ask for the missing regularization values and solver settings, as these are essential for reproducibility in a computational journal."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The core of this paper is solid: the Hessian-vector product follows the standard bordered-Hessian derivation, but the concrete JAX implementation and the Taylor remainder test make it a usable deliverable for the JAX-FEM community. The Taylor test in Fig. 5 is the strongest evidence, as the third-order convergence of the second-order expansion validates the combined gradient and Hessian-vector product against objective evaluations, independent of any finite-difference comparison. The profiling of three AD-mode compositions is also systematic and useful, and the code is linked. Credit where it is due: someone who needs implicit second-order derivatives inside a finite element code now has a working recipe and a clear sense of when it pays off.\n\nThe stress-test note is right. Section 3.3.1 shows that finite-difference Hessian-vector products agree with the AD version to around 1e-7 when the step size h is chosen around 1e-4. Yet the benchmark section reports Newton-CG (FD) failing to converge in two of four examples without specifying the step size, the tolerance, or whether SciPy's internal finite-difference approximation was used. SciPy's default relative step is around sqrt(eps) ~ 1.5e-8 for float64, which is far too small for a second-derivative estimate and explains the failures. Without a tuned finite-difference run, the claim that exact Hessians avoid 'the instability and inaccuracy of finite difference Hessian approximations' is not established. The paper's own efficiency analysis even shows that FD and AD have similar cost, so the missing FD results are not a simple cost argument.\n\nThere are a few smaller soft spots. The regularization weight alpha is never reported in the benchmarks, which weakens reproducibility of the objective functions. The paper cites dolfin-adjoint but does not acknowledge that it already computes discrete adjoint Hessians; that is a framing weakness, not a fatal one. The nonsingularity of the Jacobian of the residual is assumed silently in Eq. (15) and the surrounding derivation; standard in the adjoint literature, but worth stating explicitly.\n\nWho is this paper for? Anyone building differentiable FEM solvers or using JAX-FEM to run second-order optimization. The algorithmic core is sound and the implementation is a real contribution. The benchmark flaws are fixable and do not invalidate the central machinery. I would send this to peer review, and I would recommend the author revise the benchmarks to report the finite-difference settings, tune the FD step size, give alpha values, and soften the claims about FD instability. That revision would make this a dependable reference for the subfield.","headline":"A correct and useful Hessian-vector product primitive for implicit FEM, with benchmark claims that overreach because the finite-difference comparison is uncontrolled.","tokens_in":734,"tokens_out":799,"would_cite":true,"duration_ms":32922,"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":"This paper derives a Hessian-vector product algorithm for implicit finite-element differentiation and shows that exact Hessians rescue Newton-CG where finite-difference Hessians diverge.","keywords":["implicit differentiation","Hessian-vector product","PDE-constrained optimization","differentiable physics","finite element method","automatic differentiation","Newton-CG","inverse problems"],"falsifier":"Run the two benchmarks where Newton-CG (FD) diverges (source-field identification and traction identification) with a systematically tuned finite-difference step size, including central differences over a range of $h$ and gradient scaling. If a tuned finite-difference Hessian-vector product converges in comparable time, the paper's central contrast between exact and approximate Hessians collapses.","tokens_in":29237,"feed_emoji":"📐","tokens_out":7153,"duration_ms":70139,"temperature":0.7,"pith_summary":"This paper addresses a gap in differentiable physics: first-order gradients of implicit finite-element solvers are standard, but second-derivative (Hessian) information has been missing. It derives and implements a Hessian-vector product for PDE-constrained optimization where the state is defined implicitly by a residual equation $r(y,\\theta)=0$, using only primitive automatic-differentiation operations. Accuracy is verified by comparison with finite differences and by a Taylor remainder test showing third-order convergence of the second-order expansion. On four benchmarks, Newton-CG with the exact Hessian converges reliably for nonlinear inverse problems, while a finite-difference Hessian version fails on two of them.","feed_headline":"Exact Hessians rescue Newton-CG on nonlinear inverse problems","feed_subtitle":"The paper derives Hessian-vector products for implicit finite-element solvers and tests them on four benchmarks.","key_machinery":"The load-bearing object is the bordered system of Eq. (15), formed from the Lagrangian $L(y,\\lambda,\\theta)=g(y,\\theta)+\\lambda^T r(y,\\theta)$; its solution gives $\\partial y/\\partial\\theta$ and $\\partial\\lambda/\\partial\\theta$. Algorithm 1 rewrites that system as two linear solves for the incremental state $\\hat{y}$ and incremental adjoint $\\hat{\\lambda}$, then evaluates Eq. (22) to produce $H\\hat{\\theta}$. The essential mechanism is that second-order terms such as $(\\partial/\\partial\\theta)(\\lambda^T\\partial r/\\partial y)\\hat{\\theta}$ can be interpreted as nested Jacobian-vector and vector-Jacobian products, making the Hessian available solely from AD primitives rather than from hand-derived formulas.","core_discovery":"The central claim is that the Hessian $d^2g/d\\theta^2$ of an objective subject to an implicit finite-element constraint can be computed exactly, direction by direction, without forming the full matrix. The Hessian-vector product $H\\hat{\\theta}$ is assembled from four solves\\textemdash forward state, adjoint, incremental forward, and incremental adjoint\\textemdash combined in Eq. (22); each solve is linear even when the forward problem is nonlinear. All second-order terms in the assembly are compositions of Jacobian-vector and vector-Jacobian products, so a differentiable-programming library with those primitives can supply curvature automatically. The Taylor remainder test confirms that the residual of the second-order Taylor expansion decays as $\\epsilon^3$, consistent with a machine-precision Hessian-vector product.","pith_inferences":["Beyond the paper: because the derivation only assumes an implicit residual equation, the same Hessian-vector product recipe transfers to non-finite-element implicit solvers.","Beyond the paper: exact curvature enables Hessian-based uncertainty quantification, such as Laplace approximations of posterior covariances for PDE-constrained inverse problems.","Beyond the paper: the profiling results suggest that the choice of AD composition mode matters, and that exploiting sparsity in the second-order terms is a natural next step toward much larger problem sizes."],"forward_implications":["On the nonlinear traction identification and shape-optimization benchmarks, Newton-CG with exact Hessian-vector products converges in wall-clock time faster than L-BFGS-B.","On the linear source-identification and thermal-mechanical control benchmarks, L-BFGS-B converges faster, so the extra curvature is not always worth its cost.","Finite-difference Hessian-vector products fail to converge in two of the four benchmarks, so exact implicit Hessians buy robustness, not just speed.","Because the implicit Hessian approach solves each incremental problem once per direction and can cache the state and adjoint, repeated Hessian-vector products at a fixed parameter point are cheaper than repeated finite-difference gradient pairs."],"supporting_citations":[{"why":"Provides the first-order implicit differentiation framework that this paper extends to second order.","marker":"[12]"},{"why":"Supplies the differentiable finite-element engine and sparse assembly used to build the Hessian experiments.","marker":"[13]"},{"why":"Supplies the Jacobian-vector and vector-Jacobian product primitives from which the Hessian-vector product is composed.","marker":"[3]"},{"why":"Defines the Newton-CG method whose convergence is tested with exact and finite-difference Hessians.","marker":"[18]"},{"why":"Supplies the Hessian-vector-product-based Newton minimization variant used in the optimizer comparison.","marker":"[19]"},{"why":"Defines the L-BFGS-B optimizer that serves as the first-order baseline.","marker":"[31]"},{"why":"Provides the optimizer implementations and the finite-difference Hessian-vector product used as the approximate baseline.","marker":"[32]"},{"why":"Introduces the Taylor remainder test used to validate the Hessian-vector product accuracy.","marker":"[26]"}],"fun_headline_variants":["Implicit Hessians accelerate Newton-CG on nonlinear FE problems","Exact curvature via Hessian-vector products in finite-element simulation","Second-order implicit differentiation for PDE-constrained optimization","Newton-CG with exact Hessians: faster convergence in inverse FE problems","Hessian-vector products from primitive AD: implicit FE differentiability"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The empirical conclusion that finite-difference Hessians are unusable rests on the assumption that the observed failures are inherent to finite differences rather than artifacts of an untuned step size; the derivation also assumes the bordered system matrix $\\partial r/\\partial y$ is nonsingular at every evaluation point.","fun_headline_variants_meta":{"raw":{"variants":["Implicit Hessians accelerate Newton-CG on nonlinear FE problems","Exact curvature via Hessian-vector products in finite-element simulation","Second-order implicit differentiation for PDE-constrained optimization","Newton-CG with exact Hessians: faster convergence in inverse FE problems","Hessian-vector products from primitive AD: implicit FE differentiability"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000201,"raw_usage":{"total_tokens":1348,"prompt_tokens":886,"completion_tokens":462,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":502,"completion_tokens_details":{"reasoning_tokens":374}},"tokens_in":502,"tokens_out":462,"duration_ms":4886,"temperature":1.0,"reasoning_tokens":374,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T20:29:59.280947+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the two benchmarks where Newton-CG (FD) diverges (source-field identification and traction identification) with a systematically tuned finite-difference step size, including central differences over a range of $h$ and gradient scaling. If a tuned finite-difference Hessian-vector product converges in comparable time, the paper's central contrast between exact and approximate Hessians collapses.","supporting_citations":[{"cited_title":"Efficient and modular implicit differentia- tion","cited_arxiv_id":null,"evidence_quote":"Provides the first-order implicit differentiation framework that this paper extends to second order."},{"cited_title":"Jax-fem: A differentiable gpu-accelerated 3d finite element solver for automatic inverse design and mechanistic data science","cited_arxiv_id":null,"evidence_quote":"Supplies the differentiable finite-element engine and sparse assembly used to build the Hessian experiments."},{"cited_title":"Function minimization by conjugate gradients","cited_arxiv_id":null,"evidence_quote":"Defines the Newton-CG method whose convergence is tested with exact and finite-difference Hessians."},{"cited_title":"Newton-type minimization via the lanczos method","cited_arxiv_id":null,"evidence_quote":"Supplies the Hessian-vector-product-based Newton minimization variant used in the optimizer comparison."},{"cited_title":"dolfin-adjoint 2018.1: automated ad- joints for fenics and firedrake","cited_arxiv_id":null,"evidence_quote":"Introduces the Taylor remainder test used to validate the Hessian-vector product accuracy."}],"review_version":1}