Learning Closure of Dynamical Systems with Kernel Ridge Regression

Evan Habbershaw; John Harlim; Senwei Liang

math.NAOct 1, 2026

Learning Closure of Dynamical Systems with Kernel Ridge Regression

Authors: Evan Habbershaw, John Harlim, Senwei Liang

Organizations: Department of Mathematics, The Pennsylvania State University · Department of Mathematics and Institute for Computational and Data Sciences, The Pennsylvania State University · Department of Mathematics and Statistics, Texas Tech University

Abstract

We develop a closure modeling framework for identifying missing components of dynamical systems using Kernel Ridge Regression (KRR). The framework addresses two classes of closure problems: difference-equation closures arising in ODE and PDE settings, and algebraic closures arising from moment closure in kinetic equations. For the first class, we derive an error bound in an ODE setting that quantifies contributions from time integration, approximation of unresolved scales, and interpolation required to couple unresolved-scale effects to the resolved solver. Numerical experiments on the Lorenz-63 system and the Kuramoto-Sivashinsky equation demonstrate accurate long-horizon predictions and substantial improvements over an LSTM-based closure model. For the second class, we consider moment closure for a one-dimensional kinetic equation by modeling discrepancies between kinetic and macroscopic fluxes as a function of the resolved macroscopic variables. We compare global KRR models based on PCA coordinates with spatially local models. While the global model performs well for unimodal initial conditions, its accuracy deteriorates for bimodal initial conditions. Spatially local models with appropriate modeling inputs improve robustness and achieve higher predictive accuracy.

Figures & tables

Stage treatment	DM-KRR	RBF-KRR	LSTM
Heun / endpoint RK2	$4.1969\pm 0.0799$	$4.1783\pm 0.1165$	$4.4603\pm 1.3410$
Quadratic RK4	$5.9874\pm 0.0868$	$5.9242\pm 0.0980$	$4.4924\pm 1.2797$
Cubic RK4	$8.3815\pm 0.3321$	$7.9683\pm 0.4844$	$4.6008\pm 1.3046$
Quartic RK4	$8.4832\pm 0.4910$	$8.3477\pm 0.6612$	$4.6115\pm 1.3114$
Quintic RK4	$\mathbf{8.6505\pm 0.5239}$	$8.4533\pm 0.7001$	$4.6439\pm 1.2873$

Table 1: Mean test valid prediction time for the Lorenz–63 hybrid closure models with

N_{\textnormal{train}}=4096

. The KRR hyperparameters are selected using three rollout validations per candidate pair. Results are reported as the mean and standard deviation across independently trained-model test means.

$N_{\textnormal{train}}$	DM-KRR	RBF-KRR	LSTM
512	$8.1719\pm 0.4849$	$7.8863\pm 0.5426$	$4.2784\pm 1.3443$
1024	$7.4267\pm 1.3808$	$7.4753\pm 0.9587$	$4.2253\pm 1.3061$
2048	$8.2535\pm 0.4977$	$8.0014\pm 0.4216$	$4.2652\pm 1.2222$
4096	$\mathbf{8.3815\pm 0.3321}$	$7.9683\pm 0.4844$	$4.6008\pm 1.3046$

Table 2: Mean test valid prediction time for the Lorenz–63 hybrid models using cubic reconstruction of the closure at RK4 stage times. Results are reported as the mean and standard deviation across independently trained-model test means.

Panel of Figure 1 — Figure 1: Mean validation VPT (left) and a representative validation (right) over the diffusion-kernel hyperparameter grid for $N_{\textnormal{train}}=16384$ and $m=0$ , over ten training windows. The figure shows the mean landscape over ten independent training windows. The outlined square marker identifies the selected grid point.

Panel of Figure 2 — Figure 2: Left: Representative test rollout for the $N_{\textnormal{train}}=16384$ , $m=0$ , diffusion-kernel model. The signed physical-space difference between the full-order resolved field and the hybrid closure rollout is recorded. The dashed line marks the verified VPT, $2.608$ . Right: The time-averaged energy in the six retained Fourier modes is shown.

Figure 3: Mean test VPT as a function of training-set size for the $m=0$ KS closure models. KRR error bars denote variability across independently trained models. The LSTM values are included as a qualitative comparison because the LSTM training-data convention and noise augmentation are not identical to those used for KRR. The dashed reference line marks the reported long-memory LSTM mean, $2.494688$ , whose reported standard deviation is $0.95904$ .

Data type	Phase shift	Varying parameters
I	$\mu=\eta\Delta x$	$A,\ \sigma,\ \eta$
II	$\mu_{1,2}=\phi\pm\frac{d}{2}$	$A_{1},A_{2},\sigma_{1},\sigma_{2},$ $d$ , $\phi$

Table 3: The two physical initial condition data types used in the kinetic experiments.

Model	Training data	KRR centers	$\mathcal{E}_{1}$	$\mathcal{E}_{2}$
Global	I	$8208$	$0.001285\pm 0.000200$	$0.028925\pm 0.005158$
Global	I+II	$8208$	$0.001005\pm 0.000234$	$0.008638\pm 0.007220$
zeroth-order local	I	$4032$	$0.001455\pm 0.000252$	$0.006088\pm 0.001457$
zeroth-order local	I+II	$4032$	$0.001891\pm 0.000341$	$0.003178\pm 0.000675$
first-order local	$(\textnormal{I},\,\textnormal{I})$	$(4032,8064)$	$0.000665\pm 0.000096$	$0.008098\pm 0.002001$
first-order local	$(\textnormal{I},\,\textnormal{I+II})$	$(4032,8064)$	$0.000793\pm 0.000105$	$0.003768\pm 0.001161$

Table 4: Comparison of the kinetic closure models under different closure models and training-data compositions. The columns

\mathcal{E}_{1}

and

\mathcal{E}_{2}

report the average prediction errors for the unimodal and bimodal testing data, respectively (lower values are better). For the first-order local model, the ordered pair in the training-data column gives the trajectory types used to train

(\widetilde{\mathcal{R}}_{0},\widetilde{\mathcal{R}}_{1})

, and the number of training centers is reported as

(N_{0},N_{1})

. The zeroth-order model

\widetilde{\mathcal{R}}_{0}

is held fixed between the two first-order comparisons.

Figure 4: Final-time density profiles for a Type II validation trajectory. Each panel compares the kinetic reference with one macroscopic model: the zero-correction baseline, the global PCA–KRR closure, the zeroth-order local closure, and the first-order local residual closure.

Figure 5: Space–time evolution of temperature for the same Type II validation trajectory as in Figure 4 . Panel (a) shows the kinetic reference temperature

\theta

, while panels (b)–(d) show the absolute temperature errors, for the global model, the zeroth-order local, and the first-order local model, respectively. The three error panels use a common linear color scale so that their magnitudes are directly comparable. The global model can reproduce the final macroscopic profile reasonably well while accumulating substantially larger structured transient errors during recursive rollout; both local models substantially reduce the structured rollout error and track the kinetic dynamics more faithfully throughout the trajectory.

Appendix figures & tables6 assets

Supplementary material from the paper’s appendix.

Appendix

$N_{\mathrm{train}}$	$m=0$		$m=1$
	RBF	Diffusion	RBF	Diffusion
2048	$1.0171\pm 0.1177$	$0.9121\pm 0.1183$	$0.9885\pm 0.1562$	$0.9161\pm 0.1567$
4096	$1.1241\pm 0.1884$	$1.1101\pm 0.1780$	$1.2156\pm 0.1276$	$1.1494\pm 0.1948$
8192	$1.6903\pm 0.3016$	$1.8164\pm 0.3150$	$1.6484\pm 0.3067$	$1.8659\pm 0.2142$
16384	$2.3650\pm 0.4693$	$\mathbf{2.6516\pm 0.3819}$	$2.1400\pm 0.4798$	$2.6490\pm 0.4311$

Appendix

Table 5: Mean test VPT for the KRR closure models with memory depths

m=0

and

m=1

. For each configuration, the reported statistics are the mean and standard deviation across ten independently trained-model test means, with

50

test rollouts per trained model.

$N_{\textnormal{train}}+2N_{V}$	Memory $m$	VPTs (ETDRK4)
$2.5\times 10^{5}$	19	$2.494688\pm 0.95904$
2048+4000	0	$1.41372390\pm 0.6669579$
4096+4000	0	$1.4588309\pm 0.71275799$
8192+4000	0	$1.5381616\pm 0.725299612$
16384+4000	0	$1.7413667\pm 0.8070518383$
2048+4000	1	$1.3993662\pm 0.6346926$

Appendix

Table 6: Test VPTs of true 6 modes and predicted 6 modes for the LSTM model using ETDRK4.

Representation	Type I	Type I+II (Mixed)
Global	$72\times 114=8208$	$72\times 57+72\times 57=8208$
Zeroth-order local	$72\times 56=4032$	$72\times 28+72\times 28=4032$
First-order local	$72\times 112=8064$	$72\times 56+72\times 56=8064$

Appendix

Table 7: KRR training-pair allocations used in the kinetic-model comparisons. In each table entry, the number

72

represents the number of Type I and Type II trajectories used for training. For the global model, the second factor denotes the number of temporal full-field snapshots selected per trajectory; for the local models, it denotes the number of local space–time interface samples selected per trajectory.

Figure 6: Final-time velocity profiles for the same bimodal rollout used in the main-text kinetic comparison in Figure 4 . Each panel compares the kinetic reference with one macroscopic model: the zero-correction baseline, the global closure, the zeroth-order local closure, and the first-order local residual closure.

Figure 7: Final-time temperature profiles for the same bimodal rollout used in the main-text kinetic comparison in Figure 4 . Each panel compares the kinetic reference with one macroscopic model, using the same model ordering as Figure 6 .

Figure 8: Final-time energy profiles for the same bimodal rollout used in the main-text kinetic comparison in Figure 4 . Each panel compares the kinetic reference with one macroscopic model, using the same model ordering as Figure 6 .

PDE Surrogate Modeling

PDE surrogate modeling builds fast approximations of partial differential equation solutions for simulation, prediction, and design.

See all 322 papers

Learning dynamical systems from noisy data with Weak-form Kernel Ridge Regression
Jun 30, 2026Max Kreider, John Harlim, Daning HuangKernel Ridge Regression Kernel Regression
Department of Mathematics The Pennsylvania State University, University Park, PA 16802, USA · Department of Mathematics, Institute for Computational and Data Sciences The Pennsylvania State University, University Park, PA 16802, USA · Department of Aerospace Engineering The Pennsylvania State University, University Park, PA 16802, USA
Mamba-Assisted Non-Markovian Closure for Reduced-Order Modeling
Jun 3, 2026Zhi-Feng Wei, Saad Qadeer, Panos StinisReduced-Order Modeling
Advanced Computing, Mathematics, and Data Division, Pacific Northwest National Laboratory, Richland, WA 99354, USA · Department of Applied Mathematics, University of Washington, Seattle, WA 98195, USA · Division of Applied Mathematics, Brown University, Providence, RI 02912, USA
The Dynamic-Probabilistic Consistency Gap in Chaotic Surrogate Modeling
May 29, 2026Andre Herz, Matthijs Pals, Daniel Durstewitz +1Dynamical Systems Surrogate Modeling
1Interdisciplinary Center for Scientific Computing, Heidelberg University, Germany · Faculty of Mathematics and Computer Science, Heidelberg University, Germany · 3Dept. of Theoretical Neuroscience, Central Institute of Mental Health (CIMH), Mannheim, Germany +3

Jun 30, 2026J/K move · Enter open · S save

Learning dynamical systems from noisy data with Weak-form Kernel Ridge Regression

Max Kreider, John Harlim, Daning Huang

Department of Mathematics The Pennsylvania State University, University Park, PA 16802, USA · Department of Mathematics, Institute for Computational and Data Sciences The Pennsylvania State University, University Park, PA 16802, USA · Department of Aerospace Engineering The Pennsylvania State University, University Park, PA 16802, USA

Kernel Ridge Regression Kernel Regression Nonlinear System Identification Chaotic Dynamical Systems Scientific ML

Open paper AskPDF arXiv

Accurate prediction of complex dynamical systems from noisy measurements remains a significant challenge in scientific computing. Kernel ridge regression learning strategies are often effective when applied to clean data, but have limited success with noisy data. Recent work has observed that a weak formulation can act to filter noisy data, and different learning strategies have achieved increased noise robustness with a weak-form framework. In this manuscript, we give an overview of the filtering mechanism behind the weak formulation and provide a bias-variance error decomposition. Using these insights, we combine a weak formulation with a kernel learning strategy to propose Weak-form Kernel Ridge Regression (WKRR) for learning dynamical systems. The proposed framework is simple to implement, effective for both clean and noisy data, and outperforms several baseline methods. We demonstrate the performance of WKRR on chaotic benchmark systems in up to 64 dimensions, as well as 15,000-dimensional real-world fluid data.

Abstract

PDE Surrogate Modeling

Explore similar work

Learning dynamical systems from noisy data with Weak-form Kernel Ridge Regression

Mamba-Assisted Non-Markovian Closure for Reduced-Order Modeling

The Dynamic-Probabilistic Consistency Gap in Chaotic Surrogate Modeling