Paper deep dive
On a joint simultaneous learning of relevant feature subsets and subspaces in regression-like problems
Illia Horenko
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 8/3/2026, 9:38:18 AM
Summary
The paper introduces Entropy-Optimal Manifold Regression (EOMR), an extension of Entropy-Optimal Manifold Clustering (EOMC) for joint simultaneous identification of relevant feature subsets and subspaces in nonstationary and nonlinear regression problems. EOMR achieves linearly-scaling iteration and memory complexities. It is benchmarked against state-of-the-art AI/ML tools (gradient boosted random forests, deep neural networks, TabPFN v.03) on chaotic dynamics problems: Lorenz-96 systems and Hasegawa-Wakatani tokamak plasma data. EOMR demonstrates superior predictive accuracy (lower RMSE) and lower model complexity compared to competitors. For the Hasegawa-Wakatani example, EOMR identifies a simple 8-parameter autoregressive process for the leading Essential Orthogonal Function (EOF) dynamics.
Entities (10)
Relation Signals (8)
Illia Horenko → affiliatedwith → RPTU Kaiserslautern-Landau
confidence 95% · Illia Horenko Chair for Mathematics of AI, Faculty of Mathematics, RPTU Kaiserslautern-Landau
Entropy-Optimal Manifold Regression → appliedto → Lorenz-96
confidence 95% · on predicting the Lorenz-96 systems dynamics in strongly- and very-strongly chaotic regimes
Entropy-Optimal Manifold Regression → appliedto → Hasegawa-Wakatani model
confidence 95% · on a data from the Hasegawa-Wakatani model on the edge of the tokamak plasma
Entropy-Optimal Manifold Regression → extends → Entropy-Optimal Manifold Clustering
confidence 95% · We extend a recently introduced Entropy-Optimal Manifold Clustering (EOMC) to allow for a joint simultaneous identification... that we coin as Entropy-Optimal Manifold Regression (EOMR)
Entropy-Optimal Manifold Regression → identifies → Essential Orthogonal Function
confidence 90% · For a Hasegawa-Wakatani example, EOMR distills a very simple entropy-optimal and skilful description of the leading Essential Orthogonal Function (EOF) dynamics
Entropy-Optimal Manifold Regression → outperforms → TabPFN v.03
confidence 90% · transformer-based AI tools like TabPFN v.03... result in orders of magnitude inferior root mean squared prediction errors... when compared to the EOMR
Entropy-Optimal Manifold Regression → outperforms → Gradient Boosted Random Forests
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We extend a recently introduced Entropy-Optimal Manifold Clustering (EOMC) to allow for a joint simultaneous identification of subsets and subspaces of relevant features in nonstationary and nonlinear regression problems. It is shown that the proposed extension - that we coin as Entropy-Optimal Manifold Regression (EOMR) - allows a robust learning with linearly-scaling iteration and memory complexities. EOMR is compared to the most complete set of state-of-the-art tools from the Artificial Intelligence (AI) and Machine Learning (ML) that is available to the author, on the very challenging problems from chaotic and fluid dynamics: (i) on predicting the Lorenz-96 systems dynamics in strongly- and very-strongly chaotic regimes (with forcing parameter being $F=8$ and $F=12$, respectively); and, (ii) on a data from the Hasegawa-Wakatani model on the edge of the tokamak plasma. It is demonstrated that the proposed benchmarks (i) and (ii), indeed, are the very challenging problems for the state of the art ML and AI tools - since both the general-purpose gradient boosted random forests and deep neuronal networks, as well as transformer-based AI tools like TabPFN v.03 (more spezialised for large-dimensional small data learning problems) - result in orders of magnitude inferior root mean squared prediction errors, and orders of magnitude larger model complexities, when compared to the EOMR. For a Hasegawa-Wakatani example, EOMR distills a very simple entropy-optimal and skilful description of the leading Essential Orthogonal Function (EOF) dynamics, given by linear, causal and weakly-stationary autoregressive process described by just 8 parameters.
Tags
Links
- Source: https://arxiv.org/abs/2607.28080v1
- Canonical: https://arxiv.org/abs/2607.28080v1
Trouble viewing inline? Open PDF directly →
Full Text
48,426 characters extracted from source content.
Expand or collapse full text
On a joint simultaneous learning of relevant feature subsets and subspaces in regression-like problems Illia Horenko Chair for Mathematics of AI, Faculty of Mathematics, RPTU Kaiserslautern-Landau, Gottlieb-Daimler-Str. 48, Kaiserslautern, 67663, Germany. Contributing authors: horenko@rptu.de; Abstract We extend a recently introduced Entropy-Optimal Manifold Clustering (EOMC) to allow for a joint simultaneous identification of subsets and subspaces of relevant features in nonstationary and nonlinear regression problems. It is shown that the proposed extension - that we coin as Entropy-Optimal Manifold Regression (EOMR) - allows a robust learning with linearly- scaling iteration and memory complexities. EOMR is compared to the most complete set of state-of-the-art tools from the Artificial Intelligence (AI) and Machine Learning (ML) that is available to the author, on the very challenging problems from chaotic and fluid dynamics: (i) on predicting the Lorenz-96 systems dynamics in strongly- and very-strongly chaotic regimes (with forcing parameter beingF = 8 andF = 12, respectively); and, (i) on a data from the Hasegawa-Wakatani model on the edge of the tokamak plasma. It is demonstrated that the proposed benchmarks (i) and (i), indeed, are the very challenging problems for the state of the art ML and AI tools - since both the general-purpose gradient boosted random forests and deep neuronal networks, as well as transformer-based AI tools like TabPFN v.03 (more spezialised for large-dimensional small data learning problems) - result in orders of magnitude inferior root mean squared prediction errors, and orders of magnitude larger model complexities, when compared to the EOMR. For a Hasegawa-Wakatani example, EOMR distills a very simple entropy-optimal and skilful description of the leading Essential Orthogonal Function (EOF) dynamics, given by linear, causal and weakly-stationary autoregressive process described by just 8 parameters. 1 Introduction Data-driven approaches and tools from Artificial Intelligence (AI) and Machine Learning (ML) are gaining increasing popularity in modelling and prediction of chaotic systems [1–7]. However, applications of common data-hungry ML and AI tools to such systems - especially in practical real-life situations - are characterized by a so-called small data learning challenge, when the feature 1 arXiv:2607.28080v1 [stat.ML] 30 Jul 2026 space dimension is relatively large, and the available statistics of systems observation is relatively small, for example, due to the expense of the direct numerical simulation (e.g., in fluid mechanics), or due to the relative shortness of available observation history (e.g., in weather and climate research) [6–11]. Moreover, when applying very recent tools based on Large Language Models (LLMs), main bottleneck is represented by the unfavourable computational and memory scalings of transformer-based architectures - that scale quadratically both in the feature dimension and in the data statistics size, making training and applications of such tools very expensive for large realistic problems [11, 12]. A key to success in such situations is to identify or to learn the important feature dimensions that are most relevant for the considered problem. Existent dimension reduction and feature extrac- tion methods can be subdivided into two groups: (i) into the methods identifying the relevant subsets of original feature dimensions, and (i) the methods that identify the relevant subspaces: • Subset methods aim at finding a relevant subset of original features, and include approaches, like the methods based on 1D statistical importances (e.g., based on filtering-out all of the feature dimensions below a certain p-value treshold), the very popular Ridge- and Lasso- regression methods (using l2 and l1 parameter norms as regularizers [13]), or entropic learning methods using Shannon entropy as regularizers for probability distributions of model parameters [8–10, 12]. • Subspace methods aim at finding a linear or a nonlinear transformation of the original coordinate system in the feature space, for example, by means of an appropriate rotation and projection on a low-dimensional manifold, like in the Principal Component Analysis (PCA) [14, 15], as well as in its numerous linear and nonlinear extensions [16–26]. To address numerical issues resulting from the polynomial scalability of subspace methods - making them exceedingly expensive for large-dimensional feature spaces - another large group of methods aims at combining these two approaches into pipelines, for example, by first pre-selecting a potentially-relevant subset of original features (e.g., via the p-value tresholding), followed by an application of the subspace-based method, in the way how it is implemented in supervised PCA, sparse PCA, and in the related algorithms [27–30]. 2 In the following, we will present an extension of the recently-introduced unsupervised and linearly-scalable Entropy-Optimal Manifold Clustering method (EOMC) for feature extraction beyond global linearity and stationarity assumptions - that belongs to a family of subspace methods. We will demonstrate that this extension - that we will call Entropy-Optimal Manifold Regression (EOMR), preserves the central advantages of EOMC (like applicability to nonlinear and nonstation- ary problems, linear scalability of iteration and memory complexities, as well as the interpretable probabilistic formulation that allows finding entropy-optimal models). Moreover, we will show that beyond these advantages of EOMC, EOMR allows a supervised identification of features relevant for linear and nonlinear regression and autoregression problems - by means of a simultaneous and joint learning of relevant feature subsets and subspaces, implemented as a numerical solution of the unified, probabilistic, and entropy-regularized optimization problem. After analyzing the numerical properties of the proposed solution, and discussing the way of hyperparameter selection and relations to existing approaches, we exemplify the application of EOMR to two challenging problems: to Lorenz-96 model in strongly- and very-strongly-chaotic regimes, as well as to the Hasegawa-Wakatani model from Magnetohydrodynamics. 2 Methods 2.1 Mathematical formulation of the Entropy-Optimal Manifold Regression (EOMR) Let X ∈R D,T be a given (D× T )-dimensional real-valued data matrix, with every column X(:,t) representing a D-dimensional vector of feature values for a data instance with an index t, where t = 1,...,T (: denotes a column-extraction operation). Let Y ∈R 1,T be a given series of values, for which we want to learn a functional relation that for any given instance t should map a vector X(:,t) to a scalar value Y (t) 1 . Let there be K ”local” mappings from X(:,t) to Y (t), each of them is characterized by the orthogonal d-dimensional linear manifold projector T k ∈R D,d , k = 1,...,K, 1 Without a loose of generalaty, we consider a scalar valued regression problem, where dimensionality of the values Y (t) is equal to 1. Please note that all of the following derivations can be straightforwardly extended to multiple dimensions of the target variable Y . The reason is that in the case of m > 1 dimensions, the optimization formulation provided below will be equivalent to m independent optimization problems that can be solved separately. 3 by K of 1-by-(d + 1)-dimensional vectors of regression coefficients ̃ θ k = θ 0 k ,θ k , and by a joint diagonal matrix W = diag(w 1 ,...,w D ) with w i ≥ 0 and P D i=1 w i = 1, with a diagonal containing the probabilities that a particular feature dimension i belongs to a relevant subset of features: Y (t) = K X k=1 γ(k,t) θ 0 k + θ k T † k WX(:,t) + δ t ,(1) where δ t is the independent identically distributed scalar-valued random process with expectation zero, and γ(k,t) ≥ 0, P K k=1 γ(k,t) = 1 for all t being the probabilities that Y (t) is generated by the ”local” model k. In another words, Y (t) is predicted as the expectation over K subset-reduced (multiplication with W ) and subspace-reduced (multiplication with T k , where † is the conjugate- transpose operation) (d + 1)-dimensional linear regressions (d << D), where expectation is taken with respect to the (a priori unknown) and time-dependent model probability distributions γ(:,t). Assuming in addition that the rows of feature matrix X are not perfectly correlated with each other, from the Gauss-Markov theorem [31] it follows that the best unbiased estimator of model (1) is provided by the minimum of the following least-squares functional: L ls = T X t=1 Y (t)− K X k=1 γ(k,t) θ 0 k + θ k T † k WX(:,t) ! 2 .(2) Then, using strict convexity of the squared Euclidean norm and deploying the Jensen-inequality, we get the upper bound of (3): L ls ≤L J = T,K X t,k=1 γ(k,t) Y (t)− θ 0 k − θ k T † k WX(:,t) 2 .(3) Please note that L ls ≡L J when γ contains only zeros and ones. Moreover, let μ k ∈R D,1 for all k from 1 to K be the centroid vectors in Euclidean space. Then, following the strategy recently proposed for Entropy-Optimal Manifold Clustering (EOMC) [26], we add to L J (multiplied with a non-negative hyperparameter ε R ) a metrization term γ(k,t)(X(: ,t)− μ k ) † W (X(:,t)− μ k ) 2 , as well as the two regularization terms ε γ φ 1 (γ(:,t)) and ε W φ 2 (W ) 2 This metrization term is two-fold important: (i) as in EOMC, it removes the non-zero kernel from the manifold projec- tion, and (i) it will became indispensable when making predictions for the test data where Y (t) is unknown a priori, since it becomes essential in determining the γ(:,t) -and, hence, which regression should attain which weight in the model (1). 4 with ε γ ,ε W ≥ 0. This results in the Enropy-Optimal Manifold Regression (EOMR) learning of the least-biased version for model (1), given the fixed data X and Y , and, for the fixed selected values of hyperparameters d,K,ε R ,ε γ ,ε W - formulated and implemented as the numerical solution for the following constrained optimization problem: n W,γ ∗ ,μ ∗ 1 ,T ∗ 1 , ̃ θ ∗ 1 ,...,μ ∗ K ,T ∗ K , ̃ θ ∗ K o = arg minL EOMR ,(4) L EOMR = 1 T K,T X k,t=1 γ(k,t) ∥X(:,t)− μ k ∥ 2 W + ε R Y (t)− θ 0 k − θ k T † k WX(:,t) 2 + ε γ φ 1 (γ(:,t)) + ε W φ 2 (W ), s.t. T † k T k = I d ,(5) γ(k,t)≥ 0,and K X k=1 γ(k,t) = 1, ∀t,k,(6) W = diag(w 1 ,...,w D ), w i ≥ 0,and D X i=1 w i = 1, ∀i.(7) In the same way as it was done in EOMC [26] and in the other entropic AI methods, to achieve the entropy-optimal (i.e., least-biased) learning of probability distributions γ, a very good choice for φ 1 (γ(:,t)) would be φ 1 (γ(:,t))≡ log K (γ(:,t)), which, after a multiplication with γ(:,t), will result in a term that maximizes Shannon entropy - and, for all other variables being fixed, L EOMR will have a unique and analytically-computable minimizer for γ, given by a softmax function [9, 26]. A good choice for φ 2 (γ(:,t)) would be to use the Fast Entropy Approximation function φ 2 (W )≡ P D i=1 W i 0.6648 W i +0.2086 − 0.5754W i + 0.0206. This choice has two main reasons: (i) W appears inside of the squared Euclidean norm, resulting in the regularized Quadratic Programming (QP) problem - that, in contrast to the regularized linear problem for γ - does not have an analytic solution for W (at least, no analytic solution to the author’s knowledge); and (i) FEA was recently shown to provide a very close, numerically-cheap, smoothly-differentiable and property-preserving approximation of the Shannon entropy, allowing for orders of magnitude faster and more sparse feature extraction in regression problems, when compared to common sparsification tools like Lasso-regression [32] 5 2.2 Numerical solution of EOMR problem As in the case of EOMC and other entropic AI algorithms, mathematical structure of the optimization problem (4-7) allows deploying subspace-iteration solution, i.e., after select- ing hyperparameter values d,K,ε R ,ε γ ,ε W and initial values of EOMR optimization variables W,γ ∗ ,μ ∗ 1 ,T ∗ 1 , ̃ θ ∗ 1 ,...,μ ∗ K ,T ∗ K , ̃ θ ∗ K , one can iteratively go through all of these EOMR variables, and solve the problem (4-7) for one of these variables at a time, while keeping all other variables frozen. It appears, that for the three of these five variable types (for γ ∗ ,μ ∗ 1 , ̃ θ ∗ 1 ,...,μ ∗ K , ̃ θ ∗ K ) there exist cheap analytical solutions of the problem, that can be computed with the linear memory and time complexity scalings, whereas for the two remaining variable types (for W,T ∗ 1 ,...,T ∗ K ) we propose two linearly scalable numerical algorithms. Step 1: Minimization w.r.t. the manifold projectors T k Freezing all of the EOMR variables except of a T k for a fixed k = 1,...,K, we isolate the terms in L EOMR that depend explicitly on T k . This yields a subproblem: min T k ∈St(d,D) f (T k ) = ε R T X t=1 γ(k,t) (Y (t)− θ 0 k )− θ k T † k WX(:,t) 2 ,(8) where St(d,D) is the so-called Stiefel-manifold [33, 34], defined by the quadratic equality constraint (5). Let z(t) = Y (t)− θ 0 k . Expanding the quadratic term gives: f (T k ) = ε R T X t=1 γ(k,t)z(t) 2 − 2ε R T X t=1 γ(k,t)z(t)θ k T † k WX(:,t) + ε R T X t=1 γ(k,t) θ k T † k WX(:,t) 2 . (9) Using the cyclic property of the matrix trace (tr), the linear interaction term can be rewritten as: T X t=1 γ(k,t)z(t)θ k T † k WX(:,t) = T X t=1 tr γ(k,t)z(t)T † k WX(:,t)θ k = tr T † k " T X t=1 γ(k,t)z(t)WX(:,t)θ k #! . (10) We define the Euclidean gradient of the objective function with respect to T k by differentiating f (T k ) directly: ∇ T k f =−2ε R T X t=1 γ(k,t) Y (t)− θ 0 k − θ k T † k WX(:,t) WX(:,t)θ k .(11) 6 To minimize f (T k ) while adhering to the Stiefel manifold, we project the Euclidean gradient ∇ T k f onto the tangent space T T k St(d,D). The orthogonal projection operator yields the Riemannian gradient gradf (T k ) [33]: gradf (T k ) =∇ T k f −T k T † k (∇ T k f ) + (∇ T k f ) † T k 2 ! .(12) A descent step is taken along the negative Riemannian gradient direction V = − gradf (T k ). To project this tangent vector back onto the manifold surface, we can deploy a QR-decomposition- based retraction mapping R T k (αV), where α > 0 is a step size determined dynamically by an Armijo backtracking line search to ensure sufficient decrease [34]: T k,cand (α) =T k + αV,qrf(T k,cand (α)) =QR =⇒ T k,new =Q· diag(sign(diag(R))).(13) Step 2: Minimization w.r.t. the centroid vectors μ k Isolating the variance-like regularization block that contains μ k , we get: min μ k T X t=1 γ(k,t)(X(:,t)− μ k ) † W (X(:,t)− μ k ).(14) Taking the vector derivative with respect to μ k and setting it to zero we get: ∇ μ k L EOMR =−2 T X t=1 γ(k,t)W (X(:,t)− μ k ) =−2W T X t=1 γ(k,t)X(:,t)− T X t=1 γ(k,t) ! μ ! = 0. (15) Since W = diag(W 11 ,...,W D,D ), this decouples into D independent scalar equations: W i T X t=1 γ(k,t)X i,t − T X t=1 γ(k,t) ! μ i ! = 0, ∀i = 1,...,D.(16) For any feature dimension where W i > 0, dividing by W i yields the exact sample-weighted mean. For uninformative dimensions, where W i = 0, the coordinate has no effect on the functional value. Thus, the global minimizer simplifies to the stable sample-weighted mean vector across all dimensions: μ k = P T t=1 γ(k,t)X(:,t) P T t=1 γ(k,t) .(17) 7 Step 3: Minimization w.r.t. the reduced regression coefficients vectors ̃ θ k = θ 0 k ,θ k When all other variables except of ̃ θ k are fixed, EOMR problem reduces to a linear least-squares problem. Let us define the augmented subspace feature vector ̃x t ∈R (1+d)×1 as: ̃x t = 1 T † k WX(:,t) (18) The regression subproblem becomes: min ̃ θ k T X t=1 γ(k,t) Y (t)− ̃ θ k ̃x t 2 .(19) Taking the partial derivative with respect to ̃ θ k and setting it to zero yields the standard regularized normal equations: T X t=1 γ(k,t) ̃x t ̃x † t ! ̃ θ † k = T X t=1 γ(k,t)Y (t) ̃x t (20) Transposing this linear system gives the explicit global coordinator update for the row vector: ̃ θ k = T X t=1 γ(k,t)Y (t) ̃x † t ! T X t=1 γ(k,t) ̃x t ̃x † t ! −1 (21) Step 4: Minimization w.r.t. the probability distributions γ Applying the Euler-Lagrange principle with respect to γ in (4-7), and following the same way as in the EOMC and other entropic AI methods, one obtains a following unique analytical solution of (4-7) : γ ∗ (k,t) = exp −ε −1 γ (g(k,t)− g(j ∗ (t),t)) P K k=1 exp −ε −1 γ (g(k,t)− g(j ∗ (t),t)) ,(22) where g(k,t) =∥X(:,t)− μ k ∥ 2 W + ε R Y (t)− θ 0 k − θ k T † k WX(:,t) 2 and j ∗ (t) = argmin k g(k,t). Step 5: Minimization w.r.t. the subset feature probability matrix W Since φ 2 (W ) in (4-7) was selected as a Fast Entropy Approximation (FEA), we can directly use the efficient and linearly-scalable Sequential Projected Gradient algorithm proposed in [32] to find the optimizer with respect to W . 8 Total iteration complexity and memory scaling of the subspace algorithm Summarizing the leading order scalings for the five steps described above, we obtain O(KDT + KDdT + Kd 2 T + Kd 3 ) for iteration run time complexity, and O(KDT + KDd + KdT + Kd 2 ) as the memory complexity. Hence, if d << D, both memory and iteration costs scale in the same order as the scaling of computationally very cheap clustering methods like K-means. Steps 1 to 5 are repeated iteratively, and it is straightforward to validate that this procedure would result in the monotonic decrease of the functional value L EOMR , until, at some iteration (I), decrease of L EOMR is less then some predefined tolerance threshold tol. Total run time complexity in the leading order is then O(KDTI + KDdI + KdTI + Kd 2 I) 2.3 Selection of EOMR hyper-parameters d,K,ε R ,ε γ ,ε W Since (4-7) is a supervised learning problem, we can directly apply most of the standard hyper- parameter selection routines from ML and AI, like cross-validation and Bayesian hyper-parameter tuning [35]. Hereby, one first defines some reasonable ranges for the hyperparameters, and then subdivides the data X and Y into training, validation and testing subsets. For each of the particular values of the hyperpaparameters - chosen from the predefined ranges - one iteratively minimizes (4-7) as described in the previous chapter, until the predefined tolerance threshold tol is reached. Then, one evaluates the performance of the identified model on the validation data that was not used in training - and selects the hyperparameter combination with a best validation performance. Finally, the performance of the model with the best validation performance is further evaluated on the test data, that was not used neither in training nor in validation. In the following examples we will use this test performance to compare EOMR with other ML and AI tools, using the same cross-validation train/validate/test data splits. 2.4 Relation to Existing Algorithms EOMR derived here is closely linked to several dimensionality reduction and regression methods in the machine learning literature - that can be considered as special asymptotic cases of EOMR: 9 1. Reduced-Rank Regression (R) and Supervised PCA (SPCA): When K = 1, ε R → ∞, ε W → ∞, and ε γ = 0, the problem approaches classical Reduced-Rank Regres- sion (which fits a multivariate linear regression model under rank constraints [16]), and the standard SPCA, that selects features and projects data based on maximizing a dependency metric (such as Hilbert-Schmidt Independence Criterion) between the data X and the targets Y [27, 30]. 2. entropy-optimal Sparse Probabilistic Approximation (eSPA+): When ε R = 0, and ε γ = 0, problem (4-7) becomes equivalent to eSPA+ [8, 9]. 3 Application examples Below, we present application examples, and compare EOMR results to such AI tools as: (i) to the gradient-boosted random forests and XGBoost (abbreviated as RF/XGB, provided within a functionality of the Matlab function fitrensemble(), with the automated adaptive hyperparameter tuning); (i) to the Lasso-regularized sparse multilinear regression (abbreviated as lasso MLR, provided within a functionality of the Matlab function lasso(), with the automated hyperparameter tuning in a range of 10 −10 , 10 2 for the regularization parameter); to the Deep Neuronal Networks (abbreviated as N, provided within the Matlab Deep Learning Toolbox, with various network architectures, and with the numbers of hidden neurons ranging from 1 to 200, and the number of hidden layers ranging from 1 to 3); and transformer-based TabPFN v3.0 model - a hyperparameter free LLM, specially designed for small data learning problems [11]. For each of these models, we followed the same procedure as for EOMR (described in the Sec. 2.3), and used the same cross- validation data splits to extract the best performing representatives from each model class, to be further compared to each other and to EOMR on the test data - that was held-out of training and validation steps. We also added comparisons to the so-called ”persistent” prediction, when the next value of Y (t) is taken to be the same as its previous value. 10 Predicting leading EOF of Lorenz-96 in chaotic (F=8) and strongly-chaotic (F=12) regimes predictive skill (F=8)predictive skill (F=12) model complexity (F=8)model complexity (F=12) A.B. C.D. lasso MLR parameter sparsity (F=8)lasso MLR parameter sparsity (F=8) EOMR parameter sparsity (F=12)EOMR parameter sparsity (F=8) E.F. G.H. Fig. 1 Analysis results for the L96 in the two chaotic regimes (F = 8, left column, andF = 12, right column). Explanations can be found in the Sec. 3.1. 11 3.1 Analysis of data from Lorenz-63 in strongly-chaotic and very strongly-chaotic regimes Model description The Lorenz-96 (L96) model was introduced by Edward Lorenz in a 1996 paper (published later in 2005) as a ”toy model” of the Earth’s atmosphere for studying the fundamental issues of pre- dictability in chaotic dynamics in spatially extended systems [36, 37]. It mimics aspects of the mid-latitude atmosphere’s non-linear dynamics, such as advection, dissipation, and external forc- ing, within a computationally cheap, periodic one-dimensional domain (a latitude circle). The L96 model is widely used today as a benchmark problem for data assimilation techniques, ensemble forecasting methods, and studies on the general nature of chaos [38–41]. The L96 Type 1 represents a finite difference approximation of a partial differential equation describing a simplified 1D turbulence. Model consists of a system of N coupled ordinary differential equations (ODEs), describing the time evolution of a single scalar atmospheric quantity x j at N equally spaced grid points around a latitude circle: dx j dt = (x j+1 − x j−2 )x j−1 − x j + Ffor j = 1,...,N(23) Periodic boundary conditions are assumed, such that indices are taken modulo N (i.e., x j+N = x j and x j−N = x j ). Variables and terms in (23) have the following meaning: • x j : The value of the atmospheric quantity (e.g., temperature, vorticity) at the j-th grid point. • N : The total number of grid points in the system (system size). Common values in literature are N = 40. • t: Time. • F : A positive, constant external forcing parameter that drives the system. • (x j+1 − x j−2 )x j−1 : The non-linear advection term, which conserves energy in the absence of forcing and damping. • −x j : A linear damping (dissipation) term. Behaviour of the L96 model changes significantly with the forcing parameter F . For small values of F (e.g., F < 1), the system exhibits periodic or steady-state dynamics. As F increases, 12 the system undergoes bifurcations and transitions into chaotic regimes. A commonly studied value is F = 8, which produces robust chaotic behaviour used frequently as a standard benchmark in predictability studies. For regimes where F ≥ 7 (which includes F = 8 and F = 12 investigated below), the system is considered to be in a strong or fully turbulent chaotic state [37]. According to Andrew J. Majda (Courant Institute, deceased in 2021), development of scalable methods capable of robust predictions of L96 in these chaotic regimes represents a ”800-pound gorilla” in the area of chaotic systems. Application of EOMR to L96 output data In the following, we will use the common literature setting for N = 40, and generate long time series x ∈ R 40×2000 of L96 for the two forcing regimes F = 8, 12, with T = 10 (corresponding to 50 Earth days after rescaling of L96-units) and time step τ = 0.005 (corresponding to 36 Earth minutes after rescaling of L96-untis) . Then, we perform the PCA transformation of x 3 , and choose as target Y (t) the values of the dominant PCA mode (i.e., a mode with the largest variance), and choose as X(t) the lag delayed embedding of the ten dominant PCA modes of x at times (t− lag− nτ ), where variable lag is varying between 0.005 and 0.04, and time delay index n ranges between 0 and 9. This results in a feature space dimension D = 100 (see Fig. 1). In each of the forcing regimes we deploy the cross-validation procedure described above, and train EOMR models with hyperparameters sampled randomly from broad ranges of values: d from 1, 2,..., 7, K from ∈ 1, 2, 3, ε R from the interval [1, 100], ε γ from the interval 10 −8 , 10 −2 , ε W from the interval 10 −8 , 10 2 . As can be seen from the Figs. 1A-1B, for all of the considered prediction lag times, both Deep Learning N and RF/XGB perform worse then the trivial ”persistent” predictor, when the next value of Y (t) is taken to be the same as its previous value (solid black lines in Fig. 1). And, this is despite of the quite extensive hyperparameter tuning for N and RF/XGB. Both TabPFN and lasso MLR perform close to each other, and around 7 times better than the ”persistent” predictor. EOMR is the winner in prediction skill, achieving root mean squared error (RMSE) of around 3 In fluid mechanics and geosciences, PCA modes are called Essential Orthogonal Functions (EOFs), that is why we use this abbreviation in the Figures. 13 7× 10 −6 for a lag time 0.005 in L96 time units (36 Earth minutes after rescaling), with around 300 times smaller RMSE when compared to the next competitor lasso MLR, and approximately 2’100 smaller then the error of the trivial ”persistent” predictor. Surprisingly, both RMSE and complexity (Fig. 1C and 1D) of the EOMR models do not change noticeably, when increasing F from 8 to 12 and going from strongly-chaotic to very-strongly chaotic L96 regime. Complexities - i.e., the total numbers of tuneable model parameters that have to learned from the data - is close to the size of the training data - indicating that both of these models really struggle to learn, and just memorize the training data instead. In contrast, EOMR learns very sparse models (Figs. 1G-1H) - even more sparse than the lasso MLR (Figs. 1E-1F, one of the most widely-used sparsification tools) - with just around a hundred of tuneable parameters that should be learned from data. When making predictions for the new data points X(:,t) with the piecewise-linear regression model (1), it is sufficient to keep K of D-dimensional vectors β k = θ k T † k W . Then, making a prediction would mean computing K scalar products β k X(:,t), and, since it requires only very cheap elementary operations like multiplication and addition - and none of the much more expensive operations like divisions - this can be done extremely efficiently, even on the common hardware architectures. For example, evaluating one single EOMR prediction in this L96 example required only around 10 −8 seconds on the commodity Mac LapTop. In contrast, TabPFN v3.0 required around 10 −2 seconds for the same data instances - due to the quadratic (both in D and in T ) scaling of transformer-based LLMs [11], and their huge size. 3.2 Analysis of data from modified Hasegawa-Wakatani (mHW) model of drift-wave turbulence in the edge of a tokamak plasma. One-dimensional L96-model is considered by many to provide only a very limited, and simplified description of turbulent atmospheric dynamics - and not a model of a ”real” turbulence behaviour. To check if the findings from previous Section are induced by this oversimplification - or if they are also reflecting intrinsic properties of ”real” turbulent systems, next we will consider an application of EOMC to the output of a more realistic model from Magnetohydrodynamics (MHD). We will take the Hasegawa-Wakatani model - a seminal two-field fluid description model of drift-wave 14 turbulence in the edge of a tokamak plasma [42–44]. It reduces the complex MHD equations to the evolution of density fluctuations (n) and electrostatic potential (φ) fields. The modified Hasegawa- Wakatani (mHW) model describes the evolution of electrostatic potential φ and density fluctuations n in a two-dimensional slab geometry. In the following application example, we will deploy the mHW model version introduced by Numata et al. [45], involving the resistive coupling term acting only on non-zonal fluctuations. model complexity EOMR parameter sparsity Predicting leading PCA mode of Hasegawa-Wakatani model predictive skill predictions on test data (not used in training) A. B. C. D. E. lasso MLR parameter sparsity Causality and weak stationarity of EOMR autoregression F. Fig. 2 Analysis results for the modified Hasegawa-Wakatani model. Explanations can be found in the Sec. 3.2. Governing Equations Let the zonal average of a field f be defined as ⟨f⟩ = 1 L y R f dy, and the non-zonal fluctuation as ̃ f = f −⟨f⟩. The mHW equations are: ∂ζ ∂t +φ,ζ = α( ̃ φ− ̃n)− μ∇ 4 ζ,(24) ∂n ∂t +φ,n + κ ∂φ ∂y = α( ̃ φ− ̃n)− μ∇ 4 n,(25) 15 where ζ =∇ 2 φ is the ion vorticity anda,b = ∂ x a∂ y b−∂ y a∂ x b is the Poisson bracket representing E × B advection. Variables and parameters of mHW model equations [45] have the following meaning: • φ: electrostatic potential; • n: electron density fluctuations; • ζ: ion vorticity (∇ 2 φ); • α: adiabaticity parameter (resistive coupling); • κ: background density gradient scale length (−∂ x lnn 0 ); • μ: dissipation/viscosity coefficient. Data generation To generate the time series for analysis and comparison, we use the MATLAB code by Jean- Christophe Nave and Denis St-Onge, available at https://github.com/DenSto/HWEsolver. We use the same settings for all of the mWH model parameters as in this code, with the only change being a slightly reduced grid size (64×64 grid points). We generate a time series on the interval [600, 3000] (skipping the outputs between t = 0 and t = 600, when the model ”burns-in” and the transient states die-out), with constant time intervals δt = 0.3125 (which, after rescaling from internal model units to seconds results in a time interval of 3.16· 10 −7 seconds, or 0.316 microseconds) 4 . We follow the same protocol as in the previous example, and perform the PCA/EOF transfor- mation of the obtained simulation time series of ion vorticities and electron density fluctuations. We choose as target Y (t) the values of the dominant PCA mode (i.e., a mode with the largest variance), and choose as X(t) the lag delayed embedding of the three dominant PCA modes at times (t− nτ ), where the time delay index n ranges between 1 and 40. This results in a feature space dimension D = 120 (see Fig. 2). Application of EOMR to mHW output data In each of the forcing regimes we deploy the cross-validation procedure described above, and train EOMR models with hyperparameters sampled randomly from the following ranges of values: d 4 Light in vacuum would travel 94.73 meters during this time interval. 16 from 1, 2, 3, K from ∈ 1, 2, 3, ε R from the interval [1, 10], ε γ from the interval 10 −8 , 10 −3 , ε W from the interval 10 −10 , 10 −5 . Results of the analysis are summarized in the Fig. 2. As in the previous L96-example, despite of the low RMSE on the training data (left panel of Fig. 2A), both Deep Learning N and RF/XGB perform worse then the trivial ”persistent” predictor on the test data (right panel of Fig. 2A). Also TabPFN behaves similarly: being skilful on the training data, it performs worse than the ”persis- tent” predictor on the test data. Such a big discrepancy between the train and test performances is an indication of the overfitting phenomenon, and is also confirmed when comparing resulting model complexities to the training data size (Fig. 2C). Deep Learning N, RF/XGB and TabPFN v3.0 have around the same amount or more tuneable model parameters than the train data size, meaning that these models were ”memorizing” and not learning. Among the common tools only lasso MLR performs with RMSE of 1.5× 10 −3 , around 2 times smaller than the RMSE of the ”persistent” predictor (with RMSE of 4× 10 −3 ), and identifies a sparse model with around 122 tuneable regression parameters (see Fig. 2D), having RMSE performances on train and test data being close together. Optimal EOMR model determined from the cross-validation procedure from Sec. 2.3 appears to have d = 1, K = 1, ε R = 5, ε W = 8· 10 −9 , ε γ = 0. It results in the RMSE of around 4.5· 10 −4 , almost an order of magnitude better than the ”persistent” predictor. EOMR performances on the training and the validation data are very close together (Fig. 2A), and EOMR results in more than an order of magnitude sparser model then lasso MLR, with only around 8 tuneable parameters (compare Fig. 2D and Fig. 2E). Moreover, this very simple entropy-optimal model description for predictions of the dominant PCA mode of HsW that is ”distilled” by EOMR, requires only the lag delayed values of the same dominant PCA mode from up to lag depth 9 - and all other features from the other PCA modes are redundant (see Fig. 2E). In another words, EOMR has found a linear AutoRegressive Moving Average model (ARMA) to be the entropy-optimal regression approximation of the given time series data. This gives us an opportunity to deploy the very well established mathematical theory of ARMA time series processes, for example, the Theorem 3.1.1 on page 83 from [46]. According to this Theorem, ARMA model is weakly-stationary and causal 17 if the roots of the respective characteristic polynomial are strictly outside of the unit circle in a complex plane. As can be seen from the Fig. 2F, all of the roots are strictly outside the unit circle - since all of their absolute values a strictly-larger than 1.0 (the smallest absolute value of a root is 1.015). Hence, the model identified by EOMR is weakly-stationary and causal according to the Theorem, which means that this process will also be asymptotically-stable - and can be used as a stand-alone ”Ersatzmodel” for predictions and long-term simulations. Computing predictions for this EOMR model requires only 15 elementary operations (8 addi- tions and 7 multiplications, no divisions), and takes around 2· 10 −9 seconds - or 2 nanoseconds - on a common Mac LapTop. This is 160 times shorter time than the actual physical lag of 0.316 microseconds for the actual physical plasma process that we aim to predict here. In contrast, com- puting a single prediction instance with LLM-like TabPFN v3.0 model requires around 2· 10 −2 seconds on the same LapTop, i.e., it is around five orders of magnitude slower then the actual lag of the predicted process. 4 Discussion Motivated by the universal approximation theorems, state-of-the-art AI currently follows a track guided by the so-called ”double-descent principle”, which states that as a machine learning model’s capacity increases, its prediction error initially follows a classic U-shaped curve, spikes at the ”interpolation threshold” where it perfectly fits (or ”memorizes”) the training data, and then decreases a second time as it enters a highly overparametrized regime [47]. The downside of this evolution are data-hungry and extremely-large LLM models that tend to memorize - and not to think and to learn [48], with up to tens of trillions of tuneable parameters, and energy- as well as a money-consumption of a whole developed country. In the two challenging examples from chaotic dynamical systems and fluid mechanics provided above, we illustrated the downsides of this evolution: even setting aside the energy consumption and complexity issues, such models (i) consistently exhibited overfitting behaviour (when perfect performance on the training and validation data is followed by a very poor test data performance, being inferior even to such 18 dummy-models as the ”persistent” predictor); and (i) that the very high complexity of state-of- the-art models - and, particularly, the quadratic scaling of LLMs wrt. dimension D and sample size T - lead to extremely long computation times. For example, TabPFN transformer-based model required around 0.02 seconds for predicting a single instance of the dominant PCA mode from the modified Hasegawa-Wakatani model, for a physical prediction lead time of 0.316 microseconds, i.e., around 5-6 orders of magnitude slower than what would have been required for the online prediction of this system. Entropic AI tools [6–10, 12, 26, 32] follow a paradigm that is orthogonal to the mainstream, and aim at finding models that are simultaneously as good as possible (in terms of performance), as well as are as simple and as unbiased as possible (in terms of information theory and Shannon entropy). In this manuscript, EOMR was proposed as a supervised extension of the recently-developed unsupervised EOMC, learning functional relations between features and targets by means of the low-dimensional piecewise-linear regression splines (1). It was shown that the proposed EOMR approach has three major advantages: (i) it allows a simultaneous joint learning of low-dimensional subsets and subspaces of relevant features; (i) it was shown that the iteration computational and memory costs of EOMR scale linearly in D and T ; and (i) formulated in the probabilistic way with entropic regularizations, it aims at finding sparse entropy-optimal models (1). Surprisingly - or, may be, not surprisingly - dependent on the personal background and per- spective, EOMR was shown to be a real game-changer in both of the provided examples. It resulted in the very sparse models (see Figs. 1 and 2), that are simultaneously orders of magnitude more performant and more simple than the considered state-of-the-art AI/ML models. Evaluation of predictions for EOMR model in the Hasegawa-Wakatani example required only 15 elementary oper- ations (8 additions and 7 multiplications, no divisions), and took only around 2· 10 −9 seconds - or 2 nanoseconds - on a common Mac LapTop. This is 160 times shorter time than the actual phys- ical lag of 0.316 microseconds for the underlying physical plasma process to be predicted - and, this performance can for sure be boosted even further on a more specialized hardware. This would open new possibilities for better online predictions, optimization and control of such complex, extremely-fast and chaotic systems. 19 Availability of code. Code and data can be uploaded from https://seafile.rlp.net/u/d/ 58e2937002824805974e/. Acknowledgement. The author would like to thank Davide Bassetti and Tim Prokosch (both RPTU Kaiserslautern-Landau), Rupert Klein (FU Berlin), Lukas Pospisil (VTU Ostrava), Michael Groom and Terry O’Kane (both CSIRO), and the other members of entropic AI community, discussions with whom provided a lot of motivation for this work. This work was funded by the EU Horizons project AI4LUNGS (Grant Agreement No. 101080756). References [1] Ham, Y.-G., Kim, J.-H., Luo, J.-J.: Deep learning for multi-year ENSO forecasts. Nature 573(7775), 568–572 (2019) [2] Ramadevi, B., Bingi, K.: Chaotic time series forecasting approaches using machine learning techniques: A review. Symmetry 14(5), 955 (2022) [3] Ghadami, A., Epureanu, B.I.: Data-driven prediction in dynamical systems: recent develop- ments. Philosophical Transactions of the Royal Society A 380(2229), 20210213 (2022) [4] Jiao, L., Song, X., You, C., Liu, X., Li, L., Chen, P., Tang, X., Feng, Z., Liu, F., Guo, Y., et al.: Ai meets physics: a comprehensive survey. Artificial Intelligence Review 57(9), 256 (2024) [5] Brenowitz, N.D., Cohen, Y., Pathak, J., Mahesh, A., Bonev, B., Kurth, T., Durran, D.R., Harrington, P., Pritchard, M.S.: A practical probabilistic benchmark for ai weather models. Geophysical Research Letters 52(7), 2024–113656 (2025) [6] Groom, M., Bassetti, D., Horenko, I., O’Kane, T.J.: Entropic learning enables skilful forecasts of enso phase at up to 2 years lead time. Journal of Advances in Model- ing Earth Systems 18(1), 2025–005128 (2026) https://doi.org/10.1029/2025MS005128 https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2025MS005128.e2025MS005128 2025MS005128 20 [7] Groom, M., Bassetti, D., Horenko, I., O’Kane, T.J.: Distillation and Interpretability of Ensem- ble Forecasts of ENSO Phase using Entropic Learning (2026). https://arxiv.org/abs/2602. 16857 [8] Horenko, I.: On a scalable entropic breaching of the overfitting barrier for small data problems in machine learning. Neural Computation 32(8), 1563–1579 (2020) https://doi.org/10.1162/ necoa01296 [9] Vecchi, E., Posp ́ıˇsil, L., Albrecht, S., O’Kane, T.J., Horenko, I.: eSPA+: Scalable Entropy- Optimal Machine Learning Classification for Small Data Problems. Neural Computation 34(5), 1220–1255 (2022) https://doi.org/10.1162/necoa01490 [10] Horenko, I., Vecchi, E., Kardoˇs, J., W ̈achter, A., Schenk, O., O’Kane, T.J., Gagliardini, P., Gerber, S.: On cheap entropy-sparsified regression learning. Proceedings of the National Academy of Sciences 120(1), 2214972120 (2023) https://doi.org/10.1073/pnas.2214972120 https://w.pnas.org/doi/pdf/10.1073/pnas.2214972120 [11] Hollmann, N., M ̈uller, S., Purucker, L., Krishnakumar, A., K ̈orfer, M., Hoo, S.B., Schirrmeis- ter, R.T., Hutter, F.: Accurate predictions on small data with a tabular foundation model. Nature 637(8045), 319–326 (2025) https://doi.org/10.1038/s41586-024-08328-6 [12] Bassetti, D., Posp ́ıˇsil, L., Groom, M., O’Kane, T.J., Horenko, I.: An entropy-optimal path to humble ai. arXiv preprint arXiv:2506.17940 (2025) [13] Tibshirani, R.: Regression shrinkage and selection via the Lasso. Journal of the Royal Statistical Society. Series B (Methodological) 58, 267–228 (1996) [14] Tenenbaum, J.B., De Silva, V., Langford, J.C.: A global geometric framework for nonlinear dimensionality reduction. Science 290(5500), 2319–2323 (2000) [15] Jolliffe, I.T.: Principal Component Analysis, 2nd edn. Springer Series in Statistics. Springer, New York, NY (2002). https://doi.org/10.1007/b98835 21 [16] Izenman, A.J.: Reduced-rank regression for the multivariate linear model. Journal of Multi- variate Analysis 5(2), 248–264 (1975) [17] Horenko, I., Schmidt-Ehrenberg, J., Sch ̈utte, C.: Set-oriented dimension reduction: Local- izing principal component analysis via hidden markov models. In: R. Berthold, M., Glen, R.C., Fischer, I. (eds.) Computational Life Sciences I, p. 74–85. Springer, Berlin, Heidelberg (2006) [18] Horenko, I., Klein, R., Dolaptchiev, S., Sch ̈utte, C.: Automated generation of reduced stochas- tic weather models i: Simultaneous dimension and model reduction for time series analysis. Multiscale Modeling & Simulation 6(4), 1125–1145 (2008) https://doi.org/10.1137/060670535 https://doi.org/10.1137/060670535 [19] Maaten, L., Hinton, G.: Visualizing Data using t-SNE. Journal of Machine Learning Research (JMLR) 9, 2579–2605 (2008) [20] Giannakis, D., Majda, A.J.: Nonlinear laplacian spectral analysis for time series with intermittency and low-frequency variability. Proceedings of the National Academy ofSciences 109(7),2222–2227(2012)https://doi.org/10.1073/pnas.1118984109 https://w.pnas.org/doi/pdf/10.1073/pnas.1118984109 [21] Maaten, L.: Accelerating t-SNE using tree-based algorithms. Journal of Machine Learning Research 15(1), 3221–3245 (2014) [22] McInnes, L., Healy, J., Saul, N., Großberger, L.: UMAP: Uniform Manifold Approximation and Projection. Journal of Open Source Software 3(29), 861 (2018) https://doi.org/10.21105/ joss.00861 [23] Healy, J., McInnes, L.: Uniform manifold approximation and projection. Nature Reviews Methods Primers 4(1), 82 (2024) [24] Lotfollahi, M., Wolf, F.A., Theis, F.J.: scgen predicts single-cell perturbation responses. Nature methods 16(8), 715–721 (2019) 22 [25] Peng, D., Gui, Z., Wei, W., Li, F., Gui, J., Wu, H., Gong, J.: Sampling-enabled scalable manifold learning unveils the discriminative cluster structure of high-dimensional data. Nature Machine Intelligence, 1–16 (2025) [26] Horenko, I.: Linearly-scalable and entropy-optimal learning of nonstationary and nonlinear manifolds (2026). https://arxiv.org/abs/2512.17926 [27] Bair, E., Hastie, T., Paul, D., Tibshirani, R.: Prediction by supervised principal components. Journal of the American Statistical Association 101(473), 119–137 (2006) [28] Zou, H., Hastie, T., Tibshirani, R.: Sparse principal component analysis. Journal of Compu- tational and Graphical Statistics 15(2), 265–286 (2006) [29] Hastie, T., Tibshirani, R., Friedman, J.: The Elements of Statistical Learning: Data Mining, Inference, and Prediction, 2nd edn. Springer, ??? (2009) [30] Barshan, E., Ghodsi, A., Azimifar, Z., Zolghadri Jahromi, M.: Supervised principal component analysis: A carbonate-based approach. Pattern Recognition 44(7), 1357–1371 (2011) [31] Plackett, R.L.: A historical note on the method of least squares. Biometrika 36(3/4), 458–460 (1949). Accessed 2026-07-28 [32] Horenko, I., Bassetti, D., Pospisil, L.: Fast, close, non-singular and property-preserving approximations of entropic measures (2026). https://arxiv.org/abs/2505.14234 [33] Edelman, A., Arias, T.A., Smith, S.T.: The geometry of algorithms with orthogonality constraints. SIAM Journal on Matrix Analysis and Applications 20(2), 303–353 (1998) [34] Absil, P.-A., Mahony, R., Sepulchre, R.: Optimization Algorithms on Matrix Manifolds. Princeton University Press, ??? (2009) [35] Wu, J., Chen, X.-Y., Zhang, H., Xiong, L.-D., Lei, H., Deng, S.-H.: Hyperparameter opti- mization for machine learning models based on bayesian optimization. Journal of Electronic Science and Technology 17(1), 26–40 (2019) 23 [36] Lorenz, E.N.: Predictability: a problem partly solved, 1–18 (1996) [37] Lorenz, E.N.: Designing chaotic models. Journal of the Atmospheric Sciences 62(5), 1574–1587 (2005) https://doi.org/10.1175/JAS3430.1 [38] Lorenz, E.N., Emanuel, K.A.: Optimal sites for supplementary weather observations: Exper- iments with a small model. Journal of the Atmospheric Sciences 55(3), 399–414 (1998) https://doi.org/10.1175/1520-0469(1998)055⟨0399:OSFSWO⟩2.0.CO;2 [39] Anderson, J.L.: An ensemble adjustment kalman filter for data assimilation. Monthly Weather Review 129(12), 2884–2903 (2001) https://doi.org/10.1175/1520-0493(2001) 129⟨2884:AEAKFF⟩2.0.CO;2 [40] Houtekamer, P.L., Mitchell, H.L.: A sequential ensemble kalman filter for atmospheric data assimilation. Monthly Weather Review 133(5), 1238–1250 (2005) https://doi.org/10.1175/ MWR2955.1 [41] Bocquet, M., Brajard, J., Carrassi, A., Bertino, L.: Bayesian inference of chaotic dynamics by merging data assimilation, machine learning and expectation-maximization. Foundations of Data Science 2(1), 55–80 (2020) https://doi.org/10.3934/fods.2020004 [42] Hasegawa, A., Wakatani, M.: Plasma edge turbulence. Physical Review Letters 50(9), 682–686 (1983) [43] Wakatani, M., Hasegawa, A.: A collisional drift wave description of plasma edge turbulence. Physics of Fluids 27(3), 611–618 (1984) [44] Gottwald, G.A., Grimshaw, R.: Arakawa-like schemes for the hasegawa-wakatani equations. Journal of Computational Physics 197(1), 210–231 (2004) [45] Numata, R., Ball, R., Dewar, R.L.: A modified hasegawa-wakatani model for plasma turbulence and zonal flows. Physics of Plasmas 14(10), 102312 (2007) [46] Brockwell, P.J., Davis, R.A.: Time Series: Theory and Methods, 2nd edn. Springer Series in 24 Statistics, p. 83. Springer, New York (1991). https://doi.org/10.1007/978-1-4419-0320-4 [47] Belkin, M., Hsu, D., Ma, S., Mandal, S.: Reconciling modern machine-learning practice and the classical bias–variance trade-off. Proceedings of the National Academy of Sciences 116(32), 15849–15854 (2019) [48] Shojaee, P., Mirzadeh, I., Alizadeh, K., Horton, M., Bengio, S., Farajtabar, M.: The Illusion of Thinking: Understanding the Strengths and Limitations of Reasoning Models via the Lens of Problem Complexity (2025). https://arxiv.org/abs/2506.06941 25