Hyperparameter Optimization in the Estimation of PDE and Delay-PDE models from data

We propose an improved method for estimating partial differential equations and delay partial differential equations from data, using Bayesian optimization and the Bayesian information criterion to automatically find suitable hyperparameters for the method itself or for the equations (such as a time-delay). We show that combining time integration into an established model estimation method increases robustness and yields predictive models. Allowing hyperparameters to be optimized as part of the model estimation results in a wider modelling scope. We demonstrate the method’s performance on a number of synthetic benchmark problems of different complexity, representing different classes of physical behaviour. This includes the Allen-Cahn and Cahn-Hilliard models, as well as different reaction-diffusion systems without and with time-delay.

💡 Research Summary

The paper presents a novel framework for data‑driven identification of partial differential equations (PDEs) and delay‑PDEs that automatically tunes the hyper‑parameters governing the learning algorithm. The authors build on the sparse regression paradigm (e.g., SINDy) but address two major shortcomings of existing approaches: (i) the need to manually set sparsity thresholds and other regularisation parameters, and (ii) the lack of an integrated assessment of model quality that accounts for numerical integration errors.

The methodology proceeds as follows. A large library Θ(u) of candidate functions is constructed; it may contain polynomial terms, nonlinear functions, spatial derivatives, and optionally delayed terms with a time‑lag τ. The true dynamics are assumed to be of the form ∂ₜu = σ·Θ(u), where σ is a coefficient matrix. An initial estimate of σ is obtained by ordinary least‑squares on the time derivative of the data (computed via finite differences). Because plain least‑squares yields dense solutions, the authors employ Sequential Thresholded Least Squares (STLS) to enforce sparsity: coefficients whose absolute value falls below a threshold h are set to zero and the regression is repeated until convergence. Crucially, the thresholds h can be distinct for each state variable or for groups of terms (e.g., a separate h_τ for delayed terms).

Rather than fixing h (and any other hyper‑parameters such as τ) a priori, the paper embeds them in a Bayesian optimisation loop using the Hyperopt library and its Tree‑structured Parzen Estimator (TPE). The objective function is the Bayesian Information Criterion (BIC) defined as BIC = s·ln(Nₜ) – 2·ln(ℒ̂), where s is the number of non‑zero coefficients, Nₜ the number of temporal samples, and ℒ̂ the likelihood associated with the integrated model error. To evaluate ℒ̂, the candidate PDE is integrated forward in time (SciPy’s explicit Runge‑Kutta of order 8, with a fallback trapezoidal scheme for stability) and the L₂ distance between the simulated trajectory û and the original data u is computed. By minimising BIC, the optimisation simultaneously penalises model complexity and rewards accurate time‑integrated predictions, thereby avoiding over‑fitting that plagues pure derivative‑based residual minimisation.

The authors test the approach on several synthetic benchmarks. (1) A complex Ginzburg‑Landau reaction‑diffusion system is used to compare against the PDE‑SINDy implementation of Ref.