Learning Dynamics from Input-Output Data with Hamiltonian Gaussian Processes

1. Problem Statement

The paper addresses the challenge of learning physically consistent, non-conservative dynamical systems from input-output data only, without access to full state measurements (specifically, momentum or velocity).

• Context: Physics-informed machine learning (PIML) often incorporates Hamiltonian mechanics to ensure energy conservation and physical interpretability.
• Limitation of Existing Work: Current Hamiltonian Gaussian Process (GP) methods (e.g., Beckers et al., 2022) require training data that includes momentum ($p$) or velocity ($\dot{q}$) in addition to position ($q$).
  • In many practical applications, only position (input-output) data is available.
  • Approximating velocity via numerical differentiation is highly sensitive to noise.
  • Generalized momenta cannot always be constructed from position data alone.
• Goal: Develop a method to learn uncertainty-quantifying, energy-consistent (or dissipative) dynamics models using only input ($u$) and output ($y$, typically position) data, while simultaneously estimating unknown structural parameters (e.g., damping coefficients) and latent states.

2. Methodology

The proposed approach combines Reduced-Rank Gaussian Processes with Hamiltonian Mechanics within a Fully Bayesian Inference framework.

A. System Formulation

The system is modeled using forced Hamiltonian mechanics:
$$\dot{x} = (J(x) - R(x)) \nabla_x H(x) + G(x)u + w$$
$$y = g(x) + e$$
Where:

• $x = [q, p]^\top$ is the state (positions and momenta).
• $H(x)$ is the unknown Hamiltonian function (energy landscape).
• $J(x)$ is the skew-symmetric interconnection matrix (conservative dynamics).
• $R(x)$ is the symmetric positive semi-definite dissipation matrix (energy loss).
• $G(x)$ is the input matrix.
• $w$ and $e$ are process and measurement noises.
• Key Assumption: The structural forms of $J, R, G$ are known, but their specific parameters (e.g., damping coefficients) and the function $H$ are unknown.

B. Reduced-Rank Hamiltonian GP

To avoid the computational intractability of standard GPs (which scale cubically with data size, $O(T^3)$), the authors employ a Reduced-Rank GP (Solin and Särkkä, 2020).

• Approximation: The Hamiltonian $H(x)$ is approximated using a finite-dimensional eigenfunction expansion of the kernel's spectral density:
  $$\hat{H}(x) \approx \sum_{k=1}^M a_k \phi_k(x) = a^\top \phi(x)$$
• Gradient Modeling: Since the dynamics depend on $\nabla_x H(x)$, the gradient is modeled as a linear transformation of the basis functions: $\nabla_x \hat{H}(x) = D_\phi(x) a$.
• Advantages:
  • Linearity: The model is linear in parameters $a$, enabling closed-form updates.
  • Complexity: Computational cost is independent of the number of training data points ($O(M)$ for prediction, where $M$ is the number of basis functions), making it scalable.
  • Noise Handling: An auxiliary variable $h_t = \nabla_x \hat{H}(x_t) + v_t$ is introduced to decouple learning the Hamiltonian from inferring hidden states.

C. Fully Bayesian Inference (Particle Gibbs)

The authors propose a Particle Markov Chain Monte Carlo (PMCMC) scheme, specifically a Particle Gibbs algorithm, to estimate the joint posterior distribution of latent states ($x, h$) and parameters ($\theta, \vartheta$).

• Parameters:
  • $\theta$: GP weights and noise scale.
  • $\vartheta_K$: Kernel hyperparameters (length scale, variance).
  • $\vartheta_S$: Structural hyperparameters (e.g., damping coefficients in $R$).
• Inference Steps:
  1. State Inference (Step i): Use Conditional Sequential Monte Carlo (SMC) to sample latent trajectories $z_{0:T} = \{x_{0:T}, h_{0:T}\}$ given current parameters.
  2. Parameter Learning (Step ii):
     • Closed-form update: For $\theta$ (GP weights), the posterior is derived analytically using conjugate Normal-Inverse-Gamma (NIG) priors.
     • Metropolis-Hastings (MH): For non-conjugate hyperparameters $\vartheta_K$ and $\vartheta_S$, MH steps are used within the Gibbs sampler.
• Result: The method provides full probability density estimates for hidden states, model parameters, and structural hyperparameters, rather than just point estimates.

3. Key Contributions

1. Input-Output Only Learning: The first Hamiltonian GP framework capable of learning non-conservative dynamics solely from input-output data, eliminating the need for velocity/momentum measurements.
2. Reduced-Rank Efficiency: Introduction of a reduced-rank GP formulation that ensures computational complexity is independent of the dataset size, addressing the scalability issues of standard Hamiltonian GPs.
3. Fully Bayesian Structural Learning: A unified scheme that estimates not only the dynamics but also structural hyperparameters (like damping coefficients) and kernel hyperparameters with uncertainty quantification, moving beyond frequentist point estimates used in prior work.
4. Non-Conservative Modeling: Explicitly handles energy dissipation and external forcing, making the model applicable to real-world systems that are not perfectly conservative.

4. Results

The method was evaluated on a non-harmonic oscillator simulation ($H(q, p) = q^2/2 + p^2/2 + 2\cos q$) with unknown damping.

• Accuracy:
  • The proposed method (Input-Output) achieved a Flow Magnitude RMSE of 3.27 J and Flow Angle RMSE of 0.12 rad.
  • This is comparable to an "Exact GP" trained on full state data (Input-State), which achieved 3.09 J and 0.08 rad respectively.
  • It significantly outperformed standard approaches that fail without momentum data.
• Uncertainty Quantification:
  • The Bayesian framework successfully inferred the posterior densities of the damping coefficient and kernel hyperparameters.
  • While the damping estimate from input-output data showed slight bias (due to identifiability challenges in the specific simulation setup), the method correctly captured the multimodal nature of the system's phase space and energy dissipation.
• Physical Consistency:
  • Forward simulations using sampled models from the posterior closely matched true system trajectories.
  • The models correctly captured energy decay (dissipation) when inputs were removed.
• Computational Efficiency:
  • The reduced-rank approach maintained $O(M)$ prediction complexity, whereas exact GPs scale as $O(T^3)$.

5. Significance and Future Work

• Bridging the Gap: This work bridges the gap between physics-informed learning (which usually requires full state data) and practical system identification (where only input-output data is available).
• Control Applications: By providing uncertainty-quantified, physically consistent models from limited data, the method enables safer and more robust model-based control (e.g., stochastic MPC) for complex systems where sensors for velocity/momentum are unavailable or noisy.
• Limitations & Future Directions:
  • Identifiability: The simulation required a symmetry constraint to ensure the Hamiltonian was identifiable from input-output data. Future work needs to formalize identifiability conditions for general systems.
  • Dimensionality: The computational cost of the particle-based scheme and the basis expansion grows with the system dimension. Future research should address high-dimensional systems.
  • Symmetry Constraints: The need for manual symmetry constraints suggests a need for automated structure learning or better prior formulations.

In summary, the paper presents a robust, scalable, and fully Bayesian framework for learning complex, dissipative dynamical systems from noisy input-output data, overcoming a major bottleneck in applying Hamiltonian machine learning to real-world engineering problems.