Physics is the Best Teacher: Consistency Learning for Time-Invariant Operators of Chaotic Dynamics

Lufang Chiang; Jiachen Yao; Thomas Y. L. Lin; Anima Anandkumar

cs.LGOct 2, 2026

Physics is the Best Teacher: Consistency Learning for Time-Invariant Operators of Chaotic Dynamics

Authors: Lufang Chiang, Jiachen Yao, Thomas Y. L. Lin, Anima Anandkumar

Organizations: National Taiwan University · California Institute of Technology · University of Washington

Abstract

Accelerating the prediction of long-term behavior in chaotic systems is crucial in scientific computing. However, existing methods rely on numerical solvers or autoregressive models that advance one small step at a time, which makes long horizons expensive. We instead view this problem as learning the system's time-invariant evolution operator, which jumps the state across a large time span in a single evaluation. To this end, we derive the consistency equations a time-invariant operator must satisfy, with differential and compositional objectives in physical time. These equations also connect the learned operator to the physics-prescribed instant dynamics, enabling physics embedding in consistency learning. Across five chaotic systems, we find that physics-distilled consistency makes both short-term trajectories and long-term statistics more accurate. The learned operator survives temporal extrapolation and requires one-tenth as many evaluations as autoregressive rollout, offering an efficient route to long-term simulation of chaotic dynamics.

Figures & tables

Figure 1: Overview. (a) Autoregressive models and neural ODEs advance in small steps and need many NFEs. (b) We propose consistency conditions for a time-invariant operator that jumps any $\Delta t$ in one step. (c) The model itself or the governing equations supply the velocity in these conditions.

Evolution operator
		in distribution					out of distribution
Method	NFE	$\mathcal{E}_{L_{2}}$	$D_{\mathrm{stsp}}$	$D_{\mathrm{H}}$	$\mathcal{E}_{D_{2}}$	$\mathcal{E}_{\mathrm{trans}}$	$\mathcal{E}_{L_{2}}$	$D_{\mathrm{stsp}}$	$D_{\mathrm{H}}$	$\mathcal{E}_{D_{2}}$	$\mathcal{E}_{\mathrm{trans}}$
Data	100	0.46	0.20	3.1	2.5	1.7	0.44	0.18	2.9	1.6	1.1
Semigroup	100	0.46	0.21	3.4	3.3	1.9	0.44	0.18	2.9	1.3	1.3
Eulerian-self	100	0.45	0.22	3.5	3.8	1.6	0.43	0.17	2.6	1.6	1.4
Lagrangian-self	100	0.45	0.21	3.4	3.2	1.5	0.48	0.18	2.9	2.0	0.94

Table 1: Average prediction results over the five systems. Five metrics averaged over three seeds and over the five systems by geometric mean, with

D_{\mathrm{H}}

and

\mathcal{E}_{D_{2}}

10^{-2}

and

\mathcal{E}_{\mathrm{trans}}

10^{-3}

. Out-of-distribution means the model is trained on only half of the horizon and rolled out fully. Shading marks the physics-distilled rows, and bold the best method in each column of a block. Sections D and E define the metrics and the training setup, and Section G gives the table system by system.

Evolution operator
		Lorenz		Chen		Rössler		Thomas		Halvorsen
Method	NFE	in-dist	OOD	in-dist	OOD	in-dist	OOD	in-dist	OOD	in-dist	OOD
Data	100	0.32 $\pm$ 0.05	0.25 $\pm$ 0.04	0.26 $\pm$ 0.03	0.27 $\pm$ 0.03	0.53 $\pm$ 0.06	0.53 $\pm$ 0.02	0.90 $\pm$ 0.04	0.92 $\pm$ 0.04	0.56 $\pm$ 0.09	0.52 $\pm$ 0.01
Semigroup	100	0.30 $\pm$ 0.06	0.24 $\pm$ 0.03	0.27 $\pm$ 0.02	0.29 $\pm$ 0.04	0.54 $\pm$ 0.05	0.54 $\pm$ 0.01	0.88 $\pm$ 0.09	0.90 $\pm$ 0.07	0.55 $\pm$ 0.09	0.52 $\pm$ 0.03
Eulerian-self	100	0.27 $\pm$ 0.06	0.24 $\pm$ 0.01	0.27 $\pm$ 0.03	0.30 $\pm$ 0.04	0.56 $\pm$ 0.02	0.50 $\pm$ 0.04	0.89 $\pm$ 0.07	0.83 $\pm$ 0.05	0.50 $\pm$ 0.02	0.49 $\pm$ 0.02
Lagrangian-self	100	0.22 $\pm$ 0.03	0.31 $\pm$ 0.14	0.31 $\pm$ 0.05	0.32 $\pm$ 0.06	0.56 $\pm$ 0.06	0.54 $\pm$ 0.03	0.85 $\pm$ 0.15	0.87 $\pm$ 0.04	0.55 $\pm$ 0.07	0.56 $\pm$ 0.02

Table 2: Five-Lyapunov-time accuracy.

\mathcal{E}_{L_{2}}(5T_{\lambda})

on each of the five systems, mean

\pm

std over three seeds, with each method at its native number of evaluations, stated in the NFE column. In-distribution means every model is trained on the full horizon, and out-of-distribution means it is trained on the first half and rolled out to the full horizon. Shading marks the physics-distilled objectives and bold the best method within each block. Section G gives the same table at one Lyapunov time.

Method	Equation and distinguishing feature	NFE	Avg. $\mathcal{E}_{L_{2}}\!\downarrow$ ID / OOD
Evolution operator: elapsed-time step $h$
Data	$\mathcal{L}=\mathbb{E}\!\left[\big\\|\bar{v}_{\theta}(u_{t},h)-\tfrac{u_{t+h}-u_{t}}{h}\big\\|^{2}\right]$ . Paired-trajectory MSE ( $h>0$ ).	100	0.46 / 0.44
Semigroup	$\psi_{\theta}(u,h)=\psi_{\theta}\big(\psi_{\theta}(u,\alpha h),(1-\alpha)h\big)$ . One jump agrees with two composed jumps.	100	0.46 / 0.44
Eulerian-self	$\bar{v}_{\theta}(u,h)=\bar{v}_{\theta}(u,0)+hD\bar{v}_{\theta}(u,h)[\bar{v}_{\theta}(u,0)]-h\partial_{h}\bar{v}_{\theta}(u,h)$ . Model velocity at the starting state $u$ .	100	0.45 / 0.43
Lagrangian-self	$\bar{v}_{\theta}(u,h)=\bar{v}_{\theta}\big(u+h\bar{v}_{\theta}(u,h),0\big)-h\partial_{h}\bar{v}_{\theta}(u,h)$ . Model velocity at the predicted endpoint.	100	0.45 / 0.48
Eulerian-phys	$\bar{v}_{\theta}(u,h)=\mathcal{F}(u)+hD\bar{v}_{\theta}(u,h)[\mathcal{F}(u)]-h\partial_{h}\bar{v}_{\theta}(u,h)$ . Governing-equation velocity at the starting state $u$ .	100	0.28 / 0.27

Table 3: Method definitions for the benchmark tables. The evolution-operator rows give the training objective or consistency condition; the baseline rows give the prediction rule. NFE is the number of network evaluations used to reach five Lyapunov times in either evaluation regime. The final column reports the average relative

L_{2}

error at five Lyapunov times from Table 1 . Physics-distilled rows are shaded, and bold marks the lowest error within each group and evaluation regime.

System	Parameters	$\lambda_{\max}$	Source
Lorenz	$(\sigma,\rho,\beta)=(10,\,28,\,8/3)$	$0.9056$	Sprott [6]
Chen	$(a,b,c)=(35,\,3,\,28)$	$2.0272$	Noorani et al. [30] , Eftekhari and Jafari [31]
Rössler	$(a,b,c)=(0.2,\,0.2,\,5.7)$	$0.0714$	Sprott [6]
Thomas	$b=0.1$	$0.055$	Sprott and Chlouverakis [32]
Halvorsen	$a=1.3$	$0.6928$	Sprott [33]

Table 4: The five systems. Parameters, the largest Lyapunov exponent

\lambda_{\max}

, and its source.

Hyperparameter	Value
Training
Optimizer	Adam ( $\beta_{1}=0.9$ , $\beta_{2}=0.99$ ) [ 48 ]
LR scheduler	constant
Batch size	$2048$
Iterations	$200{,}000$
Validation interval	$5000$ iterations

Table 5: Training settings. Fixed settings above, and the ranges searched by Vizier below.

Figure 4: Sampling distribution and largest step.

\mathcal{E}_{L_{2}}(5T_{\lambda})

on Lorenz at

p=2

, averaged over three seeds, for the three distributions of the training step, and each panel pair splits by

\Delta t_{\max}

Figure 5: Adaptive weight coefficient.

\mathcal{E}_{L_{2}}(5T_{\lambda})

on Lorenz under the Beta distribution at

\Delta t_{\max}=0.05T_{\lambda}

, averaged over three seeds. The data-only model has no consistency term, so it is drawn once.

		Time (ms)
Method	NFE	Lorenz	Chen	Rössler	Thomas	Halvorsen
AR $v$ -pred	1000	404.9 $\pm$ 1.5	396.9 $\pm$ 1.3	398.4 $\pm$ 1.7	401.0 $\pm$ 1.1	401.5 $\pm$ 1.1
CFO	4000	1668.7 $\pm$ 4.9	1636.4 $\pm$ 4.3	1642.2 $\pm$ 2.7	1650.3 $\pm$ 4.0	1665.5 $\pm$ 3.2
Evolution operator	100	45.1 $\pm$ 0.1	44.3 $\pm$ 0.2	44.0 $\pm$ 0.1	45.0 $\pm$ 0.8	44.8 $\pm$ 0.2

Table 6: Wall-clock time. Time to roll out one trajectory over five Lyapunov times on each of the five systems with a single NVIDIA RTX 4090, mean

\pm

std over

100

test trajectories, with each method at its native NFE. Shading marks the evolution operator and bold the shortest time in each column.

Evolution operator
		Lorenz		Chen		Rössler		Thomas		Halvorsen
Method	NFE	in-dist	OOD	in-dist	OOD	in-dist	OOD	in-dist	OOD	in-dist	OOD
Data	20	0.99 $\pm$ 0.27	0.71 $\pm$ 0.10	0.66 $\pm$ 0.03	0.64 $\pm$ 0.09	3.1 $\pm$ 0.8	2.9 $\pm$ 0.3	2.8 $\pm$ 0.5	2.5 $\pm$ 0.4	1.4 $\pm$ 0.4	1.3 $\pm$ 0.1
Semigroup	20	0.97 $\pm$ 0.17	0.63 $\pm$ 0.09	0.70 $\pm$ 0.04	0.67 $\pm$ 0.10	3.3 $\pm$ 0.8	3.1 $\pm$ 0.9	3.1 $\pm$ 1.2	2.6 $\pm$ 0.3	1.5 $\pm$ 0.4	1.2 $\pm$ 0.2
Eulerian-self	20	0.91 $\pm$ 0.39	0.73 $\pm$ 0.35	0.64 $\pm$ 0.13	0.74 $\pm$ 0.22	3.4 $\pm$ 0.3	2.5 $\pm$ 0.6	2.5 $\pm$ 0.3	2.2 $\pm$ 0.1	1.2 $\pm$ 0.2	1.0 $\pm$ 0.1
Lagrangian-self	20	0.71 $\pm$ 0.19	0.83 $\pm$ 0.19	1.1 $\pm$ 0.4	0.86 $\pm$ 0.24	3.5 $\pm$ 1.1	2.8 $\pm$ 0.5	2.4 $\pm$ 0.7	2.4 $\pm$ 0.3	1.5 $\pm$ 0.7	1.5 $\pm$ 0.2

Table 7: One-Lyapunov-time accuracy.

\mathcal{E}_{L_{2}}(1T_{\lambda})

\times 10^{-2}

, mean

\pm

std over three seeds, with each method at its native step size, stated in the NFE column. Shading marks the physics-distilled block and bold the lowest error in each column of a block.

Evolution operator
		Thomas					Halvorsen
Method	NFE	$\mathcal{E}_{L_{2}}$	$D_{\mathrm{stsp}}$	$D_{\mathrm{H}}$	$\mathcal{E}_{D_{2}}$	$\mathcal{E}_{\mathrm{trans}}$	$\mathcal{E}_{L_{2}}$	$D_{\mathrm{stsp}}$	$D_{\mathrm{H}}$	$\mathcal{E}_{D_{2}}$	$\mathcal{E}_{\mathrm{trans}}$
Data	100	0.90	0.14	5.1	2.2	1.8	0.56	0.048	3.3	2.2	0.89
Semigroup	100	0.88	0.13	4.5	2.1	2.4	0.55	0.050	3.7	3.9	0.74
Eulerian-self	100	0.89	0.16	6.4	27	3.2	0.50	0.046	3.2	2.6	0.86
Lagrangian-self	100	0.85	0.15	5.6	1.7	2.6	0.55	0.049	3.8	6.0	1.1

Table 8: In-distribution results. Five metrics on five chaotic systems, averaged over three seeds, with

D_{\mathrm{H}}

and

\mathcal{E}_{D_{2}}

10^{-2}

and

\mathcal{E}_{\mathrm{trans}}

10^{-3}

. Shading marks the physics-distilled rows and bold the best method in each column of a block. Sections D and E define the metrics and the training setup.

Evolution operator
		Thomas					Halvorsen
Method	NFE	$\mathcal{E}_{L_{2}}$	$D_{\mathrm{stsp}}$	$D_{\mathrm{H}}$	$\mathcal{E}_{D_{2}}$	$\mathcal{E}_{\mathrm{trans}}$	$\mathcal{E}_{L_{2}}$	$D_{\mathrm{stsp}}$	$D_{\mathrm{H}}$	$\mathcal{E}_{D_{2}}$	$\mathcal{E}_{\mathrm{trans}}$
Data	100	0.92	0.13	4.3	1.1	0.88	0.52	0.047	2.2	1.7	0.90
Semigroup	100	0.90	0.14	5.2	0.90	2.1	0.52	0.047	2.5	1.1	1.1
Eulerian-self	100	0.83	0.13	4.8	1.2	1.9	0.49	0.048	3.0	3.1	0.53
Lagrangian-self	100	0.87	0.14	4.1	0.77	1.5	0.56	0.049	4.0	5.8	0.72

Table 9: Out-of-distribution results. The same five metrics on the five systems, averaged over three seeds, with

D_{\mathrm{H}}

and

\mathcal{E}_{D_{2}}

10^{-2}

and

\mathcal{E}_{\mathrm{trans}}

10^{-3}

. Every model is trained on the first half of the horizon and rolled out to all of it. Shading marks the physics-distilled rows and bold the best method in each column of a block. Sections D and E define the metrics and the training setup.

Figure 6: Error curves.

\mathcal{E}_{L_{2}}

against time on the five systems, averaged over three seeds with one-standard-deviation bands, and the dotted line is the relative error

0.3

. Top: in distribution. Bottom: out of distribution, with the half beyond the training window shaded.

Self-distillation
	Lorenz		Chen		Rössler		Thomas		Halvorsen
Consistency term	in-dist	OOD	in-dist	OOD	in-dist	OOD	in-dist	OOD	in-dist	OOD
Eulerian	0.27 $\pm$ 0.06	0.24 $\pm$ 0.01	0.27 $\pm$ 0.03	0.30 $\pm$ 0.04	0.56 $\pm$ 0.02	0.50 $\pm$ 0.04	0.89 $\pm$ 0.07	0.83 $\pm$ 0.05	0.50 $\pm$ 0.02	0.49 $\pm$ 0.02
Lagrangian	0.22 $\pm$ 0.03	0.31 $\pm$ 0.14	0.31 $\pm$ 0.05	0.32 $\pm$ 0.06	0.56 $\pm$ 0.06	0.54 $\pm$ 0.03	0.85 $\pm$ 0.15	0.87 $\pm$ 0.04	0.55 $\pm$ 0.07	0.56 $\pm$ 0.02
Eulerian $+$ Lagrangian	0.21 $\pm$ 0.04	0.31 $\pm$ 0.13	0.27 $\pm$ 0.04	0.27 $\pm$ 0.03	0.55 $\pm$ 0.03	0.50 $\pm$ 0.02	0.81 $\pm$ 0.05	0.84 $\pm$ 0.07	0.49 $\pm$ 0.03	0.46 $\pm$ 0.02
Physics-distillation

Table 10: Combined differential identities.

\mathcal{E}_{L_{2}}(5T_{\lambda})

on each of the five systems, mean

\pm

std over three seeds, for each differential identity on its own and for the combination of Equation 52 . Shading marks the physics-distilled block and bold the lowest error in each column of a block.

Figure 7: Long-term statistics on all five systems. One column per method on seed $0$ , showing the attractor cloud of the $50T_{\lambda}$ free run with the truth in gray behind it, and the power spectrum of one coordinate on one trajectory. The frequency $f$ is in units of the Lyapunov rate $\lambda_{\max}$ .

Chaotic Time-Series Prediction

Chaotic time-series prediction studies how to forecast nonlinear systems whose sensitive dependence on initial conditions limits long-range accuracy.

See all 54 papers

Conservation Buys Stability and Factoring Buys Counterfactuals in Physical World Models
Sep 17, 2026Yufeng Wang, Parivesh Priye, Lu Wei +1Neural Surrogate Modeling Dynamical Systems
Stony Brook University · Georgia Institute of Technology · Westlake University
Symbolic Neural ODEs: Learning interpretable models from time-series data
Aug 22, 2026Nibodh Boddupalli, Jeff MoehlisNonlinear System Identification System Identification
Dept. of Mechanical Engineering, University of California, Santa Barbara, CA 93106
Physics-Integrated Operator Learning via Gaussian Splatting Representations
Aug 25, 2026Jihao Zhang, Junyi Guo, Jian-Xun WangPDE Operator Learning PDE Surrogate Modeling
Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY, USA

Sep 17, 2026J/K move · Enter open · S save

Conservation Buys Stability and Factoring Buys Counterfactuals in Physical World Models

Yufeng Wang, Parivesh Priye, Lu Wei, Haibin Ling

Stony Brook University · Georgia Institute of Technology · Westlake University

Neural Surrogate Modeling Dynamical Systems World Models Counterfactual Prediction Structure-Preserving Neural Networks

Open paper AskPDF arXiv

A learned simulator can reproduce its training conditions accurately yet fail in two distinct ways once those conditions change. Over long rollouts, small errors accumulate until the trajectory drifts away from physically plausible behavior; under an intervention on a physical parameter, the model may continue to follow the law seen during training rather than the intervened one. We show that these two failures require different structural remedies. Evolving a learned energy with a symplectic integrator preserves the geometry of the conservative dynamics and keeps rollouts bounded and physically meaningful for up to

100\times

the training horizon, while equal-capacity predictors, an energy-regularized predictor, and a tuned neural ODE diverge. By contrast, encoding the physical coupling through an explicit linear factorization enables the model to follow a never-seen sign of that coupling, whereas an unrestricted parameterization remains locked to the training law. Crucially, the two mechanisms are separable: removing the structure responsible for long-horizon stability leaves counterfactual transfer intact, while removing the factorized coupling destroys counterfactual transfer without eliminating stability. This double dissociation, established with matched controls that remove or replace one structural component at a time, persists beyond the headline three-body system and remains visible when the physical state must be inferred from pixels rather than provided directly. The result is a concrete design principle for physical world models: long-horizon stability and changed-law generalization arise from distinct structural commitments, and each can be imposed deliberately without requiring the other.

Abstract

Chaotic Time-Series Prediction

Explore similar work

Conservation Buys Stability and Factoring Buys Counterfactuals in Physical World Models

Symbolic Neural ODEs: Learning interpretable models from time-series data

Physics-Integrated Operator Learning via Gaussian Splatting Representations