Bayesian identification of nonseparable Hamiltonians with multiplicative noise using deep learning and reduced-order modeling

The space of real numbers is denoted by $\mathbb{R}$ and the set of positive integers by $\mathbb{Z}_{+}$. Vectors are written in lowercase, non-italic, bold font, e.g., $\mathbf{x}$, and matrices in uppercase, non-italic, bold font, e.g., $\mathbf{A}$. The transpose is denoted by the ^⊤ symbol. If a vector varies in space and/or time, it is spatially indexed as $\mathbf{x}^{s}$ and/or temporally indexed as $\mathbf{x}_{t}$ for discrete space and time $s,t\in\mathbb{Z}_{+}$.

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability triple where $\Omega$ is a sample space, $\mathcal{F}$ is a $\sigma$-algebra, and $\mathbb{P}$ is a probability measure. Random variables are denoted in lowercase, italic, bold font, e.g., $\bm{z}$, and their realizations are denoted as their non-italic counterparts, e.g., $\mathbf{z}$. We assume that for a continuous random variable $\bm{z}$, the probability measure $\mathbb{P}({\rm d}\bm{z})$ admits a probability density function (pdf) $\pi(\mathbf{z})$. A $d$-dimensional uniform distribution with lower and upper bounds $\mathbf{a},\mathbf{b}\in\mathbb{R}^{d}$ is denoted as $\mathcal{U}[\mathbf{a},\mathbf{b}]$. A normal distribution with mean $\mathbf{m}$ and covariance $\mathbf{C}$ is denoted as $\mathcal{N}(\mathbf{m},\mathbf{C})$. If $\bm{z}\sim\mathcal{N}(\mathbf{0},\mathbf{C})$, then $\lvert\bm{z}\rvert$ follows a half-normal distribution denoted as $\text{half-}\mathcal{N}(\mathbf{0},\mathbf{C})$.

The symbol $\odot$ represents element-wise multiplication. The $\odot$ operation is defined when the dimensions of the operands match or when the operands are a matrix and a vector with length equal to the number of matrix columns. In the latter case, $\odot$ multiplies the $i$th column of the matrix by the $i$th element of the vector. This operator has the useful property that $\bm{z}\odot\mathbf{a}\sim\mathcal{N}\left(\mathbf{m}\odot\mathbf{a},\mathbf{C}\odot(\mathbf{a}\mathbf{a}^{\top})\right)$ for the Gaussian random vector $\bm{z}\sim\mathcal{N}(\mathbf{m},\mathbf{C})$ and a constant vector $\mathbf{a}$ of the same length.

In this section, we describe a probabilistic problem formulation for Bayesian system ID for arbitrary noise models. We model the states $\bm{x}_{k}(\omega)\in\mathbb{R}^{n}$ and the outputs $\bm{y}_{k}(\omega)\in\mathbb{R}^{m}$ for $\omega\in\Omega$ as discrete-time stochastic processes indexed by $k\in\mathbb{Z}_{+}$. The dynamics are modeled as a hidden Markov model (HMM)

where $\mathcal{T}:\mathbb{R}^{n}\times\mathbb{R}^{\ell}\times\Omega\mapsto\mathbb{R}^{n}$ is the state-transition function and $\mathcal{M}:\mathbb{R}^{n}\times\mathbb{R}^{\ell}\times\Omega\mapsto\mathbb{R}^{m}$ the measurement function. These functions are parameterized by the random variable $\bm{\theta}(\omega)\in\mathbb{R}^{\ell}$, whose realizations we denote with $\bm{\uptheta}$, and both operators are functions of $\omega\in\Omega$ to represent that their outputs are random variables. The system is therefore characterized by the sequences of transitional pdfs $\pi(\mathbf{x}_{k+1}|\mathbf{x}_{k},\bm{\uptheta})$ and conditional output pdfs $\pi(\mathbf{y}_{k}|\mathbf{x}_{k},\bm{\uptheta})$ induced by $\mathcal{T}$ and $\mathcal{M}$, respectively.

We seek to represent the posterior distribution $\pi(\bm{\uptheta}|\mathcal{Y}_{N})$ characterizing the uncertainty in system parameters $\bm{\theta}$ given a collection of measurements $\mathcal{Y}_{N}\coloneqq(\mathbf{y}_{1},\ldots,\mathbf{y}_{N})$. Bayes’ rule expresses the posterior in a computable form via the likelihood $\mathcal{L}(\bm{\uptheta};\mathcal{Y}_{N})\coloneqq\pi(\mathcal{Y}_{N}|\bm{\uptheta})$, the prior $\pi(\bm{\uptheta})$, and a normalizing constant $\pi(\mathcal{Y}_{N})$ known as the evidence according to

The HMM, however, has the additional collection of uncertain variables $\bm{x}_{1},\ldots,\bm{x}_{N}$ about which we do not intend to make inferences. These uncertain states are marginalized out of the inference problem within the likelihood as $\int\pi(\mathcal{Y}_{N},\mathbf{x}_{1},\ldots,\mathbf{x}_{N}|\bm{\uptheta}){\rm d}\mathbf{x}_{1},\ldots,{\rm d}\mathbf{x}_{N}$. At first, this marginalization appears to require a costly $nN$-dimensional integration. However, the recursive structure of the HMM can be exploited to break the high-dimensional integral into $N$ integrals of the more manageable dimension $n$ by the decomposition

Each term in this product can be efficiently computed using recursion, as shown in Algorithm 1 from [37, Th. 12.3].

Integral evaluation methods within filtering can be divided into two classes: those that estimate the exact integral and those that estimate an approximation of the integral. Estimation of the exact integral is typically performed using sequential Monte Carlo algorithms such as particle filtering [38]. The efficiency of this approach, however, is strongly dependent on the ability to draw uncorrelated samples from the appropriate distributions, which is, in general, nontrivial. When the efficiency of Monte Carlo sampling is prohibitive, approximations are used instead. The most common class of approximation methods for these integrals is Gaussian filtering. These methods approximate the marginal pdf $\mathcal{L}(\bm{\uptheta};\mathcal{Y}_{k})$ as Gaussian, which requires tracking only the first two moments of all other pdfs in Algorithm 1. Tracking the mean and covariance of the prediction pdf $\pi(\mathbf{x}_{k+1}|\mathcal{Y}_{k},\bm{\uptheta})$ and output pdf $\pi(\mathbf{y}_{k}|\mathcal{Y}_{k-1},\bm{\uptheta})$ can be achieved through linearization using Taylor series expansion or through Gaussian integration. These techniques are possible because the prediction and output pdfs are defined in terms of available functions, in this case $\mathcal{T}$ and $\mathcal{M}$. There is no such function, however, for the update pdf $\pi(\mathbf{x}_{k}|\mathcal{Y}_{k},\bm{\uptheta})$. The mean and covariance of this pdf are instead approximated with the Kalman update, which delivers the minimum mean squared error (MMSE) estimate of these quantities.

Here we consider the choice of model parameterization with a focus on embedding geometric structure into the model. The benefits of physics-informed parameterization are that it leads to models that generalize better beyond the training data and reduces the number of data required for training. There are two main ways to enforce the system physics within the dynamics model $\mathcal{T}$. First, the parameterization of $\mathcal{T}$ should be designed to only admit models whose dynamics possess the same physical structure of the system. Second, the model dynamics should be evaluated in such a way that the resulting flow preserves the structure of the dynamics. Since this work considers specialization of methods to nonseparable Hamiltonian systems, we briefly review aspects of the Hamiltonian structure that will inform the proposed parameterization method.

Finite-dimensional canonical Hamiltonian systems are defined by a scalar function $H$, known as the Hamiltonian, of the canonical position $\mathbf{q}$ and momentum $\mathbf{p}$, both in $\mathbb{R}^{d}$. For these systems, the state is defined as $\mathbf{x}=\begin{bmatrix}\mathbf{q}^{\top}&\mathbf{p}^{\top}\end{bmatrix}^{\top}$ such that $n=2d$. The governing equations of the system, known as Hamilton’s equations, are derived from this function

Many Hamiltonian systems of interest to engineers and scientists (e.g., [32, 34, 35, 33, 36]) are not additively separable with respect to functions of the position and momentum. Such Hamiltonian systems are said to be nonseparable. Unlike the separable Hamiltonian systems which can be written as $H(\mathbf{q},\mathbf{p})=T(\mathbf{p})+U(\mathbf{q})$ with a kinetic energy function $T(\mathbf{p})$ and a potential energy function $U(\mathbf{q})$, Eq. (5) cannot be further simplified for nonseparable Hamiltonians. A distinctive feature of Hamilton’s equations is that they possess physically meaningful geometric properties that can be described in the form of symplecticity, invariants of motion, and energy conservation. We discuss how these properties can be embedded in the estimation procedure in Section 4.2.

In this section, we extend the filter for multiplicative scalar noise from [19] to models with statistically-dependent, vector-valued noise. Consider the HMM with additive and multiplicative noise in Eq. (4). Let the vector-valued multiplicative noise terms $\bm{w}_{k}$ and $\bm{v}_{k}$ each be i.i.d. with means $\bar{\mathbf{w}}$ and $\bar{\mathbf{v}}$ and covariances $\mathbf{\Sigma}_{\odot}$ and $\mathbf{\Gamma}_{\odot}$, respectively. Similarly, let the additive noise terms $\bm{\xi}_{k}$ and $\bm{\eta}_{k}$ also be i.i.d. with zero means and covariances $\mathbf{\Sigma}_{+}$ and $\mathbf{\Gamma}_{+}$, respectively. Since we are interested in Gaussian filtering, the higher order moments of these noise terms are not needed.

Recall from Section 3.1 that the goal of Gaussian filtering within Algorithm 1 is to compute (approximations of) the means and covariances of the distributions $\pi(\mathbf{x}_{k}|\mathcal{Y}_{k-1},\bm{\uptheta})$, $\pi(\mathbf{y}_{k}|\mathcal{Y}_{k-1},\bm{\uptheta})$, and $\pi(\mathbf{x}_{k}|\mathcal{Y}_{k},\bm{\uptheta})$. Many of these evaluations require statistics of a function output that can be computed in closed form if the function is linear or approximated with Gaussian approximation if the function is nonlinear. To represent these function outputs, we denote $\Psi(\bm{x}_{k},\bm{\theta})$ as $\bm{\psi}_{k}$ and $h(\mathbf{y}_{k},\bm{\theta})$ as $\bm{h}_{k}$. The following equations outline the filtering procedure with line numbers denoting where each group of equations are evaluated within Algorithm 1. The mean $\mathbf{m}_{k}$ and covariance $\mathbf{C}_{k}$ of $\mathcal{T}(\bm{x}_{k-1},\bm{\theta},\omega)$ with respect to $\pi(\mathbf{x}_{k-1}|\mathcal{Y}_{k-1},\bm{\uptheta})$ (line 6):

The mean $\bm{\upmu}_{k}$ and covariance $\mathbf{S}_{k}$ of $\mathcal{M}(\bm{x}_{k},\bm{\theta},\omega)$ with respect to $\pi(\mathbf{x}_{k}|\mathcal{Y}_{k-1},\bm{\uptheta})$ and the covariance $\mathbf{U}_{k}$ between $\mathcal{T}(\bm{x}_{k-1},\bm{\theta},\omega)$ and $\mathcal{M}(\bm{x}_{k},\bm{\theta},\omega)$ with respect to $\pi(\mathbf{x}_{k},\mathbf{x}_{k-1}|\mathcal{Y}_{k-1},\bm{\uptheta})$ (line 3):

The mean $\mathbf{m}_{k}^{+}$ and covariance $\mathbf{C}_{k}^{+}$ of $\bm{x}_{k}$ with respect to $\pi(\mathbf{x}_{k}|\mathcal{Y}_{k},\bm{\uptheta})$ (line 5):

These last three equations represent the linear MMSE estimator with gain matrix $\mathbf{K}_{k}\in\mathbb{R}^{n\times m}$. Therefore, $\mathbf{m}_{k}^{+}$ and $\mathbf{C}_{k}^{+}$ are only approximations unless $\bm{x}_{k}$ and $\bm{y}_{k}$ are jointly Gaussian. Also notice that $\mathbb{E}\left[\bm{\psi}_{k}\bm{\psi}_{k}^{\top}\right]$ and $\mathbb{E}\left[\bm{h}_{k}\bm{h}_{k}^{\top}\right]$ are required when multiplicative noise is present. Again, these could be computed with linearization or Gaussian integration, but for computational efficiency, we assume that the function outputs are Gaussian. This assumption allows Eq. (7) to be computed in the following form:

A similar form can be used to compute Eq. (9).

In a past work [15], we analyzed the computational complexity of the filtering procedure for additive noise models and found it to be on the order $\mathcal{O}(N(n^{3}+m^{3}))$. Here, we analyze the increase in computational complexity of this filtering procedure when multiplicative noise is present in addition to additive noise. Because we assumed that the noise is stationary, the outer products $\bar{\mathbf{w}}\bar{\mathbf{w}}^{\top}$ and $\bar{\mathbf{v}}\bar{\mathbf{v}}^{\top}$ need only be evaluated once. All other operations scale linearly with the number of data $N$. The added computational cost is summarized in Table 1. The order of the added expense is $\mathcal{O}(N(n^{2}+m^{2}))$, so the additional computation required when multiplicative noise is added to the model does not affect the order of complexity of the overall algorithm.

Next we describe a parameterization strategy for embedding symplectic structure within the learning process that is shown graphically in Fig. 1. The main idea behind this strategy is that rather than parameterizing the time derivatives or the propagator $\Psi$ directly, we parameterize the Hamiltonian $H(\mathbf{q},\mathbf{p},\bm{\theta})$. From this Hamiltonian, the time derivatives are derived from Hamilton’s equations and then integrated with a symplectic integrator. This process ensures that the estimated model will be Hamiltonian and also that its flow will be symplectic. This approach has shown utility in various other works [41, 17, 7, 42].

From the parameterized Hamiltonian, the continuous-time dynamics are given by Hamilton’s equations (5). However, the HMM (4) that we are using is formulated in discrete time. To resolve this discrepancy, we require a propagator $\Psi$ that can map the state forward in time using the estimated time derivatives without violating the physical properties mentioned in Section 3.2. For this purpose, we utilize a recently-developed explicit symplectic integrator for nonseparable Hamiltonian systems [43] that we refer to as “Tao’s integrator” throughout the paper. We briefly describe the integration method here.

Consider an arbitrary Hamiltonian $H(\mathbf{q},\mathbf{p})$ and fictitious position and momentum vectors $\bar{\mathbf{q}}$ and $\bar{\mathbf{p}}$ corresponding to $\mathbf{q}$ and $\mathbf{p}$. Next, define the augmented Hamiltonian

where $H_{a}:=H(\mathbf{q},\bar{\mathbf{p}})$ and $H_{b}:=H(\bar{\mathbf{q}},\mathbf{p})$ correspond to two copies of the original nonseparable Hamiltonian system with mixed-up positions and momenta; $H_{c}$ is an artificial constraint; and $\lambda$ is a constant that controls the binding of the two copies. Unlike the original Hamiltonian $H$, the extended Hamiltonian $\bar{H}$ is amenable to explicit symplectic integration.

From $\bar{H}$, a propagator $\Psi$ is derived using a second-order explicit symplectic method based on Strang splitting

where $\psi^{\Delta t}_{H_{a}},\psi^{\Delta t}_{H_{b}},$ and $\psi^{\Delta t}_{\lambda H_{c}}$ are the time-$\Delta t$ flow of $H_{a},H_{b},$ and $\lambda H_{c}$. Each flow is evaluated with explicit symplectic Euler substeps, so the result of the composition is an explicit and symplectic propagator as desired.

Reducing the dimension of a dynamical system hinges upon a key hypothesis: that the system of interest is a high-dimensional realization of underlying low-dimensional dynamics. Under this assumption, it is theoretically possible to estimate dynamics with a considerably lower state dimension than the ambient dimension without any loss of accuracy. Estimating the hypothesized low-dimensional dynamics involves two primary components. The first is identifying a low-dimensional subspace on which a significant portion of the state evolution is contained. The other component is estimating a dynamical model whose state evolves within this low-dimensional space. Here we first describe linear dimension reduction for dynamical systems and then expound on alterations that can be made to preserve symplectic structure.