{"id":"72e2dc4c-59d6-444a-bfd4-121e50ddc7a7","arxiv_id":"2607.17577","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":2,"one_line_summary":"A new open-source harmonic-balance solver applies automatic differentiation only to reduced localized nonlinear forces, enabling efficient periodic-response and stability analysis of high-dimensional nonlinear systems without user-derived Jacobians.","lead":"pyHB is an open-source solver that uses automatic differentiation inside the harmonic balance method to find periodic responses of nonlinear vibration systems without requiring users to manually derive Jacobians. It can trace stable and unstable solution branches on high-dimensional models, reporting a several-hundred-fold speedup over direct time integration in a 2000-degree-of-freedom benchmark.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Scalability and the 'general' claim rest on Eq. (47) localization (q,r<<n); no evidence or caveat that the AD/GPU-memory advantages survive distributed nonlinearities.","rationale":"The paper delivers a working, well-engineered tool with real validation on four examples and a reproducible open-source implementation. The central technical contribution, localized AD for HB Jacobians, is sound and clearly demonstrated for the cases presented. The reader's conditional verdict already marks the main risk: the strongest claims of generality and memory efficiency depend on the localization assumption in Eq. (47), which is not flagged as a limitation. I agree that this is the most load-bearing concern. My proposed test directly checks whether the AD/GPU-memory advantage degrades when q,r ~ n. If it does, the abstract's 'general' language is too strong, but the core method and its value for localized-nonlinearity systems remain intact. Thus the reader's CONDITIONAL verdict is appropriate, and my stress-test does not change it. Secondary issues such as the theoretical status of AD subgradients for non-smooth forces are real but less central to the claimed scalability result; they would also be worth addressing, but the localization premise is the one on which the headline numbers depend.","tokens_in":36790,"tokens_out":5764,"duration_ms":53275,"concrete_test":"Re-run the 2000-DOF Bernoulli beam benchmark with the same n=2000 and m=5, 20, 50, but replace the single local nonlinear spring at the free end with a nonlinear force acting on every DOF (q=r=n), keeping all other settings identical. If per-point time or GPU memory scales with n (e.g., GPU peak no longer flat at 243.5 MB) or the run exceeds 16 GB, the localized-nonlinearity condition in Eq. (47) is essential to the headline scalability, and the 'general' claim must be qualified.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The headline performance and memory numbers (Table 1: 637.8 MB RAM / 243.5 MB GPU at 202,000 HB unknowns; 0.44 s/point) are achieved by applying AD only to the reduced nonlinear force f_nl^(r) in Eq. (47), with q, r << n. This is explicit in Sections 3.4 and 3.5: the computational graph 'depends mainly on the local dimensions q and r' and the GPU peak stays 243-265 MB across all examples. Nothing in the paper establishes that the same efficiency holds when nonlinearities are distributed, i.e. when q and r are O(n). In that case the Jacobian block H_nl^(r) in R^{(2p+1) x r x q} and the per-sample AD graph scale with n, so the 'controllable GPU memory' and 'linear scaling in unknowns' claims collapse. The abstract nevertheless advertises 'general user-defined nonlinear systems' and a 'one-stop benchmark platform for HB-based nonlinear dynamics analysis'; Section 6 lists quasi-periodic and DAE extensions as future work but does not flag this localization requirement as a scope limitation. The gap between the demonstrated capability (localized nonlinearities) and the claimed generality is the load-bearing weakness. This is a scope-limitation issue rather than an internal inconsistency, but it is central because the abstract's strongest quantitative claims are conditional on it.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper presents pyHB, an open-source Python implementation of the Harmonic Balance (HB) method for periodic response analysis of nonlinear dynamical systems. The central methodological contribution is the use of PyTorch automatic differentiation (AD) to compute the Jacobian of the nonlinear force, applied only to a reduced, localized representation of that force (Eq. (47), Sections 3.4 and 3.5). The authors supplement this with FFT-based residual projection, sparse assembly of the linear Jacobian part, a blocked solver for the arc-length continuation system, weighted arc-length scaling, and Floquet stability analysis via an implicit trapezoidal integrator. Four numerical examples are presented: a quasi-zero-stiffness isolator, a piezoelectric-magnetic energy harvester, a 284-DOF aeroengine dual-rotor model, and a 2000-DOF Bernoulli beam with a localized nonlinear spring. The results are verified against explicit Runge-Kutta or Newmark-β time integration for stable branches, and unstable branches are classified via computed Floquet multipliers. The headline performance claim is that the AD-enhanced solver reaches about 0.44 s per continuation point at 202,000 HB unknowns with 637.8 MB additional RAM and 243.5 MB GPU memory, while HB-AD (the authors' prior method) runs out of memory.","tokens_in":37101,"tokens_out":4644,"duration_ms":45083,"significance":"If the claims hold, pyHB would be a useful and reproducible engineering tool: it removes the need for manual derivation of nonlinear-force Jacobians in HB analysis, provides a modular workflow from model definition through stability assessment, and demonstrates good scalability on a genuinely high-dimensional example (2000 DOFs, 202,000 HB unknowns). The comparison with the authors' previous HB-AD implementation is informative and shows a substantive improvement in GPU memory footprint. The open-source release, the inclusion of executable examples, and the verification against independent time-integration baselines are concrete strengths that support the credibility of the reported results. The main significance risk is that the generality and scalability claims are conditioned on a localization assumption that is not stated as a limitation and is not tested in the distributed-nonlinearity regime.","major_comments":[{"comment":"The scalability and memory claims rest on the localization assumption f_nl = Ξ f_nl^(r)(Θx'', Θx', Θx, τ, ω) with q, r << n. The paper states in §3.5 that the AD computational graph 'depends mainly on the local dimensions q and r', and Table 1 shows GPU memory staying near 243 MB precisely because the nonlinearity is localized. For distributed nonlinearities with q, r = O(n), the reduced block H_nl^(r) and the AD graph scale with n, so the claimed 'controllable GPU memory' and near-linear scaling would not hold. This is not flagged as a scope limitation in §6; the abstract's 'general user-defined nonlinear systems' and 'one-stop benchmark platform' therefore overstate the demonstrated capability. A caveat and a small distributed-nonlinearity experiment or complexity estimate are needed.","section":"§3.4, Eq. (47), and §5/Table 1"},{"comment":"The paper advertises support for non-smooth nonlinear forces (gaps, piecewise stiffness, Heaviside functions) and asserts that PyTorch's subgradients provide a 'robust solution'. However, for Newton-HB iterations the relevant object is a consistent tangent of the projected residual; for discontinuous or non-differentiable forces a subgradient is not automatically a valid Newton derivative. The aeroengine example uses a Heaviside function and a fractional exponent, yet no convergence diagnostics, mesh-refinement study, or comparison of subgradient-based Jacobians against a smoothed approximation is reported. As it stands, the claim that non-smooth forces are handled 'without additional manual processing' is not supported by evidence. At minimum, the paper should report behavior at switching points and discuss conditions under which the AD subgradient yields a convergent HB iteration.","section":"§3.5 and §4.3 (Eqs. (62)-(63))"},{"comment":"The performance comparison with Newmark-β is presented as evidence of a 'several-hundred-fold speedup'. The comparison is per continuation point versus a 10-period Newmark trajectory initialized directly from the HB solution. This is a legitimate pointwise cost comparison, but it does not compare full workflows (e.g., obtaining a complete response curve with continuation versus a swept time-integration run), and the text in §4.4 even contrasts '15 s for the complete pyHB curve' with '>80 s for 10 periods at one ω'. The abstract's 'several-hundred-fold' number comes from the per-point ratio, which should be stated explicitly. I recommend clarifying the comparison basis in the text so that readers do not infer a total-workflow speedup of that magnitude.","section":"§4.4 and §5, Table 1"}],"minor_comments":[{"comment":"In the expression for ̄C(τ), the term ∂f_nl/∂x'' should presumably be ∂f_nl/∂x'. As written, the perturbation damping term contains a derivative with respect to x'', which appears twice and is inconsistent with the preceding linearization.","section":"Eq. (17)"},{"comment":"Several typographical errors should be corrected: 'baed' (p. 10), 'Simlarly' (p. 8), 'discribed' (p. 7), 'anbitrary' (p. 9), 'utlized' (p. 17), 'omｅga' (p. 17), and inconsistent spellings of 'preprocessing/preparation'.","section":"Throughout"},{"comment":"The weighted arc-length scaling factors q_s and ω_s are introduced as user choices with no guidance or sensitivity study. Since the continuation behavior and even convergence can depend on these factors, an automatic or recommended scaling procedure would strengthen reproducibility.","section":"§3.6"},{"comment":"The frequency resolution parameter F_r is described as enabling subharmonic analysis, and Example 4.1 uses F_r = 1/3. The paper should explain more concretely how F_r, the period T, and the number of FFT samples N1 interact, and whether the harmonic basis is always commensurate.","section":"§3.9 and §4.1"},{"comment":"The comparison with COCO and NLvib cites the earlier computational environment in reference [128]. This is acceptable as a reference point, but the sentence 'the computational cost of both methods exceeds 24 hours' should be accompanied by the original hardware/software conditions to avoid implying a same-machine benchmark.","section":"§5"},{"comment":"The statement that 'the explicit Runge-Kutta method adopted in the previous examples is unstable for this system' is imprecise: DOP853 is not globally unstable; rather it would require excessively small steps for this stiff (contact-bearing) problem. A more careful wording would avoid an unfair characterization of explicit integrators.","section":"§1 and §4.3"}],"recommendation":"major_revision","confidential_remarks":"The manuscript is competent and the open-source contribution is valuable, but the advertised generality and the headline performance numbers are conditioned on the localized-nonlinearity assumption (Eq. (47)). I would ask the authors to temper the abstract and Section 6 claims, add an explicit limitation, and provide at least one experiment or complexity analysis for a case with distributed nonlinearity. The non-smooth AD claim also needs a concrete validation before the paper can be accepted as a general-purpose solver statement."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: this is a real tool, and the central engineering idea is a genuine improvement. The paper applies PyTorch-based AD only to the reduced nonlinear force of a localized nonlinearity, so the AD computational graph and GPU memory are set by the local DOF count (q,r) rather than the global DOF count (n). That is a clean, useful trick, and their own prior HB-AD code fails on the same problems (out of memory on the aeroengine model, ~1GB GPU on small ones), while pyHB sits around 250MB. The code is open source and modular, and the four examples are validated against independent time integration — DOP853 and Newmark-β. The unstable branches, the 1/3 subharmonic response, the combination resonances in the dual-rotor system, and the mixed first/second-order harvester are all checkable results. That is credible, reproducible evidence.\n\nThe main soft spot is real and central: everything — the small memory, the scaling, the 0.44 s/point number — rests on the localization assumption in Eq. (47). If nonlinearities couple many DOFs, q and r are O(n), the AD graph and the H_nl tensors grow with system size, and the claimed memory and speed advantages collapse. The paper states the assumption when it derives the localization, but the abstract's 'general user-defined nonlinear systems' and 'one-stop benchmark platform' do not flag that the headline efficiency only holds for localized nonlinearities, and the future-work discussion doesn't list it as a scope limit. That is a scope overclaim, not an internal inconsistency. Since most engineering nonlinearities really are local (bearings, joints, supports), the tool is still highly useful; it just needs an honest scope statement near the front.\n\nThe other issues are smaller. The 'several-hundred-fold speedup' compares a steady-state solver against a transient integrator. The protocol is disclosed fairly (10 periods, initialized from the HB solution), so this is not dishonest, but the headline number is not an apples-to-apples benchmark. The use of PyTorch subgradients for non-smooth forces (Heaviside, clearance) is pragmatic and works in the tested cases, but there is no theory for Newton convergence at non-smooth points; that deserves a caveat. And there is a typo in Eq. (17): the damping Jacobian repeats ∂f_nl/∂x'' instead of ∂f_nl/∂x'.\n\nWho this is for: people computing periodic responses of mechanical or structural systems with local nonlinearities, and anyone building HB software who wants a reference baseline. It deserves a serious referee. I would send it, treating the scope-claim revision and the comparison-framing fix as conditions, not as reasons to reject.","headline":"A real and well-tested HB tool whose headline numbers come from localized nonlinearities; the 'general/one-stop' wrapper overclaims that scope, but the core is worth serious engagement.","tokens_in":37599,"tokens_out":6798,"would_cite":true,"duration_ms":59681,"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":"A solver that needs only the system equations, not hand-derived Jacobians, maps full periodic-response landscapes—stable and unstable—of million-scale vibration problems.","keywords":["harmonic balance","automatic differentiation","nonlinear dynamics","arc-length continuation","periodic orbit stability","localized nonlinearities","open-source solver","sparse linear algebra"],"falsifier":"Run pyHB on a system with n=2000 where the nonlinear force couples every DOF to every other (e.g., a fully populated cubic stiffness matrix) with m=50 harmonics. If per-iteration time or GPU memory scales with n rather than staying near 0.44 s and 250 MB, the central scalability claim fails. Reproducing the reported 2000-DOF beam benchmark and measuring total wall time for a full branch sweep would also confirm or refute the several-hundred-fold speedup.","tokens_in":36657,"feed_emoji":"⚡","tokens_out":6397,"duration_ms":56361,"temperature":0.7,"pith_summary":"pyHB claims that automatic differentiation, applied only to the small nonlinear-force block of a large mechanical system, removes the main obstacle to using harmonic balance analysis on realistic engineering models: the manual derivation of Jacobian matrices. By exploiting the fact that most nonlinearities are localized at joints, bearings, or supports, the solver keeps the automatic-differentiation computational graph independent of global degrees of freedom, so memory stays flat even as the number of harmonics grows. The paper demonstrates the complete workflow—periodic-response solution, arc-length continuation, and periodic-orbit stability—on four systems, ending with a 2000-degree-of-freedom beam with 202,000 harmonic unknowns solved in about 0.44 seconds per continuation point with a few hundred megabytes of extra memory. If the localization premise holds for a given model, the user supplies only the equations of motion; pyHB returns the full stable and unstable response landscape, including subharmonic and combination resonances.","feed_headline":"0.44 s per continuation point on 202,000 unknowns","feed_subtitle":"Open-source tool computes full stable and unstable response branches from equations alone—no hand-derived Jacobians.","key_machinery":"The load-bearing object is the localized-nonlinearity decomposition f_nl = Ξ f_nl^(r)(Θx'', Θx', Θx, τ, ω) with q,r << n. It is what lets automatic differentiation operate on a low-dimensional computational graph; it also makes the nonlinear Jacobian a sparse scatter of a small dense tensor. Around it, three pieces do the heavy lifting: precomputed Fourier-basis tensors for the nonlinear Jacobian update; a Kronecker-structured linear Jacobian with sparse LU factorization; and a blocked solve of the arc-length-augmented linear system that reuses one sparse factorization per continuation step. Periodic-orbit stability is computed by propagating a monodromy matrix via an implicit trapezoidal in","core_discovery":"The central claim is that the entire harmonic balance workflow can be driven by user-defined equations alone. The key move is rewriting the nonlinear force as a selection-matrix product that restricts it to a few degrees of freedom, so automatic differentiation builds the required tangent matrices for this reduced force, and the Jacobian for the full system is assembled from precomputed linear blocks plus a sparse nonlinear contribution obtained by FFT-based harmonic projection. The paper reports that this makes a 2000-DOF discretized beam with 50 harmonics (about 202,000 unknowns) solvable at roughly 0.44 s per continuation point, several hundred times faster than stepping through ten perio","pith_inferences":["The localization assumption (Eq. 47) is untested for globally coupled nonlinearities; if a system has dense nonlinear coupling (e.g., a fully populated cubic stiffness matrix), the automatic-differentiation graph would grow with the global DOF count and the flat GPU-memory plateau would break. A direct test would run pyHB on an all-to-all coupled 2000-DOF chain and watch memory.","The reported speedup is against a 10-period time integration that already starts from the harmonic balance solution; for a full operating-domain sweep the time-integration baseline would likely be far slower, so the practical speedup may exceed the paper's stated number, but a fair comparison would report total wall-clock per complete branch, not just per continuation point.","The modular architecture suggests pyHB's localized Jacobian could serve as a preconditioner in matrix-free Newton-Krylov implementations, which the paper itself lists as future work but does not test.","The implicit-trapezoidal monodromy integration requires inverting a 2n-by-2n matrix per time segment; whether that remains efficient at lower sparsity or higher DOF counts is untested and could become a bottleneck for stability analysis of very large systems."],"forward_implications":["Users can change a nonlinear force model by editing one function, without re-deriving tangent stiffness or damping matrices, so design iterations in isolators, harvesters, and rotors become fast.","Full stable and unstable branches become routinely available for high-dimensional finite-element models, allowing prediction of jump phenomena, hysteresis, and hidden resonances that time integration misses.","The memory footprint of automatic-differentiation-enhanced harmonic balance no longer grows with the number of degrees of freedom or harmonics in the observed cases, making single-GPU analysis of 2000-DOF systems with 50 harmonics feasible.","A standard benchmark platform emerges: the same modular pipeline reproduces low-DOF strongly nonlinear isolator responses and high-DOF aeroengine and beam responses, facilitating comparison of future harmonic balance methods.","The automatically computed tangent blocks could feed other Newton variants, such as quasi-Newton or matrix-free implementations, as preconditioners, extending the approach toward systems with millions of unknowns."],"fun_headline_variants":["0.44s/point on 202k unknowns via AD harmonic balance","No user Jacobians needed: pyHB does the tangents automatically","Open-source pyHB: AD harmonic balance for any nonlinear system","From equations alone: pyHB computes full response branches"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"The paper assumes nonlinear forces act on only a few degrees of freedom (local contact, joint, or support), so the automatic-differentiation graph stays small; if a nonlinearity couples many degrees of freedom at once, the claimed speed and flat memory would not hold.","fun_headline_variants_meta":{"raw":{"variants":["0.44s/point on 202k unknowns via AD harmonic balance","No user Jacobians needed: pyHB does the tangents automatically","Open-source pyHB: AD harmonic balance for any nonlinear system","From equations alone: pyHB computes full response branches"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001262,"raw_usage":{"total_tokens":5055,"prompt_tokens":846,"completion_tokens":4209,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":590,"completion_tokens_details":{"reasoning_tokens":4136}},"tokens_in":590,"tokens_out":4209,"duration_ms":25577,"temperature":1.0,"reasoning_tokens":4136,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-01T17:33:53.262564+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run pyHB on a system with n=2000 where the nonlinear force couples every DOF to every other (e.g., a fully populated cubic stiffness matrix) with m=50 harmonics. If per-iteration time or GPU memory scales with n rather than staying near 0.44 s and 250 MB, the central scalability claim fails. Reproducing the reported 2000-DOF beam benchmark and measuring total wall time for a full branch sweep would also confirm or refute the several-hundred-fold speedup.","supporting_citations":[],"review_version":1}