Bayesian Regularization for Dynamical System Identification: Additive Noise Models

Consider the dynamical system <inline-formula> <mml:math id="mm1"> <mml:semantics> <mml:mrow> <mml:mover accent="true"> <mml:mi mathvariant="bold-italic">x</mml:mi> <mml:mo>˙</mml:mo> </mml:mover> <mml:mo>=</mml:mo> <mml:mi>f</mml:mi> <mml:mrow> <mml:mo>(</mml:mo> <mml:mi mathvariant="bold-italic">x</mml:mi> <mml:mo>)</mml:mo> </mml:mrow> </mml:mrow> </mml:semantics> </mml:math> </inline-formula>, where <inline-formula> <mml:math id="mm2"> <mml:semantics> <mml:mrow> <mml:mi mathvariant="bold-italic">x</mml:mi> <mml:mo>∈</mml:mo> <mml:msup> <mml:mi mathvariant="double-struck">R</mml:mi> <mml:mi>n</mml:mi> </mml:msup> </mml:mrow> </mml:semantics> </mml:math> </inline-formula> is the state vector, <inline-formula> <mml:math id="mm3"> <mml:semantics> <mml:mover accent="true"> <mml:mi mathvariant="bold-italic">x</mml:mi> <mml:mo>˙</mml:mo> </mml:mover> </mml:semantics> </mml:math> </inline-formula> is the time or spatial derivative, and <i>f</i> is the system model. We wish to identify unknown <i>f</i> from its time-series or spatial data. For this, we propose a Bayesian framework based on the maximum a posteriori (MAP) point estimate, to give a generalized Tikhonov regularization method with the residual and regularization terms identified, respectively, with the negative logarithms of the likelihood and prior distributions. As well as estimates of the model coefficients, the Bayesian interpretation provides access to the full Bayesian apparatus, including the ranking of models, the quantification of model uncertainties, and the estimation of unknown (nuisance) hyperparameters. For multivariate Gaussian likelihood and prior distributions, the Bayesian formulation gives a Gaussian posterior distribution, in which the numerator contains a Mahalanobis distance or “Gaussian norm”. In this study, two Bayesian algorithms for the estimation of hyperparameters—the joint maximum a posteriori (JMAP) and variational Bayesian approximation (VBA)—are compared to the popular SINDy, LASSO, and ridge regression algorithms for the analysis of several dynamical systems with additive noise. We consider two dynamical systems, the Lorenz convection system and the Shil’nikov cubic system, with four choices of noise model: symmetric Gaussian or Laplace noise and skewed Rayleigh or Erlang noise, with different magnitudes. The posterior Gaussian norm is found to provide a robust metric for quantitative model selection—with quantification of the model uncertainties—across all dynamical systems and noise models examined.