
==== Front
PLoS One
PLoS One
plos
PLOS ONE
1932-6203
Public Library of Science San Francisco, CA USA

PONE-D-23-14436
10.1371/journal.pone.0309661
Research Article
Physical Sciences
Mathematics
Algebra
Polynomials
Computer and Information Sciences
Systems Science
Dynamical Systems
Physical Sciences
Mathematics
Systems Science
Dynamical Systems
Physical Sciences
Mathematics
Probability Theory
Random Variables
Covariance
Computer and Information Sciences
Neural Networks
Recurrent Neural Networks
Biology and Life Sciences
Neuroscience
Neural Networks
Recurrent Neural Networks
Physical Sciences
Mathematics
Optimization
Physical Sciences
Mathematics
Operator Theory
Kernel Functions
Computer and Information Sciences
Systems Science
Nonlinear Dynamics
Physical Sciences
Mathematics
Systems Science
Nonlinear Dynamics
Physical Sciences
Mathematics
Discrete Mathematics
Combinatorics
Permutation
Forward variable selection enables fast and accurate dynamic system identification with Karhunen-Loève decomposed Gaussian processes
Dynamic systems identification with fast Gaussian processes
Hayes Kyle Data curation Investigation Software Validation Visualization 1 2
Fouts Michael W. Data curation Investigation Software Validation Visualization 2
Baheri Ali Investigation Methodology Writing – original draft 2 ¤
https://orcid.org/0000-0002-2198-7276
Mebane David S. Conceptualization Formal analysis Funding acquisition Investigation Methodology Project administration Software Supervision Writing – original draft Writing – review & editing 1 2 *
1 National Energy Technology Laboratory, Morgantown, WV, United States of America
2 Department of Mechanical and Aerospace Engineering, West Virginia University, Morgantown, WV, United States of America
Zhou Yu Editor
Inner Mongolia University, CHINA
Competing Interests: The authors have declared that no competing interests exist.

¤ Current address: Department of Mechanical Engineering, Rochester Institute of Technology, Rochester, NY, United States of America

* E-mail: david.mebane@mail.wvu.edu
2024
20 9 2024
19 9 e030966116 5 2023
16 8 2024
© 2024 Hayes et al
2024
Hayes et al
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

A promising approach for scalable Gaussian processes (GPs) is the Karhunen-Loève (KL) decomposition, in which the GP kernel is represented by a set of basis functions which are the eigenfunctions of the kernel operator. Such decomposed kernels have the potential to be very fast, and do not depend on the selection of a reduced set of inducing points. However KL decompositions lead to high dimensionality, and variable selection thus becomes paramount. This paper reports a new method of forward variable selection, enabled by the ordered nature of the basis functions in the KL expansion of the Bayesian Smoothing Spline ANOVA kernel (BSS-ANOVA), coupled with fast Gibbs sampling in a fully Bayesian approach. It quickly and effectively limits the number of terms, yielding a method with competitive accuracies, training and inference times for tabular datasets of low feature set dimensionality. Theoretical computational complexities are O(NP2) in training and O(P) per point in inference, where N is the number of instances and P the number of expansion terms. The inference speed and accuracy makes the method especially useful for dynamic systems identification, by modeling the dynamics in the tangent space as a static problem, then integrating the learned dynamics using a high-order scheme. The methods are demonstrated on two dynamic datasets: a ‘Susceptible, Infected, Recovered’ (SIR) toy problem, along with the experimental ‘Cascaded Tanks’ benchmark dataset. Comparisons on the static prediction of time derivatives are made with a random forest (RF), a residual neural network (ResNet), and the Orthogonal Additive Kernel (OAK) inducing points scalable GP, while for the timeseries prediction comparisons are made with LSTM and GRU recurrent neural networks (RNNs) along with the SINDy package.

http://dx.doi.org/10.13039/100013165 National Energy Technology Laboratory Subcontract No. P010220883 Task 30 https://orcid.org/0000-0002-2198-7276
Mebane David S. http://dx.doi.org/10.13039/100000001 National Science Foundation 2119688 https://orcid.org/0000-0002-2198-7276
Mebane David S. Award 1: DSM and KH, U.S. Department of Energy / National Energy Technology Laboratory site support contract (Leidos), Subcontract No. P010220883 Task 30. http://netl.doe.gov. The funders approved the decision to publish and the overall study design, data collection and analysis, and the final manuscript. Award 2: DSM and MWF, U.S. National Science Foundation award 2119688. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Disclaimer: This project was funded by the United States Department of Energy, National Energy Technology Laboratory, in part, through a site support contract. Neither the United States Government nor any agency thereof, nor any of their employees, nor the support contractor, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. Data AvailabilityAll relevant data are within the manuscript and its Supporting information files.
Data Availability

All relevant data are within the manuscript and its Supporting information files.
==== Body
pmcIntroduction

Gaussian processes (GPs) are stochastic functions that are engines for nonparametric regression. Initially developed for modeling and interpolation in geographic information systems datasets, applications have multiplied across many fields of data science. A key advantage of the GP is its broad, continuous nonparametric support and the amenability of different GP kernels to precise analysis. They are widely recognized as powerful vehicles for static dataset modeling and interpolation with uncertainty quantification.

A GP is Gaussian in that it is a covariance model linking pairs of points on functional draws. As such a GP is completely described by a mean function (often zero in the prior) and covariance kernel. The most famous and perhaps simplest of the covariance kernels is the squared exponential: κ(x,x′)=ς2exp[-(x-x′)2ξ] (1)

where the sill ς2 and range ξ parameters determine the scale and smoothness of the draws. In a typical implementation modeling a static dataset Z, the statistical model Z=δ(x|ς2,ξ)+ϵ (2)

with δ∼N(0,κ) the GP and ϵ an observation error process, is first used to infer the hyperparameters, after which predictions conditioned on the training dataset can be made. The draws on the squared exponential GP—a limiting case of the Matérn covariance family—are infinitely differentiable.

From a practical standpoint the training of the above GP is O(N3), where N is the number of training data points, requiring a Cholesky decomposition of the full covariance matrix. This limits the use of the GP to moderately-sized datasets, generally of a thousand instances or fewer. An important avenue of research is accelerating the speed of both training and inference, such that the GP’s superior modeling features can be trained on large datasets and deployed as machine learning vehicles within other types of models.

Scalable Gaussian processes with inducing points

Liu, et al. [1] provide a thorough overview of efforts that aim to improve scalability while maintaining prediction accuracy using global kernel approximations derived in some sense from a set of M << N inducing points [2–6]. Generally the goal is to approximate the full-rank kernel matrix with local approximations. Of particular note is a O(N) method that directly estimates the covariance with training and inference times that limits the increase in M for large N developed by Wilson, et al. [7]. Some methods employ Analysis of Variance (ANOVA) decompositions to the full kernel which break out contributions in terms of features and their combinations: κ(x,x′)=∑i=1nκi(xi,xi′)+∑i=1n-1∑j=i+1nκi(xi,xi′)κj(xj,xj′)+⋯ (3)

where n is the feature set dimensionality and the individual terms are not necessarily orthogonal. This presents opportunities for variable selection [8]; of particular note is the recent Orthogonal Additive Kernel (OAK) which orthogonalizes the kernels in (3) in order to minimize overlap between main effects and higher-order interactions [9].

Inducing points acceleration of GPs opens up GP regression to large datasets. However, inference is still O(M2) per point, limiting its usefulness in contexts where inference speed is important, such as those where the GP model will be used in the context of control or optimization applications. This motivates a search for alternative methods that are fast in both training and evaluation.

Karhunen-Loève decomposition and BSS-ANOVA

Another approach to scalability in GPs that is distictive to the inducing points approach is the Karhunen-Loève (KL) expansion, in which the kernel is expressed in terms of a sum over its eigenfunctions: δ(x;β)=∑i=1Pβiϕi(x) (4)

where ϕi(x)=λiui(x)∫κ(x,x′)ui(x′)dx′=λiu(x)βi∼N(0,1)

Such methods have the potential to be fast: O(NP2) in training and O(P) per point for inference, where P is the number of terms in the expansion. However such kernels have not been the subject of much research in machine learning contexts generally. The main issues are tractable calculation of the basis functions {ϕi} and a curse of dimensionality [10].

In 2009 Reich and collaborators [11] introduced the Bayesian Smoothing Spline ANOVA (BSS-ANOVA) kernel, which is subject first to an ANOVA decomposition, followed by a KL decomposition. The core of the BSS-ANOVA kernel is: κ1(x,x′)=B1(x)B1(x′)+B2(x)B2(x′)-124B4(|x-x′|) (5)

where Bk is the kth Bernoulli polynomial, defined by the generating function tetxet-1=∑i=0∞Bi(x)tii! (6)

yielding B1(x)=x-12B2(x)=x2-x+16B4(x)=x4-2x3+x2-130

This kernel is effectively a sum of a non-stationary quadratic response surface—corresponding to the first two terms in (5)—and a stationary deviation (the final term). Covariances for higher-order interactions are constructed with dyadic products of the main effect covariance given in (3): κ2([xj,xk],[xj′,xk′])=κ1(xj,xj′)κ1(xk,xk′) (7)

and so on for higher-order interactions. Terms are then multiplied by scaling hyperparameters and added together to produce the full kernel: κ=σ02τ02+σ12τ12∑i=1nκ1,i+σ22τ22∑i=1n-1∑j=i+1nκ2,ij+⋯ (8)

The kernel so constructed is supported on a second-order Sobolev space [11], a broad support which is one of its primary advantages. The significance of separate σ2 and τ2 scaling parameters will become clear below in the context of Bayesian linear regression.

Building the kernel in this fashion effectively addresses the problem of generating the eigenfunctions from the KL decomposition: because all of the terms in (8) are based on the generative kernel (5), the KL decomposition of (8) will depend only on eigenfunctions of κ1. Additionally if all input features are normalized to an [0, 1] interval (we restrict the discussion to continuous input features for now), then it is only necessary to compute a single set of basis functions {ϕi}. The decomposed BSS-ANOVA GP is written: δ(x;β)=β0+∑i=1n∑k=1∞βikϕk(xi)+∑i=1n-1∑j=i+1n∑k=1∞∑l=1∞βik,jlϕk(xi)ϕl(xj)+⋯ (9)

Given the assumption σ02τ02=σ12τ12=σ22τ22=⋯=σ2τ2 (10)

then the priors for the coefficients β are iid normal β·,·∼N(0,σ2τ2) (11)

where the notation ⋅ indicates an arbitrary index. While any kernel is amenable to a KL decomposition and as such to the variable selection routine described below, the form of the BSS-ANOVA kernel as a sum of scaled terms (which separates the basis decomposition from the hyperparameters) and its broad support make it particularly suitable for many applications, and we restrict our analysis to this kernel for the remainder of the article.

Following [11] we generate the set {ϕi} by producing κ1 for a dense grid consisting of 500 intervals on [0, 1], eigendecompose and fit to cubic splines. Fig 1 shows the first 6 basis functions. These basis functions are nonparametric, pairwise orthogonal, and ordered: note the increase in frequency and decrease in amplitude as the orders increase.

10.1371/journal.pone.0309661.g001 Fig 1 The first six basis functions of the KL-decomposed BSS-ANOVA kernel.

The basis is nonparametric, spectral, pairwise orthogonal and ordered.

Variable selection

It’s clear from (9) that the number of terms in the expansion can increase rapidly, even for low-dimensional input spaces. A key component of applying the GP to a modeling problem is thus the selection of terms. Effectively we seek to minimize the objective function Φ(β)=||Z-δ(x;β)||2+ζ(β) (12)

where ζ is a penalty function which leads to a sufficiently sparse solution.

Indicator variable methods

Reich, et al. [11] took a hierarchical Bayesian approach to the problem, estimating a separate variance τ2 for each term in the expansion, which is in turn expressed in terms of an indicator variable with a Bernoulli prior. This approach, like other ‘indicator variable’ methods, accomplishes the variable selection and the training simultaneously and comprehensively, at the cost of requiring a large number of variables in the prior model and a computationally onerous Markov chain Monte Carlo (MCMC) sampling procedure.

Other sparse optimization methods such as ridge regression or LASSO share the limitation that many high-order terms must be included in the initial model before downselection occurs.

Forward variable selection

The ordered and orthogonal nature of the basis functions suggests a forward variable selection approach. Rewriting the model (9) for a basis function set of maximum order q, δ(x;β)=β0+∑i=1n∑k=1qβikϕk(xi)+∑i=1n-1∑j=i+1n∑k=1q∑l=1qβik,jlϕk(xi)ϕl(xj)+⋯ (13)

then considering a model building procedure which increases q stepwise starting with q = 1 reveals that each subsequent step adds n main effect terms (each depending on a single input—n is the number of inputs), n2[2(q−1)+1] two-way interactions, and n3[3(q−1)2+3(q−1)+1] three-way interactions. As the model order increases the L2 truncation error for the full kernel decreases as (for the case of a single input) [10]: ||κ(x,x′)-∑i=1qϕi(x)ϕi(x′)||<(∑i=q+1∞λi2)1/2 (14)

Since the eigenvalues of the BSS-ANOVA kernel decomposition decrease quickly with increasing order, an approach to the optimization problem (12) focusing on low-order models will sacrifice little in the way of accuracy while realizing significant advantages in both training and inference times.

The design and implementation of such an approach is the main contribution of this work. It approaches the optimization of (12) with an iterative process, finding the most efficient truncation of the system while evaluating the cost function only for candidate models with fewer terms than the optimum truncation. The method is fully Bayesian, with a fast linear Bayes sampling procedure at its core. As such the form of the cost function is also Bayesian in nature, taking the form of the Bayesian or Akaike information criteria (BIC/AIC), which incorporate L0 penalties: BIC=PlnN-2lnL^ (15)

AIC=2P-2lnL^ (16)

where L^ is the maximum of the likelihood function on the parameter space.

The following sections describe the components of the optimizer in detail and present the implemented variable selection algorithm.

Linear Bayesian regression

Given a statistical model with a fixed number of terms zi=δi(xi;β)+ϵ(σ2) (17)

with ϵ a white noise observation error of variance σ2 and δ is a BSS-ANOVA model corresponding to a given truncation to the KL expansion (9), the model is linear in the coefficients β. Using priors that are conjugate to the likelihood for the coefficients and the observation error variance, a Gibbs sampling methodology can quickly obtain the posterior for the parameters.

The conjugate prior for β is zero-mean iid normal: β∼N(0,σ2τ2I) The conjugate prior for σ2 is inverse gamma: σ2 ∼ IG(a, b), with a and b the shape and scale parameters, respectively; likewise the conjugate prior for τ2 is inverse gamma: τ2 ∼ IG(aτ, bτ). The Gibbs sampler functions iteratively, such that for fixed {σ2, τ2}, β∼N(μ,Σ), with μ=(XTX+1/τ2I)-1XTZ (18)

Σ=σ2(XTX+1/τ2I)-1 (19)

where X∈RN×P is a matrix constructed from the basis functions appearing in the expansion. Its rows correspond to instances and columns to terms in the expansion. For fixed {β, τ2}, σ2 ∼ IG(a*, b*), with a*=a+N/2+P/2 (20)

b*=b+12[(μ-β)T(XTX+1/τ2I)(μ-β)+ZTZ-μTXTZ] (21)

For fixed {β, σ2}, τ2∼IG(aτ*,bτ*), with aτ*=aτ+P/2 (22)

bτ*=bτ+12σ2βTβ (23)

For more details see [12].

Optimization via forward search

The algorithm constructs models with terms having up to three-way interactions. The algorithm runs sequentially through stages labeled by an integer index I initialized at 1. At each stage, a set of trial terms come in to the expansion for which the sum over all basis functions appearing in the term add up to I. Since there are many permutations of basis functions that form this sum, and for any given set of function orders in a term, different ways they can be permuted among the inputs, adding terms at each stage occurs in substages. Each substage adds terms corresponding to a particular set of basis functions, and includes all input permutations. For example: stage 1 adds only first order main effects: ϕ1(x1) and ϕ1(x2) for a two-input dataset. Stage 2 adds second order main effects and first order two way interactions—ϕ2(x1), ϕ2(x2) and ϕ1(x1)ϕ1(x2)—in two substages. The substages occur such that terms involving lower-order basis functions come first. For the example in the case of stage 2, this is the first order two-way interactions ϕ1(x1)ϕ1(x2). Each substage adds at once all combinations of inputs and all permutations among each combination, such that each substage adds 2n2 terms for two-way interactions and 6n3 terms for three-way interactions (for the case where all function orders are different; fewer when two or more are the same). With trial terms added at a given substage, the sampler estimates model parameters and calculates the BIC or AIC. The routine terminates at an optimum BIC or AIC. Because there can be local minima, a “tolerance” setting controls how many substages the algorithm can iterate through without finding a new minimum BIC or AIC before it terminates. The terminated algorithm returns the optimum model.

Algorithm 1 BSS-ANOVA forward variable selection algorithm

1: procedure FwdVarSelect(x, z, ϕ, tol)

2:    ⊳ x is a matrix of inputs, columns are features and rows are instances

3:       ⊳ z is a column vector of data

4:      ⊳ ϕ is an ordered set of basis functions

5:       ⊳ tol is an integer tolerance

6:  Form a column vector of ones RN∍X=1

7:  Set I=1

8:  Set count = 0

9:  while count < tol do

10:   Set Z={j∈Z:0≤j≤I}

11:   Set Q={(i∈Z,j∈Z,k∈Z):i+j+k=I}

12:   Order q∈Q s.t. qm ≺ qn when max(qm) < max(qn)

13:   for l={1,2,⋯|Q|} do

14:    Set md = {ql, 0, 0, ⋯} s.t. |md| = |x|

15:    Form Md=S(md)    ⊳ Md contains all permutations of md

16:        ⊳ Each element of Md is a term in the expansion

17:    Build Xd where Xd,ij=∏k:Md,jk≠0ϕMd,jk(xik)

18:    Recursively concatenate: X = [XXd], M=[M;Md]

19:    Call the sampler: β, BIC = gibbs(X, z, ϕ)

20:    if the BIC is a minimum for all models then

21:     count = 0

22:    else

23:     count = count + 1

24:   I=I+1

25:  Return M, β, BIC

Experiments: Dynamic system identification

Procedure

BSS-ANOVA regression—as is the case for other GPs—is most effective for tabular datasets with continuous inputs and targets of moderate dimensionality. This suggests an application in dynamic systems identification. Indeed BSS-ANOVA GPs have been utilized as components of other models (“intrusively”) for this purpose in a number of applications [13–15]. We demonstrate here that they may also be used directly to identify dynamics in more general cases, without the aid of an accompanying model. In treating the dynamic systems application we should be clear to make a distinction between the present approach in which the GP models the static relationship between the time derivatives of states and the states and forcing functions themselves, and other “dynamic GP” applications, e.g. [16].

The procedure is a concurrent one, in that time derivatives estimated from the datasets are modeled directly using BSS-ANOVA with forward variable selection, using the concurrent values of the system states and other inputs; for example a two-state system is modeled using two separate GPs: x˙1=δ1(x1,x2,u) (24)

x˙2=δ2(x1,x2,u) (25)

The identified system is then integrated to yield predictions with uncertainty.

The procedure was demonstrated on two nonlinear dynamic datasets: a synthetic dataset derived from the susceptible, infected, recovered model (SIR model) for infectious disease, and the ‘Cascaded Tanks’ experimental benchmark dataset. In both cases comparisons were made to long short term memory (LSTM) and gated recurrent unit (GRU) neural networks, along with the sparse identification of nonlinear dynamical systems (SINDy) package [17] for timeseries prediction. In the case of the cascaded tanks benchmark comparisons were made against random forest (RF), a residual neural network (ResNet) and the state-of-the-art OAK inducing points scalable GP [9] for the static derivative estimation problem.

Experimental benchmark: Cascaded tanks

The cascaded tanks nonlinear benchmark dataset is an experimental nonlinear dynamic system [18]. The experiment consists of a set of two tanks and a reservoir of water. An upper tank is filled by a pump from the reservoir. An outlet in the upper tank empties into the lower tank, which in turn empties through an outlet back into the reservoir. A signal sent to the pump serves as the forcing function for the system, with the tank water level heights the two states of the system.

Hyperparameters are an important component of any comparison of numerical methods. While sometimes it is advantageous to perform hyperparameter optimization via sweeps or other methods, it is also important to compare methods using reasonable values that are in the neighborhood of commonly used defaults, with trial-and-error adjustments, as this is how most users will approach the task. This is therefore the approach we took in the following comparisons. An exception is for the comparisons with the SINDy package, for which a sweep was performed as described below.

We first compared the performance of BSS-ANOVA with RF, ResNet and OAK static regressors. Derivatives were calculated via direct finite differences for the relatively noise-free dataset, yielding 10000 instances. Each method was trained on concurrent values of both states and the forcing function for each derivative. For the GP we used hyperparameters of a = 1000, b = 1.001, aτ = 4 and bτ = 55 for h˙1 and 69.1 for h˙2, with tolerances of 3 for h˙1 and 5 for h˙2, and the AIC as discriminator. Inputs were normalized to [0, 1] using min-max scaling. Of 2000 draws the first 1000 were discarded. Only two-way interactions were required. For the RF 100 trees were used with a leaf size of 5. The ResNet had a depth of 6 (filter sizes ranging from 16 to 64) and in between each fully connected layer is a batch normalization and relu layer. The mini batch size is 16, initial learn rate is 0.001, the data was shuffled every epoch for a total of 30 epochs, and the validation frequency was 1000. OAK was applied at a maximum dimension of 3 and with the default value of 200 inducing points. The 5-fold cross-validated results appear in Table 1. OAK performed best for both outputs, followed closely by BSS-ANOVA. Both GPs outperformed the RF and the ResNet by clear margins.

10.1371/journal.pone.0309661.t001 Table 1 Cascaded tanks 5-fold cross validated accuracies: Derivatives, mean absolute error (MAE).

Method	h˙1 (MAE/10−4)	h˙2 (MAE/10−4)	
OAK	17±4.7	36±2.4	
BSS-ANOVA	18±6.5	39±3.6	
ResNet	36±14	61±15	
RF	30±9.4	49±4.9	

The nature of the decomposed GP as a sum over terms pertaining explicitly to certain inputs is a distinct advantage of the KL decomposition. The results indicate that the most important factors in the timeseries model for tank 1 are the water levels in both tanks, with tank 2 slightly more important. This is counter-intuitive, since the leve in tank 2 (the lower tank) does not even appear in a naive physical model of the system [18]. The analysis shows that the actual system contains significant feedback. For the dynamics of the level in tank 2, the most important term is an interaction between the tank 2 level (second order basis function) and the pump signal (first order).

Timeseries predictions follow for the GP via a 4th-order Runge-Kutta integration routine. These were compared with LSTM and GRU recurrent neural networks (RNNs), along with the SINDy package. For the LSTM there was one LSTM layer and a total of 128 hidden layers, the data was shuffled every epoch for a maximum of 125 epochs, verbose was equal to 0, and the sequence was padded to the left. The GRU had one GRU layer and 150 total hidden layers, the data was shuffled every epoch for a total of 150 epochs, verbose was equal to zero and the sequence was padded to the left. The 5-fold cross-validated results (datapoints were not randomized before creating the folds so as to preserve the timeseries order) appear in Table 2. BSS-ANOVA is most accurate, followed by the LSTM and the GRU. Fig 2 shows the predictions of the GP and the LSTM for the upper tank for one of the test folds. The GP predictions are superior near the sharp inflection and critical points where nonlinearities are strongest. Note that the first 50 points of each test set, which were provided to the LSTM and GRU as a start-up set in the prediction phase, were removed from the calculation of error for both methods.

10.1371/journal.pone.0309661.g002 Fig 2 (a) BSS-ANOVA and (b) LSTM predictions vs. test set data for the water level height in tank 1 of the cascaded tanks dataset. Shaded regions in (a) are 95% confidence bounds estimated from a draw of 40 curves.

10.1371/journal.pone.0309661.t002 Table 2 Cascaded tanks 5-fold cross validated accuracies: Timeseries, mean absolute error (MAE).

Method	h1 (MAE)	h2 (MAE)	
BSS-ANOVA	0.1167±0.0382	0.1577±0.0334	
SINDy	0.1391±0.0631	0.1768±0.0695	
LSTM	0.2345±0.1006	0.2296±0.0378	
GRU	0.3243±0.1092	0.2481±0.0402	

To evaluate the SINDy performance on the Cascading Tanks task, an exploratory analysis was performed over the hyperparameters (optimizer, threshold, alpha, and basis function libraries). Optimizers were first tested using default alpha (0.05) and basis functions (2nd order polynomials), and a threshold of 0.001 due to the small coefficients of the model terms. Of the 8 optimizers tested, 2 produced errors and were not able to be evaluated. Of the 6 remaining optimizers, all produced similar results. STLSQ, SR3, and Constrained SR3 were marginally the best performing models and were used for basis function evaluation.

1st, 2nd, 3rd, and 5th order polynomials, and 3rd order polynomials with the addition of Fourier functions were tested successfully. The largest improvement was seen moving to 3rd order polynomials with an additional small improvement adding in the Fourier functions. 4th order polynomials and 2nd/4th order polynomials with Fourier functions all failed to converge.

The three chosen optimizers were used with 3rd degree polynomials with Fourier functions in a threshold parameter search to find the best results with a default alpha of 0.05. The optimum occurred at a threshold of 0.0001 using the STLSQ optimizer. Lastly, a search was performed over alpha to determine the optimal value which occurred at 0.05 (default).

While it is reasonable to expect that OAK with 200 inducing points would outperform BSS-ANOVA in the time integration, it was not practical to make this comparison for reasons of computing time. A comparison with a reduced number of inducing points and increased time step in the integrator was made—results are discussed below.

Comparing the results of the sparse KL-decomposed BSS-ANOVA method with other published treatments of the same benchmark is difficult because of the lack of rigorous separation between training and test data in most applications. One study of note is that of La Cava, et al. [19], who used an approach to reconstructing the dynamic system from a predetermined set of parametric terms using a version of genetic programming. Assuming clear separation in this study between test and train sets, the approach described in this paper outperforms the “Epigenetic Linear Genetic Programming” approach by over an order of magnitude in mean square error, while also outperforming the NARX-NN approach (a lightly parameterized recurrent neural network) offered as comparison in [19] by approximately one order of magnitude.

Synthetic benchmark: Susceptible, infected, recovered model

The susceptible, infected, recovered model (SIR model) is a common simulation for infectious disease. Though there are several versions, the simplest is three states, only two of which are independent. The system is written S˙=-BIS/NP (26)

I˙=BIS/NP-γI (27)

R˙=γI (28)

where S(t) is the susceptible population, I(t) the infected, R(t) the recovered, B(t) is the transmissibility (which we utilize as a forcing function), γ is the recovery rate (which we leave fixed at 0.5) and NP is the total population. Because NP is fixed and S + I + R = NP, only two states are independent, so the system dynamics can be captured by modeling only two of the three. We chose I(t) and R(t).

The training data consists of 58 curves. All curves in the training set have a fixed B value ranging from 0.5 to 9, in six intervals of 1.7. For each value of B there are 8–10 siumulations corresponding to different initial conditions designed in such a way to provide coverage of the state space. (Exact initial conditions used appear in the supplement.) Each simulation used NP = 1000.

The test data consists of 24 curves, each of which features a temporally changing transmissibility B(t). There are three initial B0 values: 1.35, 4.75 and 8.15. For each starting point there are two types of transmissibility curves: a ramp and a sinusoid. The B0 = 1.35 and B0 = 4.75 starting points have ramps with a positive slope of 1, while the B0 = 8.15 curves have a slope of -1. All ramps run from t = 0 to t = 4, where they level off. The sinusoids have amplitudes between 0.5 and 3 and a period of 1.

Hyperparameters for BSS-ANOVA were: a = aτ = 4 for both states, bτ,R = 8.95 and bτ,I = 72.1, while bI = 1.25 and bR = 20. 2000 draws were taken and the first 1000 discarded. The tolerance was 6. Hyperparameters for SINDy, LSTM and GRU were the same as for the Cascaded Tanks.

Results for BSS-ANOVA and SINDy are shown in Table 3. SINDy performed better (in average) on both state predictions. Statistics were not calculated for the GRU and LSTM as each failed to replicate the dynamics in most test cases and were obviously inferior to both BSS-ANOVA and SINDy in every instance. A graphical comparison for a selected number of curves from the test set are shown in Fig 3. The largest overall term in the model for the infected population was a two-way interaction between the infected state and recovered state, consisting of second and third order basis functions.

10.1371/journal.pone.0309661.g003 Fig 3 BSS-ANOVA and SINDy predicted results for I in three representative test datasets.

(a) and (b) have sinusoidal and and (c) upward ramp control function dynamics that were not present in the training set. Red circles are test data, the blue curve is the mean and grey shading the 95% confidence bounds for BSS-ANOVA, while the green curve is SINDy with 3rd-order polynomials plus Fourier terms.

10.1371/journal.pone.0309661.t003 Table 3 SIR test set results for BSS-ANOVA and SINDy: Mean absolute error (MAE).

Method	I (MAE)	R (MAE)	
BSS-ANOVA	5.2739±4.0183	11.8345±21.7337	
SINDy	1.8722±4.1736	8.0659±22.2364	

Training and inference times

Training and inference times for BSS-ANOVA were fast, with a mean total train time of 6.3 seconds for the cascaded tanks and 10.8 seconds for the SIR, with 8,000 and 20,000 training data points, respectively, on a 2019 6-core i7 processor with 16 GB of RAM. The routines were implemented in MATLAB, but not parallelized or optimized for speed. Models for h˙1 contain between 23 and 41 terms, while h˙2 has between 38 and 57 terms. Prediction times for 2000 static points for the cascaded tanks averages 0.5437 s, and the time for evaluating integrals over the test set averages 20.22 s. For the SIR model the I˙ model had 81 terms and the R˙ model 9 terms, with a mean integration time of 5.3 s. Analyses have shown that the rate limiting step in BSS-ANOVA build algorithms are the O(NP) construction of the X matrix from the inputs and basis functions. The neural networks were native MATLAB functions, parallelized and optimized for speed. Nonetheless train times were considerably longer, with mean train times of 130s for the ResNet and 175 and 123 s, respectively, for training the LSTM and GRU for the cascaded tanks. This is to be expected given that the number of weights in the neural nets are on the order of 104.

It was not feasible to integrate OAK at the level of 200 inducing points to the same standard as that of BSS-ANOVA because of time considerations. A reduced set of 40 inducing points yielded accuracies in the static estimation problem that were approximately the same as BSS-ANOVA. A reduced time step (500 vs. 20,000 integration steps) brought the integration time down to 51 minutes for OAK, with MAE/MAPE of 0.1554/6.3 for h1 and 0.2378/9.1 for h2. Reducing the integration step to the same level as BSS-ANOVA (where we could expect comparable integration accuracies) would require approximately 33 hours.

Discussion

The results show that the forward variable selection methodology makes the KL-decomposed GP a viable option for dynamic systems—competitive in these preliminary results with state-of-the-art routines in both static and dynamic modeling tasks—due to its combination of speed and accuracy. In addition the Bayesian nature of the method yields estimates of uncertainty in the predictions, which can be useful in design of experiments and optimization. It also creates opportunities for fast Bayesian model updates in a control context, wherein only new data need be taken into account due to the presence of a strong prior probability distribution constructed from previously-utilized data.

The neural networks were not able to obtain the dynamics in the SIR test because the RNNs map the recurrent inputs to outputs directly on the Hilbert space of time-dependent states and control functions. That is, if the specific type of time-dependent control function behavior found in the test set does not appear in the training set—as is the case in the SIR task—then the test set is out-of-distribution for the RNN. SINDy and BSS-ANOVA by contrast construct static models of the system in the (much more tractable) Euclidean state/control function space. Since the test set forcing function never exceeds the bounds of the training set in that Euclidean space, the test set is in-distribution for both methods.

More speculatively, the performance difference between SINDy and BSS-ANOVA on the different tasks may arise from the relative suitability of the basis to the different types of datasets represented. The nonparametric GP basis was better suited to the experimental benchmark (where the dynamics are nonparametric) while the 3rd-order polynomial SINDy basis excelled in the synthetic case that arises from dynamics which are in fact generated from first, second and third-order polynomials. It will be interesting to explore this hypothesis further in future experiments.

There are several ways in which future work might improve on the current algorithm. A more discerning selection of terms at each substage in the variable selection routine—adding some, but not necessarily all terms at each stage—may improve performance by limiting overfitting while also automatically selecting features. Here we have been exploring the use of simulated annealing and Markov Chain Monte Carlo approaches modified to retain the speed advantages of fast Gibbs sampling as potential innermost optimization processes. Such potential improvements along with code optimization and GPU acceleration will be included in future versions of the method.

Conclusion

A new forward variable selection algorithm has made the scalable Gaussian process BSS-ANOVA a fast and accurate method for nonparametric regression of tabular data on continuous input spaces. The speed and accuracy for this type of dataset makes it an advantageous method for dynamic system identification. Favorable comparisons with successful and popular sparse basis and neural network approaches for timeseries problems were made in prediction tasks for a pair of nonlinear synthetic and experimental dynamic datasets.

Supporting information

S1 File Code and data used in the production of this manuscript: MATLAB results for all methods except SINDy.

(ZIP)

S2 File Code and data used in the production of this manuscript: SINDy results for Cascading Tanks.

(ZIP)

S3 File Code and data used in the production of this manuscript: SINDy results for SIR.

(ZIP)

Many thanks to Charlie Harmison and Josh Caswell for producing some of the neural network results.

10.1371/journal.pone.0309661.r001
Decision Letter 0
Landenmark Hanna Staff Editor
© 2024 Hanna Landenmark
2024
Hanna Landenmark
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version0
21 Jul 2023

PONE-D-23-14436Forward variable selection enables fast and accurate dynamic system identification with Karhunen-Loève decomposed Gaussian processesPLOS ONE

Dear Dr. Mebane,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

Please find the reviews from two reviewers below. Please note that both reviewers have expressed concerns that the algorithm and motivations behind its structure and implementation is not clear, and we ask that you thoroughly address this in order to ensure that the manuscript is fully reproducible.

Please submit your revised manuscript by Sep 04 2023 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Hanna Landenmark

Staff Editor

PLOS ONE

Journal requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. Thank you for stating the following in the Acknowledgments Section of your manuscript:

“Many thanks to Charlie Harmison and Josh Caswell for producing some of the neural

network results.

DSM and KH were supported by the U.S. Department of Energy, National Energy

Technology Laboratory through a site support contract: Subcontract No. P010220883

Task 30

DSM and MWF were supported by the National Science Foundation under award

2119688”

We note that you have provided funding information that is currently declared in your Funding Statement. However, funding information should not appear in the Acknowledgments section or other areas of your manuscript. We will only publish funding information present in the Funding Statement section of the online submission form.

Please remove any funding-related text from the manuscript and let us know how you would like to update your Funding Statement. Currently, your Funding Statement reads as follows:

 “Award 1: DSM and KH, U.S. Department of Energy / National Energy Technology Laboratory site support contract (Leidos), Subcontract No. P010220883 Task 30. http://netl.doe.gov. The funders approved the decision to publish and the overall study design, data collection and analysis, and the final manuscript.

Award 2: DSM and MWF, U.S. National Science Foundation award 2119688. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.”

Please include your amended statements within your cover letter; we will change the online submission form on your behalf.

3. We noted in your submission details that a portion of your manuscript may have been presented or published elsewhere. [There is a preprint copy on arXiv: https://arxiv.org/abs/2205.13676] Please clarify whether this publication was peer-reviewed and formally published. If this work was previously peer-reviewed and published, in the cover letter please provide the reason that this work does not constitute dual publication and should be included in the current manuscript.

4. Please remove your figures from within your manuscript file, leaving only the individual TIFF/EPS image files, uploaded separately. These will be automatically included in the reviewers’ PDF.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: I Don't Know

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: No

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: This paper focuses on developing scalable Gaussian process methods with applications to accurate dynamic system identification. To do so, the authors adopt the idea of Karhunen-Loeve decomposition and propose a forward variable selection approach to incrementally select the terms in the decomposition. Techniques in Bayes model selection are used to accomplish the goal. The authors provide numerical experiments to showcase the performance of their method.

I found the general idea natural to understand, but the algorithm described in the section "Optimization via forward search" and "Algorithm 1" are VERY DIFFICULT for me to comprehend. The authors should improve the clarity of the writing since the algorithm is the paper's main contribution. For example, can the authors be more explicit about the different stages of the algorithm and how different terms are added? Moreover, can the authors specify how the BIC and AIC are calculated and used? I found these parts poorly explained in the paper, making readers hard to follow what the authors have exactly done.

Other comments:

1. In the explanation of Karhunen-Loeve decomposition on page 4, I am confused about the notations phi_i and u_i. Is u_i normalized? If so, the covariance of (4) seems not equal to k. Please explain. The same question applies to equation (14). Please check.

2. For equations (16)-(21), I would suggest adding a reference or some details so readers can follow more easily.

3. Please be more explicit about the terminology used. For example, the author never explained the meaning of "MAE" (e.g., in Table 1)

4. Notations and grammar:

(1) Page 6, equation (11), beta_{dot, k}. What is the meaning of "dot"?

(2) Page 5, line 59. "supported by" -> "supported on"?

5. The quality of the figure should be improved. Currently, it is hard to read and compare different lines.

Reviewer #2: The authors propose a new approach for variable selection for Gaussian processes by making use of the Bayesian smooting spline ANOVA kernel (BSS-ANOVA) a (numerically determined) Karhunen-Loeve decomposition.

I have a list of comments and suggestions to the authors that should be considered before re-committing the paper:

1. I think it should be pointed out more clearly in the abstract already what exactly is new in the method. What was there before and what exactly is your innovation?

2. I think it should be added in the abstract how efficient the proposed method is, in terms of complexity w.r.t. size of the data set, dimensionality etc.

3. Concerning Eq. (3): Some authors mean by the ANOVA decomposition already a uniquely defined decomposition whichs terms are orthogonal w.r.t. to an inner product. In your setting, the single parts in (3) do not neccessarily have to fulfill some orthogonality condition?

4. A lot of parameters are introduced and re-used some pages later wihtout naming them or refering to them. Reading the paper for the first one will face difficulties. Some notations are a somewhat mixed up and not so clear, example: in Eq. (2) it is delta(x|zeta^^2,xi), then in (4) it is delta(x;beta) ... is it the same?

5. Equation (5): what is the motivation to look at such kernels? What are special properties or advantages in comparison to other frequently applied kernels. Could or should it be used in other frameworks, too?

6. Equation (8): It is not so clear to my, why two types of parameters (taus and sigmas) are needed, they are just multiplied.

7. Equation (12): This function will be minimized with respect to the betas. Which regularization are you coing to be using? Which optimization algorithms and implementations do you use?

8. Equation (14): What are the lambdas?

9. Equations (16) and (17): Using the eigenfunctions and constructing the matrix X, you could switch from the kernel formulation to the primal formulation of a Gaussian process. Is this what is expressed by (16) and (17)? Equation (16) is just like kernel ridge regression or linear regression with non-linear feature map.

10. The Section "Variable Selection" is not clearly structured. What is the relation between the different subsections and sub-subsections (not nicely visible)? What is the final result? Please, refer to the Algorithm below. Complexity is not adressed.

11. In general, I was wondering why sections and subsections are not numbered.

12. ANOVA means analysis of variance. I think, the found betas should give some interpretability (which dimensions are important and which are coupled)!? What insights do you gain from that in the considered numerical examples?

13. Numerical Experiments / Data sets: Please comment a bit on data pre-processing. Your x's have to be scaled to [0,1]. How was that achieved? Already given or via min-max scaling? What is the dimension of the data in the examples? The used hyper parameters have been tuned automatically or tuned/set by hand? What is the outcome of the variable selection procedure? How do the betas look like and which information can you gain from that?

14. Table 1: The results are self-produced. How comparable are the results actually? What are the best results for the data set published in the literature so far?

15. Training and inference times: Tables summarizing and comparing the results would be nice.

16. I think, a lot of research is going on around Gaussian processes at the moment. I was wondering why there is not so much literature cited (only 16 references). There are also recent papers concerning hyperparameter optimization for Gaussian processes (minimizing the negative log-likelihood).

17. I also want to draw your attention to publications by Potts and Schmischke, who make use of ANOVA-decompositions of functions in terms of Fourier partial sums (https://doi.org/10.1137/20M1354921 and follow-up studies) and using it within least-squares regression / Machine Learning applications.

Overall I have the feeling that a very up-to-date toping is discussed with an innovation that is worth beeing published. However, some things should be clarified better, making the paper easier to read and understand. Advantages, strengths, novelty and innovation should be better discussed and highlighted.

I recommend accepting the paper after a revision.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

**********

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

10.1371/journal.pone.0309661.r002
Author response to Decision Letter 0
Submission Version1
24 Dec 2023

Dear Hanna Landenmark,

On behalf of the authors, thanks to you and the reviewers for assessing our manuscript. We carefully considered all comments and made substantive changes to the manuscript in response. This document summarizes those changes and provides additional explanations where needed.

Journal requirements:

1. We are using a PLoS One LaTeX template. Apologies for missing any details on the submission requirements; we have tried to find and follow all of them in the resubmission. (Where possible: we could not find guidance for file naming in the documents appearing at the links you sent or in the general submission guidelines.)

2. Funding information has been removed from the acknowledgements section. The funding information you presented is correct.

3. An unreviewed preprint has been submitted to arXiv.

4. There are no figures in the manuscript file, only placeholders. However we have removed these as well.

Reviewer 1:

Comment: I found the general idea natural to understand, but the algorithm described in the section "Optimization via forward search" and "Algorithm 1" are VERY DIFFICULT for me to comprehend. The authors should improve the clarity of the writing since the algorithm is the paper's main contribution. For example, can the authors be more explicit about the different stages of the algorithm and how different terms are added? Moreover, can the authors specify how the BIC and AIC are calculated and used? I found these parts poorly explained in the paper, making readers hard to follow what the authors have exactly done.

Response: We have attempted to clarify the algorithm using plain English in addition to the more technical and precise symbolic notation. We note that plain descriptions of the algorithm also appear in the text of the manuscript, and the algorithm itself is designed to add precision to those descriptions. There were indeed several aspects of the algorithm that needed clarification and we the revised version is in fact much clearer. We have also reproduced expressions for BIC and AIC – which, as noted in the text, are the objective functions of the optimization.

Comment: In the explanation of Karhunen-Loeve decomposition on page 4, I am confused about the notations phi_i and u_i. Is u_i normalized? If so, the covariance of (4) seems not equal to k. Please explain. The same question applies to equation (14). Please check.

Response: u_i is normalized by the square root of the eigenvalue, a common approach that renders the variances of the coefficients (betas) iid. However there is a mistake here in the specification of the variance of beta. Thanks for catching it – that mistake has been corrected.

Comment: For equations (16)-(21), I would suggest adding a reference or some details so readers can follow more easily.

Response: This is (mostly) textbook Bayesian linear regression. There are many external resources, including Wikipedia. We added one: a blog post written by one of the authors. That post provides more detail and commentary for readers seeking that.

Comment: Please be more explicit about the terminology used. For example, the author never explained the meaning of "MAE" (e.g., in Table 1)

Response: Sorry for the oversight. MAE is “mean absolute error” and we now define that in the table caption. We also added it to subsequent tables.

Comment: Page 6, equation (11), beta_{dot, k}. What is the meaning of "dot"?

Response: In statistics notation it is common to use a dot as a stand-in for an argument or index when it is arbitrary. Here the index i indicates the input, and the equation indicates that the given relationship holds irrespective of input. However it is also true that the variance of the coefficients is independent of the order of the basis function (as was explained in an earlier response) – we corrected this and now have two dots instead of one. A parenthetical explanation of the notation was also given.

Comment: Page 5, line 59. "supported by" -> "supported on"?

Response: Sure.

Comment: The quality of the figure should be improved. Currently, it is hard to read and compare different lines.

Response: Which figure? In any event we have done our best, producing high-resolution graphics. Since the figures are illustrations of the data appearing in the tables we think these are sufficient.

Reviewer 2:

Comment: I think it should be pointed out more clearly in the abstract already what exactly is new in the method. What was there before and what exactly is your innovation?

Response: From the abstract: “This paper reports a new method of forward variable selection, enabled by the ordered nature of the basis functions in the KL expansion of the Bayesian Smoothing Spline ANOVA kernel (BSS-ANOVA), coupled with fast Gibbs sampling in a fully Bayesian approach.” It goes on to describe what was done to demonstrate this and comparisons to methods that have come before. We are always happy to address comments as best we are able, but this one leaves us at a bit of a loss – we are not sure how we could better accomplish what the reviewer is asking for.

Comment: I think it should be added in the abstract how efficient the proposed method is, in terms of complexity w.r.t. size of the data set, dimensionality etc.

Response: Done.

Comment: Concerning Eq. (3): Some authors mean by the ANOVA decomposition already a uniquely defined decomposition whichs terms are orthogonal w.r.t. to an inner product. In your setting, the single parts in (3) do not neccessarily have to fulfill some orthogonality condition?

Response: In our case, it does not. This was made clear in the revision.

Comment: A lot of parameters are introduced and re-used some pages later wihtout naming them or refering to them. Reading the paper for the first one will face difficulties. Some notations are a somewhat mixed up and not so clear, example: in Eq. (2) it is delta(x|zeta^^2,xi), then in (4) it is delta(x;beta) ... is it the same?

Response: Equations (2) and (4) refer to different cases – equation (2) refers to the conventional squared exponential kernel introduced in (1), whereas (4) refers to the decomposed kernel. The only commonality is the use of delta to represent the GP. The new section and the language introducing it should make clear this is a new case.

Comment: Equation (5): what is the motivation to look at such kernels? What are special properties or advantages in comparison to other frequently applied kernels. Could or should it be used in other frameworks, too?

Response: This is a good question. While any kernel is amenable to a KL decomposition, the kernel hyperparameters don’t always come through the decomposition intact in the prior specification for the expansion coefficients – e.g. the range parameter \\xi of the squared exponential kernel given in (1) – as they do for the BSS-ANOVA kernel. There are undoubtedly other kernels that have this property also. However the support of the BSS-ANOVA kernel is a Sobolev space, as noted in the manuscript, which makes it broad enough for many applications. We do expect that different kernels may be better suited to different types of datasets. The variable selection method works for an arbitrary kernel basis. We’ve added related language to the manuscript.

Comment: Equation (8): It is not so clear to my, why two types of parameters (taus and sigmas) are needed, they are just multiplied.

Response: This is common in Bayesian linear regression, where sigma^2 is the variance of the observation error and sigma^2tau^2 is the variance of the coefficients. This becomes clearer in the section on sampling – we’ve added a note near equation 8.

Comment: Equation (12): This function will be minimized with respect to the betas. Which regularization are you coing to be using? Which optimization algorithms and implementations do you use?

Response: Regularization is L0, as embodied by the BIC / AIC, and the algorithm is precisely the topic of the paper: the forward variable selection routine. We have written out the formulas for the BIC and AIC in response to another reviewer’s comment.

Comment: Equation (14): What are the lambdas?

Response: Eigenvalues of the KL decomposition, as defined in (4).

Comment: Equations (16) and (17): Using the eigenfunctions and constructing the matrix X, you could switch from the kernel formulation to the primal formulation of a Gaussian process. Is this what is expressed by (16) and (17)? Equation (16) is just like kernel ridge regression or linear regression with non-linear feature map.

Response: We are not sure what is meant by the ‘primal formulation of a Gaussian process’, but (16) is linear regression for a fixed number of terms. Ridge regression is variable selection with L2 penalties, which does not make use of a forward variable selection method, and would therefore defeat the purpose here.

Comment: The Section "Variable Selection" is not clearly structured. What is the relation between the different subsections and sub-subsections (not nicely visible)? What is the final result? Please, refer to the Algorithm below. Complexity is not addressed.

Response: See above; we modified the algorithm and the corresponding text to try and make it easier to follow.

Comment: In general, I was wondering why sections and subsections are not numbered.

Response: We used the PLoS LaTeX template.

Comment: ANOVA means analysis of variance. I think, the found betas should give some interpretability (which dimensions are important and which are coupled)!? What insights do you gain from that in the considered numerical examples?

Response: This is a good point. We have added a bit more information (beyond what already appears in the experimental section) regarding the regression results that highlights the interpretability of the KL-decomposed GP.

Comment: Numerical Experiments / Data sets: Please comment a bit on data pre-processing. Your x's have to be scaled to [0,1]. How was that achieved? Already given or via min-max scaling? What is the dimension of the data in the examples? The used hyper parameters have been tuned automatically or tuned/set by hand? What is the outcome of the variable selection procedure? How do the betas look like and which information can you gain from that?

Response: We used min-max scaling. We added that to the manuscript, along with more information about the regression results as mentioned previously. The data dimensionality already appears in the manuscript in the experimental part. For the test cases considered, the model is not overly sensitive to hyperparameters – we added this to the manuscript. They were set by hand – the values used already appear in the manuscript. We added some interpretation of the results based on the estimated coefficients as mentioned above.

Comment: Table 1: The results are self-produced. How comparable are the results actually? What are the best results for the data set published in the literature so far?

Response: At the time the original manuscript was written, we were not aware of previous results for this particular version of the Cascaded Tanks dataset. There are a number of results for another version of the benchmark, which we did not use because of the lack of initial conditions and lack of information for the second state. However, this comment created the incentive to try the search again. We found a number that use the same version of the benchmark as we did, however comparisons are difficult because of the lack of consistent standards for validation in dynamic systems. We note that while all of these studies compared a new method against other methods as we did, none of them included results from other studies. Several do not make predictions against a validation set, reporting only results from training, which is not interesting. One of the studies that used a validation set performed very poorly. Of the two that remained, one (Zabiri et al., 2018, who implement a Support Vector Regression) used a means of validation that separated regions of extrapolation (test sets often overshoot the extremes in the state space encountered in training) from the rest of the test, while also including some training data in the test set. Since the extrapolation regions are where most of the error occurs, and including training data in the test set “pads the stats,” this makes quantitative comparison difficult.

The remaining paper (La Cava, et al., 2016) does not describe how the test set was chosen, but even if one assumes that there is a clear separation between test and training, comparisons are still difficult because of the different methods of testing. As described in the manuscript, we took the standard cross-validation approach more common in the machine learning community, which partitions the whole dataset (from a combination of two separate experiments provided by the benchmark) into five equal subsets, running five train-test studies using each of the partitions as test data in turn, then averaging the test statistics. We would be glad to include any study that used the same method in our tabulated results, however such a study seems not to exist. We did note in the text of the revision that we seem to easily outperform the new method reported in this paper -- Epigenetic Linear Genetic Programming (ELGP), in which an optimization procedure selects parametric terms to use in constructing a dynamic system – by significant margins: our MSE values are over an order of magnitude lower than ELGP, and approximately an order of magnitude lower than a very thinly parameterized recurrent neural network (NARX-NN), which outperformed ELGP.

Comment: Training and inference times: Tables summarizing and comparing the results would be nice.

Response: Though we report and compare training and inference times in the manuscript, simple numerical comparisons are again difficult because the methods we compare against are highly optimized, to include GPU acceleration. In our opinion, tabulated results are only appropriate for “apples-to-apples” comparisons, as many (likely most) readers will at least initially skim papers looking for tables and figures for quick comparisons. Any comparison with significant caveats is therefore best provided in the text of the manuscript. In this case, an at least equally compelling comparison regarding execution time is the number of parameters used in the model, which for us is two orders of magnitude lower than the RNNs, though higher than SINDy.

Comment: I think, a lot of research is going on around Gaussian processes at the moment. I was wondering why there is not so much literature cited (only 16 references). There are also recent papers concerning hyperparameter optimization for Gaussian processes (minimizing the negative log-likelihood).

Response: We have cited recent studies and included one new method (the OAK GP) in our comparisons. However as we note in the manuscript, most scalable GPs are better suited to static problems involving large datasets, rather than for dynamic systems. Since the original manuscript was written we have been made aware of a GP application to dynamic systems. This method uses applies a pair of Gaussian process kernels, one to link a high-dimensional observation space to a low-dimensional latent space, and another to describe dynamics on the latent space (Wang et al., 2005). While the dimensionality reduction is quite interesting it is not something we are addressing in this contribution, and the dynamic modeling seems to have more in common with the RNN methods than with our approach or SINDy. Nonetheless, since some may make the mistake of thinking we are using GPs in this sense, we included the reference.

Incidentally this paper, published at NeurIPS, cites 15 references.

Comment: I also want to draw your attention to publications by Potts and Schmischke, who make use of ANOVA-decompositions of functions in terms of Fourier partial sums (https://doi.org/10.1137/20M1354921 and follow-up studies) and using it within least-squares regression / Machine Learning applications.

Response: Thanks for the suggestion. This is an interesting application to which – like the dynamic GP – seems best suited to learning in situations where there i

Attachment Submitted filename: PLoS_response.docx

10.1371/journal.pone.0309661.r003
Decision Letter 1
Zhou Yu Academic Editor
© 2024 Yu Zhou
2024
Yu Zhou
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version1
5 Jun 2024

PONE-D-23-14436R1Forward variable selection enables fast and accurate dynamic system identification with Karhunen-Loève decomposed Gaussian processesPLOS ONE

Dear Dr. Mebane,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

Please submit your revised manuscript by Jul 20 2024 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Yu Zhou

Academic Editor

PLOS ONE

Journal Requirements:

Please review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript. If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #2: All comments have been addressed

Reviewer #3: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #2: Yes

Reviewer #3: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #2: Yes

Reviewer #3: No

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #2: Yes

Reviewer #3: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #2: Yes

Reviewer #3: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #2: I would like to thank the authors of the paper very much for their insightful responses to my numerous comments and the changes made.

Some of my questions were quite petty, for which I apologize. Despite my scientific proximity to the field, some of my questions arose from a certain lack of knowledge concerning some used basics details as well as pure interest. I hope that my comments have nevertheless been of help in revising the paper.

In my opinion, the authors have responded to all the comments made in the best possible way and have revised the paper accordingly. I recommend accepting the paper in this form.

Just a few small comments (that might be considered for the final revision for publication):

- in the abstract the complexity is stated, but it is not explained what N and P stand for (same page 3, line 18 and following. N=number of points, n=number of dimensions/features?)

- equation (4): one could write in the sum i=1 ... P, then it is immediately clear that P is the number of terms in the expansion

- the style of Figure 1 (size, visability, axis ticks) could be improved

Reviewer #3: The overall presentation of the work is good. Authors presented most of their funding in detail. However a small problem is on the definition of hyper parameters (HP), in some cases a form of optimization is applied to models and in some cases authors stated that “(Results were not particularly sensitive to the hyperparameter values”. It would be a good idea to present a detailed section or at least a table where the author provided full information on HP. HP is a crucial parameter that can significantly affect the overall performance of algorithms. Although authors mentioned that their models' HP are tuned, they did not provide a fully clear information on how they tuned it and whether the HP of counterpart methods are tuned optimally or not. Authors should add a section about this matter and novel approaches in this topic for the possible readers and clearly explain that if the hyper-parameters are not selected accurately or properly even the state-of-the-art methods would fail when compared with a tuned traditional model. 

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #2: No

Reviewer #3: No

**********

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

10.1371/journal.pone.0309661.r004
Author response to Decision Letter 1
Submission Version2
25 Jun 2024

See attached document, 'ResponseToReviewers2.docx'

Attachment Submitted filename: ResponseToReviewers2.docx

10.1371/journal.pone.0309661.r005
Decision Letter 2
Zhou Yu Academic Editor
© 2024 Yu Zhou
2024
Yu Zhou
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version2
16 Aug 2024

Forward variable selection enables fast and accurate dynamic system identification with Karhunen-Loève decomposed Gaussian processes

PONE-D-23-14436R2

Dear Dr. Mebane,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice will be generated when your article is formally accepted. Please note, if your institution has a publishing partnership with PLOS and your article meets the relevant criteria, all or part of your publication costs will be covered. Please make sure your user information is up-to-date by logging into Editorial Manager at Editorial Manager® and clicking the ‘Update My Information' link at the top of the page. If you have any questions relating to publication charges, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Yu Zhou

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #3: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #3: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #3: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #3: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #3: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #3: Authors had made an effort to improve the overall quality of the work. The reviewer has no further comments for this work. The current form of the work can be accepted for publication.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #3: No

**********

10.1371/journal.pone.0309661.r006
Acceptance letter
Zhou Yu Academic Editor
© 2024 Yu Zhou
2024
Yu Zhou
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
11 Sep 2024

PONE-D-23-14436R2

PLOS ONE

Dear Dr. Mebane,

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now being handed over to our production team.

At this stage, our production department will prepare your paper for publication. This includes ensuring the following:

* All references, tables, and figures are properly cited

* All relevant supporting information is included in the manuscript submission,

* There are no issues that prevent the paper from being properly typeset

If revisions are needed, the production department will contact you directly to resolve them. If no revisions are needed, you will receive an email when the publication date has been set. At this time, we do not offer pre-publication proofs to authors during production of the accepted work. Please keep in mind that we are working through a large volume of accepted articles, so please give us a few weeks to review your paper and let you know the next and final steps.

Lastly, if your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

If we can help with anything else, please email us at customercare@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Yu Zhou

Academic Editor

PLOS ONE
==== Refs
References

1 Liu H , Ong YS , Shen X , Cai J . When Gaussian process meets big data: A review of scalable GPs. IEEE Transactions on Neural Networks and Learning Systems. 2020;31 (11 ):4405–4423. doi: 10.1109/TNNLS.2019.2957109 31944966
2 Chalupka K , Williams CK , Murray I . A framework for evaluating approximation methods for Gaussian process regression. Journal of Machine Learning Research. 2013;14 :333–350.
3 Quinonero-Candela J , Rasmussen CE . A unifying view of sparse approximate Gaussian process regression. The Journal of Machine Learning Research. 2005;6 :1939–1959.
4 Deisenroth M, Ng JW. Distributed Gaussian processes. In: International Conference on Machine Learning. PMLR; 2015. p. 1481–1490.
5 Rasmussen C , Ghahramani Z . Infinite mixtures of Gaussian process experts. Advances in Neural Information Processing Systems. 2001;14 .
6 Wang L , Yerramilli S , Iyer A , Apley D , Zhu P , Chen W . Scalable Gaussian processes for data-driven design using big data with categorical factors. Journal of Mechanical Design. 2022;144 (2 ):1–36. doi: 10.1115/1.4052221
7 Wilson AG, Dann C, Nickisch H. Thoughts on massively scalable Gaussian processes. arXiv preprint arXiv:151101870. 2015.
8 Duvenaud D , Nickisch H , Rasmussen CE . Additive Gaussian processes. Advances in Neural Information Processing Systems. 2011;24 .
9 Lu X , Boukouvalas A , Hansman J . Additive Gaussian processes revisited. Proceedings of Machine Learning Research. 2022;162 :14358–14383.
10 Greengard P, O’Neil M. Efficient reduced-rank methods for Gaussian process with eigenfunction expansions. arXiv preprint arXiv:210805924. 2021.
11 Reich BJ , Storlie CB , Bondell HD . Variable selection in Bayesian smoothing spline ANOVA models: Application to deterministic computer codes. Technometrics. 2009;51 (2 ):110–120. doi: 10.1198/TECH.2009.0013 19789732
12 https://scienceforthepeople.substack.com/p/bayesian-linear-calibration-with
13 Bhat KS , Mebane DS , Mahapatra P , Storlie CB . Upscaling uncertainty with dynamic discrepancy for a multi-scale carbon capture system. Journal of the American Statistical Association. 2017;112 (520 ):1453–1467. doi: 10.1080/01621459.2017.1295863
14 Lei Y , Chen TL , Mebane DS , Wen YH . Reduced-order model for microstructure evolution prediction in the electrodes of solid oxide fuel cell with dynamic discrepancy reduced modeling. Journal of Power Sources. 2019;416 :37–49. doi: 10.1016/j.jpowsour.2019.01.046
15 Ostace A , Kocan KX , Mebane DS , Schmal JP , Bhattacharyya D . Probabilistic model building with uncertainty quantification and propagation for a dynamic fixed bed CO2 capture process. Energy & Fuels. 2020;34 (2 ):2516–2532. doi: 10.1021/acs.energyfuels.9b03250
16 Wang J , Hertzmann A , Fleet DJ . Gaussian process dynamical models. Advances in Neural Information Processing Systems. 2005;18 .
17 Brunton SL, Proctor JL, Kutz N. Discovering governing equations from data: Sparse identification of nonlinear dynamical systems. arXiv preprint arXiv:150903580. 2015.
18 Wigren T, Schoukens J. Three free data sets for development and benchmarking in nonlinear system identification. In: 2013 European Control Conference (ECC). European Control Assoc; IFAC; IEEE; CSS; 2013. p. 2933–2938.
19 La Cava W , Danai K , Spector L . Inference of compact nonlinear dynamic models by epigenetic local search. Engineering Applications of Artificial Intelligence. 2016;55 :292–306. doi: 10.1016/j.engappai.2016.07.004
