{"id":"f9b6618e-548f-4b7f-8da4-ed2d7105cc66","arxiv_id":"2507.09652","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"Polynomial regression with 512-bit arithmetic reaches machine-precision forecasting of low-dimensional chaotic systems from noise-free data, far exceeding previous valid prediction times.","lead":"The paper shows that a simple curve-fitting method, polynomial regression run at very high numerical precision, can forecast the Lorenz-63 chaotic system about as well as a numerical solver that knows the exact equations. This result suggests that today's noise-free benchmarks for learning chaotic dynamics are too easy, so the field can move on to harder, more realistic problems.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The record VPTs are measured against a fixed 512-bit RK4 discrete map that also generates the training data; the paper never tests whether the learned propagator tracks a different integrator or a more accurate reference solution.","rationale":"The reader's weakest_assumption identifies the same load-bearing issue: the benchmark uses one fixed RK4 trajectory as ground truth, and both training and testing come from that same numerical scheme. This is the most important vulnerability because every quantitative headline in Table 2 is measured against that discrete map. The paper's own Appendix G supports the concern: at dt=dt0, the learned target is exactly a degree-8 polynomial, and the 322-Lyapunov-time result is essentially exact recovery of the RK4 map. If the target were changed to a different integrator or a more accurate reference solution, the learned map would be different, and the reported VPT values could change materially. The paper does not include any such cross-integrator or refinement check, so the claim that noise-free low-dimensional chaotic forecasting is 'a solved problem' is conditional on the specific numerical surrogate. The reader also flagged test-set selection over degree, dt, and normalization, which is a real secondary concern, but the ground-truth issue is deeper and is the one that most directly threatens the central claim. The proposed experiment with DOP853 or a smaller RK4 step would settle whether the method is tracking the continuous system or merely reproducing one integrator's discrete map. Since the concern is real but not yet demonstrated to falsify the empirical results, the appropriate verdict remains conditional; the reader's CONDITIONAL verdict is unchanged.","tokens_in":65145,"tokens_out":6957,"duration_ms":83777,"concrete_test":"Regenerate the L63 experiments of Table 2 and Figure 4 using a high-order adaptive integrator (e.g., DOP853 with rtol/atol=1e-12) or RK4 with dt0=2^-12 as the only change; keep 64-bit storage, the same training sizes n, sampling steps dt, and the same grid of polynomial degrees. Record the best VPT0.5 (mean over 100 repetitions, selected as in the paper) relative to this new reference trajectory. If the method still reaches roughly 35.6 Lyapunov times and exceeds the 64-bit solver baseline for that reference, the integrator-relative concern is resolved. If the best VPT drops materially (for example, below the 64-bit RK4 baseline or below about 20 Lyapunov times), the headline claim is specific to the RK4 map and should be reframed.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central quantitative claim (Abstract; Section 3.1, Table 2) is that a polynomial propagator trained on noise-free data exceeds the accuracy of a standard 64-bit RK4 solver run from the true equations. But the benchmark target is not the ODE flow: all trajectories are generated by the authors' fixed RK4 solver at dt0=2^-10 for L63/L96 and dt0=2^-6 for TCSA, and the same solver output is used for training and for VPT evaluation (Section 2.4, Table 6). What is actually learned is the discrete map RK4^k, not the propagator of the continuous-time system. Appendix G makes this explicit: for L63 with dt=dt0, the target is exactly a degree-8 polynomial, and Section 3.2 reports more than 322 Lyapunov times in that setting 'because one RK4 solver step of the L63 system amounts to a polynomial propagator of degree 8'. That is a statement about recovering the integrator's map, not about machine-precision forecasting of the continuous chaotic system. Nothing in the paper varies the integrator or the solver step for the target while keeping data precision fixed, so the headline values of 35.6 and 104.9 Lyapunov times could be consequences of matching this particular RK4 discretization error pattern. A different high-accuracy reference trajectory would change the discrete map and, potentially, the measured VPT considerably.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a system-agnostic forecasting method that fits high-degree polynomial features to the discrete-time propagator of a chaotic system using ordinary least squares, executed in 512-bit arithmetic. On Lorenz-63, it reports valid prediction times of about 35.6 Lyapunov times in the default setup and up to 104.9 (or 322 at the solver step size) with 512-bit data, far exceeding prior ML results and slightly exceeding the 64-bit RK4 solver baseline. The method is also applied to Thomas' Cyclically Symmetric Attractor and Lorenz-96 in dimensions 5 through 9, where it matches the respective solver baselines. All training and test trajectories are generated by the authors' fixed RK4 integrator, and the VPT is evaluated against that same discrete reference.","tokens_in":65414,"tokens_out":7293,"duration_ms":78295,"significance":"If the claims are read as \"the method reproduces the data-generating discrete map at machine precision,\" the results are substantial: the empirical design is thorough, with 100 repetitions per configuration, confidence intervals, a Zenodo data release, and a formal proof in Appendix G that k RK4 steps of L63 are a polynomial of Fibonacci degree. The method is simple, data-efficient, and makes falsifiable predictions that prior ML benchmarks consistently fail to reach. However, the current abstract and title claim accuracy relative to \"the true underlying dynamical systems,\" which is not established because the ground truth is always the output of one fixed RK4 integrator. This overstatement materially tempers the significance and needs to be addressed.","major_comments":[{"comment":"The load-bearing claim that the method \"exceeds the accuracy of standard 64-bit numerical ODE solvers of the true underlying dynamical systems\" is not supported by the experiments, because the ground truth is always the output of the authors' fixed RK4 integrator with dt0=2^-10 (L63/L96) or 2^-6 (TCSA), and the same solver output is used for training and for VPT evaluation (Section 2.4, Table 6). The polynomial propagator therefore fits the discrete RK4 map, not the continuous-time flow; Appendix G makes this explicit by proving that one RK4 step of L63 is exactly a degree-8 polynomial, and Section 3.2 reports more than 322 Lyapunov times in the dt=dt0 case \"because one RK4 solver step of the L63 system amounts to a polynomial propagator of degree 8.\" I recommend either reframing the abstract and title claims as matching or exceeding 64-bit solvers on the same discrete data-generating map, or adding an experiment that trains on the RK4 data and evaluates against a different high-accuracy reference (e.g., a much smaller-step 512-bit integration or a different integration scheme) to show that the learned propagator tracks the continuous flow itself.","section":"Abstract; §1; §2.4; §3.1"},{"comment":"The Lorenz-96 claims of \"single-precision machine accuracy\" and of results \"statistically indistinguishable\" from the solver share the same discrete-map limitation: both the training data and the RK4 reference forecast are generated by the same RK4 solver at the same dt0, so the comparison quantifies how well the polynomial propagator emulates the RK4 map given single-precision initial data, not how well it approximates the continuous-time Lorenz-96 flow. The equivalence test is also informal: \"statistically indistinguishable\" is defined as the solver's VPT lying inside the 95% confidence interval of the method, which with 100 repetitions and interval widths around ±0.8 Lyapunov times provides only weak evidence of equivalence; a formal test or a stated equivalence margin would be more appropriate.","section":"§3.4; Table 3"},{"comment":"The sentence \"chaotic systems can be predicted to arbitrary accuracy—and thus arbitrary valid times\" overreaches the evidence. All reported VPT values are for a fixed polynomial degree bound (p <= 16, or 25 for TCSA) and for a fixed data precision; the paper's own analysis (Figure 5 and Appendix C) shows the forecast error floor is set by the 64-bit data rounding, and no convergence result in p and n is proved. The claim should be narrowed to \"arbitrary accuracy in reproducing the discrete-time propagator, up to the precision of the stored data and within the chosen polynomial feature family.\"","section":"§4 (Discussion); §3.2 (Figure 7)"},{"comment":"The statement that the dt=dt0 result \"is not a concern for larger time steps\" is insufficient, because every VPT value in Table 2, including the headline 35.6 and 104.9 values, is measured against the same RK4-generated reference trajectory. The distinction between (a) exactly recovering the RK4 map when dt=dt0 and (b) approximating the map of the continuous-time system needs to be drawn in the methods and reflected in the interpretation of the headline results; otherwise readers will take \"machine precision\" to refer to the ODE flow rather than to the fixed numerical scheme.","section":"§3.2; Appendix G"}],"minor_comments":[{"comment":"The heading \"Layapunov Exponent\" is a typo and should read \"Lyapunov Exponent.\"","section":"Appendix A"},{"comment":"Several captions and labels contain stray spaces such as \"V alid\", \"T raining\", \"T est\", and \"T ruth\"; these should be cleaned up.","section":"Table 1; Figure 1; Appendix I captions"},{"comment":"The diagonal normalization is written with \\hat{\\sigma}^2_i inside a diagonal matrix, but the notation is ambiguous; please define \\hat{S} explicitly as diag(\\hat{\\sigma}^2_1, ..., \\hat{\\sigma}^2_d) and clarify how \\hat{S}^{-1/2} is computed.","section":"§2.4, Eq. (13)"},{"comment":"The plateau of constant error for about 8 Lyapunov times is an artifact of comparing a 512-bit forecast against 64-bit stored ground truth; the text should state this explicitly so readers do not attribute the plateau to the physics of the Lorenz system.","section":"§3.1, Figure 5"},{"comment":"The statement \"there is no noise to overfit to\" is too categorical, since finite-sample variance and multicollinearity of high-degree monomial bases are present even for noise-free data; the poor performance at small n and high degree visible in the Appendix I tables suggests these effects matter and deserve a sentence of qualification.","section":"§4.2"}],"recommendation":"major_revision","confidential_remarks":"The manuscript is likely to attract attention because of the strong headline claims. In my view the core empirical result is real, but the abstract and discussion overstate the connection to the actual ODE flow, and the discrete-map interpretation should be made central. I would encourage the editor to request the reframing and the additional integrator-robustness experiment described in Major Comment 1. I see no citation or novelty concerns; the prior-literature table is useful and the data release is a clear strength."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Plainly: this paper delivers a real result. An OLS polynomial regression with 512-bit arithmetic posts valid prediction times of about 36 and up to 105 Lyapunov times on Lorenz-63, against a prior best of 13. It also matches or exceeds a 64-bit RK4 solver on Thomas' attractor and Lorenz-96. That is new, and the empirical work is unusually solid: 100 repetitions per cell, confidence intervals, a Zenodo dataset, and an appendix that formally proves one RK4 step of L63 is a degree-8 polynomial. Credit where due.\n\nThe soft spots are real but not deformational. First, the ground truth is the output of a fixed 512-bit RK4 integrator at step 2^-10, and the same integrator generates the training data. So the learned map is the discrete propagator of that integrator, not the continuous flow. The paper concedes this in Appendix G, and it does not invalidate the comparison against a 64-bit RK4 solver, but it should temper claims about machine precision for the underlying system. Second, the headline VPTs are the best over polynomial degree, time step, and normalization, selected by looking at the test set. That inflates expected performance; a nested or hold-out selection would be fairer. The abstract's 'solved problem' is too strong beyond the low-dimensional, noise-free, integrator-relative setup, and the Discussion's own limitation paragraph says real measurements are noisy.\n\nWho this is for: anyone benchmarking data-driven forecasting of low-dimensional chaos. It resets the bar. I'd bring it to reading group and would cite it if I were comparing methods. It deserves a serious referee. With honest reframing of the claims and test-set-independent model selection, it would be a solid accept; in current form, conditional.","headline":"Record VPTs from a high-precision polynomial propagator, but the benchmark is relative to a fixed RK4 integrator and the best numbers are selected on the test set; still deserves a serious referee.","tokens_in":65973,"tokens_out":3392,"would_cite":true,"duration_ms":36916,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["37D45","37M10","65P20"],"pacs":["05.45.-a","05.45.Pq"],"model":"deepseek-v4-flash","headline":"Noise-free forecasting of low-dimensional chaos is effectively solved: polynomial regression in 512-bit arithmetic matches or beats 64-bit solvers that know the true equations, reaching 105 Lyapunov times on Lorenz-63.","keywords":["chaos forecasting","valid prediction time","polynomial regression","machine precision","multi-precision arithmetic","Lorenz-63","Lorenz-96","Thomas cyclically symmetric attractor"],"falsifier":"Generate training and test data from a qualitatively different high-accuracy reference, for example an adaptive-step solver or a Taylor-method integrator with error control, and rerun the default configuration: if the 32 to 105 Lyapunov-time valid predictions persist, they reflect the underlying dynamics, while a collapse toward the 13 Lyapunov-time regime would show the result is specific to imitating the RK4 map. A second, cheaper check is to inject deterministic rounding noise into the last stored bit of the 64-bit training data and measure the valid prediction time as a function of that noise level.","tokens_in":64914,"feed_emoji":"🌀","tokens_out":9741,"duration_ms":98052,"temperature":0.7,"pith_summary":"Chaotic forecasting is usually framed as hard because tiny errors double on the Lyapunov timescale. The paper claims that with perfect (noise-free) data the difficulty mostly evaporates: fitting the system's one-step propagator by ordinary least squares on high-degree monomials, executed in 512-bit arithmetic, yields forecasts that are limited only by how accurately the training data itself represents the true trajectory. On Lorenz-63 the method reaches valid prediction times of 32 to 105 Lyapunov times, far beyond the prior ceiling of 13, and matches or exceeds the forecast accuracy of RK4 solvers that are given the governing equations. The same recipe works for the non-polynomial Thomas attractor and for Lorenz-96 in dimensions 5 through 9. The authors conclude that noise-free low-dimensional chaotic benchmarks are effectively solved problems, and too easy for evaluating machine-learning dynamics methods.","feed_headline":"Data-only model beats true-equation solvers on chaos benchmarks","feed_subtitle":"A 512-bit polynomial fit forecasts Lorenz-63 up to 105 Lyapunov times, versus 13 for past methods.","key_machinery":"The load-bearing object is the polynomial propagator: the map $\\hat\\Phi_{\\Delta t}(s) = \\hat\\beta^\\top \\xi_p(s)$ obtained by ordinary least squares, where $\\xi_p$ stacks all monomials of degree $\\le p$ in the $d$ state coordinates and $\\hat\\beta = (X^\\top X)^{-1}X^\\top Y$ solves the normal equations. Its two supporting mechanisms are high-precision linear algebra (512-bit arithmetic) to handle the ill-conditioned design matrix, and an auxiliary structural fact: the exact propagator of the RK4 integrator applied $k$ times to the Lorenz-63 equations is a polynomial of degree $F_{4k+2}$, the $(4k+2)$-nd Fibonacci number, so a single RK4 step is exactly a degree-$8$ polynomial. That identity explains why moderate-degree fits saturate the solver's own accuracy when the data step equals the solver step, and why sub-sampling the solver output (two steps already give degree 55) makes the learning task intrinsically harder.","core_discovery":"The paper's central claim is that a system-agnostic forecaster, namely ordinary least squares regression of the propagator $u(t)\\mapsto u(t+\\Delta t)$ on the multivariate monomials of the current state, reaches the numerical accuracy of the RK4 integrator that generated the data, and in the best configurations slightly exceeds it. The achievement is not a clever statistical model but numerical hygiene: high-degree polynomial features make the normal equations extremely ill-conditioned, and the paper shows that 512-bit arithmetic with full data normalization removes that bottleneck, while a 32-bit or 64-bit pipeline caps out near 7 or 21 Lyapunov times. Because the training data is noise-free, regularization is counterproductive and OLS is the right estimator. In the default setup (512-bit solver, 64-bit stored data, 512-bit fitter) the learned propagator reaches a valid prediction time of 35.6 Lyapunov times on Lorenz-63 versus 32.1 for a 64-bit RK4 solver, and with 512 bits used end to end it reaches 105 Lyapunov times, a ceiling the authors attribute to their chosen degree bound rather than to the method.","pith_inferences":["If the claim holds, the relevant question shifts from 'can we learn the system' to 'where does the precision floor sit': the error budget is set by the least precise component of the data pipeline, so adding deterministic rounding noise at the level of the last stored bit should visibly cut the achievable valid prediction time, a cheap testable extension of the paper's own logic.","Because the fitted object is the discrete integrator map, the method's machine-precision result is relative to whatever solver generated the data; training on RK4 output and validating against a radically different high-accuracy integrator would separate 'learning the continuous dynamics' from 'imitating the solver', which the paper does not test.","The observed preference for a medium data step ($\\Delta t \\approx 2^{-5}$) over finer sampling suggests a trade-off between propagator simplicity and state-space exploration; sweeping $n$ beyond $2^{15}$ and degree beyond 16 could map where the 105-Lyapunov-time ceiling actually moves.","The authors' own caveat that dimensions of order 100 are computationally infeasible implies the method's success is tied to low dimension; applying it to moderately high-dimensional or noisy settings is where the 'solved problem' conclusion is most likely to break."],"forward_implications":["On Lorenz-63, previously published methods top out near 13 Lyapunov times; this method reaches 32 to 36 in realistic precision configurations and 105 when data and computation are both carried at 512-bit precision, averaged over 100 initial conditions.","The learned propagator outperforms a 64-bit RK4 solver of the true equations in the default setup (35.6 versus 32.1 Lyapunov times), so a data-only forecast can be more accurate than a standard solver that is handed the exact vector field.","The results transfer to the non-polynomial Thomas' Cyclically Symmetric Attractor, where degree-25 polynomials match a 512-bit solver at 38.0 Lyapunov times, showing the success does not depend on the vector field being polynomial.","For the spatiotemporally chaotic Lorenz-96 model in dimensions 5 through 9, the method reaches single-precision machine accuracy, with valid prediction times statistically indistinguishable from the reference solver's.","Forecasting from a single state (random test mode) performs as well as continuing from the training trajectory, because the estimated propagator depends only on the current state.","The authors draw the practical conclusion that noise-free Lorenz-63-style benchmarks are too easy to discriminate among methods, and that future work on learning dynamical systems should focus on noisy data."],"supporting_citations":[{"why":"Supplies the Lorenz-63 equations, the benchmark system on which the 32 to 105 Lyapunov-time results are measured.","marker":"[Lor63]"},{"why":"Provides Thomas' Cyclically Symmetric Attractor, the non-polynomial system used to show the method is not exploiting polynomial vector fields.","marker":"[Tho99]"},{"why":"Provides the Lorenz-96 model, the higher-dimensional test case with dimensions 5 through 9.","marker":"[Lor95]"},{"why":"Holds the prior state-of-the-art valid prediction time of 13 Lyapunov times for Lorenz-63, the baseline this method must beat, and exemplifies a non-system-agnostic alternative.","marker":"[BPK16]"},{"why":"Introduced polynomial-feature forecasting with Ridge regularization, the approach this paper extends by dropping regularization and raising the degree and the arithmetic precision.","marker":"[Gau+21]"},{"why":"The MPLAPACK library supplies the 512-bit linear-algebra routines that make the ill-conditioned normal equations solvable.","marker":"[Nak22]"},{"why":"The MPFR library provides the multi-precision arithmetic on which the 512-bit computations rest.","marker":"[Fou+07]"},{"why":"The Armadillo solver options carry the single- and double-precision variants and the normalization comparisons.","marker":"[SC25]"},{"why":"The predecessor study that introduced the polynomial-propagator idea at degree 6 in 64-bit precision, extended here to higher degrees and 512-bit arithmetic.","marker":"[Sch+25]"},{"why":"Provides the reference value of the Lorenz-63 largest Lyapunov exponent used to convert system time to Lyapunov times.","marker":"[Vis04]"}],"fun_headline_variants":["512-bit polynomial fit beats exact chaos equations","Data-only forecaster reaches 105 Lyapunov times on Lorenz-63","Polynomial regression matches 64-bit ODE solver precision","Chaos prediction from data hits machine precision ceiling"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The whole comparison treats the output of one particular numerical solver, run at one fixed time step, as the true system: the polynomial learns to imitate that solver's discrete map rather than the exact continuous-time equations, so the claimed machine-precision ceiling is relative to that solver.","fun_headline_variants_meta":{"raw":{"variants":["512-bit polynomial fit beats exact chaos equations","Data-only forecaster reaches 105 Lyapunov times on Lorenz-63","Polynomial regression matches 64-bit ODE solver precision","Chaos prediction from data hits machine precision ceiling"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000232,"raw_usage":{"total_tokens":1507,"prompt_tokens":980,"completion_tokens":527,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":596,"completion_tokens_details":{"reasoning_tokens":461}},"tokens_in":596,"tokens_out":527,"duration_ms":5684,"temperature":1.0,"reasoning_tokens":461,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T17:52:09.966560+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Generate training and test data from a qualitatively different high-accuracy reference, for example an adaptive-step solver or a Taylor-method integrator with error control, and rerun the default configuration: if the 32 to 105 Lyapunov-time valid predictions persist, they reflect the underlying dynamics, while a collapse toward the 13 Lyapunov-time regime would show the result is specific to imitating the RK4 map. A second, cheaper check is to inject deterministic rounding noise into the last stored bit of the 64-bit training data and measure the valid prediction time as a function of that noise level.","supporting_citations":[],"review_version":1}