MS-PHS GP: Multistep Port-Hamiltonian GP

• MS-PHS GP is a framework for learning physically consistent continuous-time dynamics by placing a Gaussian process prior directly on the Hamiltonian surface.
• It employs a structured vector-field kernel induced from the Hamiltonian and incorporates variable-step multistep integrator constraints as finite linear functionals.
• The method achieves closed-form Gaussian conditioning without latent states, ensuring energy balance, passivity, and providing finite-sample error guarantees.

Searching arXiv for papers on MS-PHS GP and related port-Hamiltonian Gaussian process methods. Multistep Port-Hamiltonian Gaussian Process (MS-PHS GP) is a framework for learning physically consistent continuous-time dynamics and a posterior over the Hamiltonian surface from noisy, irregularly-sampled trajectories. It places a Gaussian process (GP) prior on the Hamiltonian surface and encodes variable-step multistep integrator constraints as finite linear functionals, which enables closed-form conditioning of both the vector field and the Hamiltonian surface without latent states, while enforcing energy balance and passivity by design (Leung et al., 1 Oct 2025).

1. Port-Hamiltonian setting and problem formulation

MS-PHS GP is formulated for port-Hamiltonian systems (PHS), which model energy-conserving and dissipative dynamics with input/output ports. In the notation used for the method,

$$\dot{x} = [J({x}) - R({x})] \nabla_x H({x}) + G({x}) u, \qquad y_{\text{out}} = G(x)^\top \nabla_x H(x),$$

where $x \in \mathbb{R}^n$ is the state, $H$ is the Hamiltonian, $J(x)$ is a skew-symmetric interconnection matrix, $R(x)$ is a positive semi-definite dissipation matrix, and $G(x)$ is the port mapping (Leung et al., 1 Oct 2025).

The learning problem is posed from noisy trajectory data of the form

$$\mathcal{D} = \{(\tilde{x}(t_k), u(t_k), t_k)\}_{k=1}^K, \qquad \tilde{x}(t_k) = x(t_k) + \varepsilon_k.$$

The specific objective is not only to recover a vector field, but also to infer the Hamiltonian surface itself, together with posterior uncertainty. Within this formulation, energy balance and passivity are not imposed a posteriori; they are encoded in the model class.

A useful antecedent is the GP-PHS framework, which models the unknown Hamiltonian $H(x)$ nonparametrically using a GP and proves that all samples from the resulting model define systems that are passive with respect to the input-output pair. GP-PHS also preserves the compositional nature of Port-Hamiltonian systems (Beckers et al., 2023). MS-PHS GP inherits that port-Hamiltonian perspective, but changes the observation model by embedding multistep discretization constraints directly into the GP conditioning procedure.

2. Gaussian-process prior on the Hamiltonian and induced vector-field kernel

The central modeling choice is to place the GP prior directly on the Hamiltonian: $$H(x) \sim \mathcal{GP}\bigl(0, k_{\rm base}(x, x')\bigr).$$ By closure of GPs under linear, and in particular differential, transformations, this prior induces a matrix-valued GP prior on the vector field: $$k_{\text{phs}}(x, x') = \sigma_f^2 J_R(x)\nabla_x\nabla_{x'} k_{\rm base}(x, x') J_R(x')^\top,$$ where $x \in \mathbb{R}^n$0. The induced vector-field prior is then

$x \in \mathbb{R}^n$1

This construction ensures that sampled vector fields remain in port-Hamiltonian form (Leung et al., 1 Oct 2025).

The significance of this choice is structural rather than merely statistical. Because the unknown object is the Hamiltonian surface, and because the vector field is obtained through the port-Hamiltonian differential operator, the model couples inference and physical admissibility at the prior level. In the earlier GP-PHS formulation, the same design principle was used to guarantee passivity and compositionality for learned finite-dimensional systems (Beckers et al., 2023). MS-PHS GP retains that physics-informed prior, but augments it with a multistep observation operator tailored to irregular time grids.

3. Variable-step multistep constraints as finite linear functionals

The distinctive feature of MS-PHS GP is the encoding of a variable-step multistep integrator as a finite linear functional over the unknown dynamics. For a multistep scheme such as Adams-Bashforth, the approximate ODE flow over the data is expressed as

$x \in \mathbb{R}^n$2

where $x \in \mathbb{R}^n$3 and $x \in \mathbb{R}^n$4 are integrator-dependent matrices that depend on step sizes and method order, $x \in \mathbb{R}^n$5, and $x \in \mathbb{R}^n$6. For system dimension $x \in \mathbb{R}^n$7, the construction is extended via Kronecker products, yielding $x \in \mathbb{R}^n$8 and corresponding block operators (Leung et al., 1 Oct 2025).

With noisy observations, the vectorized data constraint becomes

$x \in \mathbb{R}^n$9

The key insight is that linear constraints on the vector field and derivative observations translate into finite-dimensional linear functionals on the GP prior. Since GPs are closed under conditioning on linear functionals, including derivative information, this yields closed-form Gaussian posteriors over both the vector field and the Hamiltonian without the need to explicitly introduce or marginalize over latent state variables (Leung et al., 1 Oct 2025).

The multistep kernel projection is written as

$H$0

The posterior mean and covariance for the vector field at a test point are then

$H$1

$H$2

All quantities are available in closed form; no sampling, filtering, or latent ODE solvers are required (Leung et al., 1 Oct 2025).

A common misconception is to regard MS-PHS GP as a derivative-preprocessing pipeline followed by ordinary GP regression. The framework is explicitly designed to handle irregular sampling and noise natively, with no need for data filtering or estimation of derivatives as a separate preprocessing step (Leung et al., 1 Oct 2025).

4. Posterior over the Hamiltonian surface

MS-PHS GP provides a posterior not only over the vector field, but also over the Hamiltonian surface itself. To fix the integration constant, the construction includes an anchor point, for example $H$3. One then forms a joint Gaussian distribution over

$H$4

with covariance matrix $H$5 built from $H$6, its derivatives, and the projected kernel terms (Leung et al., 1 Oct 2025).

At a test point $H$7, the Hamiltonian posterior is

$H$8

with

$H$9

The resulting posterior uncertainty is intended to be physically interpretable as uncertainty over the energy landscape itself, not merely over short-horizon derivatives (Leung et al., 1 Oct 2025).

This construction is central to the method’s interpretation. The model does not first infer latent derivatives and then integrate them into an energy function. Instead, the Hamiltonian prior and the multistep linear functional are combined in one Gaussian conditioning problem. This suggests that MS-PHS GP is best understood as a port-Hamiltonian kernel method in which the numerical discretization is part of the observation operator.

5. Finite-sample guarantee and empirical behavior

The method states a finite-sample vector-field bound. Under the assumption that $J(x)$0 lies in the base RKHS induced by $J(x)$1, for any $J(x)$2, with probability at least $J(x)$3,

$J(x)$4

where $J(x)$5, $J(x)$6, and $J(x)$7 are problem/integrator-dependent constants, $J(x)$8 is the number of samples, $J(x)$9 is the time step, and $R(x)$0 is the integrator order (Leung et al., 1 Oct 2025).

The interpretation given for this result is that the error splits into a statistical estimation term, which improves with more data and scales as $R(x)$1 in the worst case, and a discretization bias term, which improves with smaller step sizes and higher-order integrators, scaling as $R(x)$2. Increasing data alone cannot eliminate discretization bias; accurate, calibrated vector field and Hamiltonian inference requires embedding the numerical integrator jointly with the GP (Leung et al., 1 Oct 2025).

Empirical evaluation was carried out on a linear mass-spring system, the Van der Pol oscillator, and the Duffing oscillator, all under noisy, irregular sampling schemes. Comparative baselines included MS-PHS-ab with Adams-Bashforth orders 1, 2, and 3; MS-ODE-ab, which performs multistep GP inference without a PHS prior; and GP-PHS-savgol/loess, which relies on Savitzky-Golay or LOESS preprocessing for derivative estimation (Leung et al., 1 Oct 2025).

On these benchmarks, MS-PHS matches or outperforms other methods in vector-field MSE and directional alignment, and higher-order integrators help further. The physics-informed prior regularizes extrapolation and prevents spurious behaviors in data-sparse regions. For Hamiltonian recovery, the posterior mean more accurately reconstructs the underlying energy landscape, and the predictive variance of $R(x)$3 tracks empirical error well. By contrast, preprocessing-based baselines can significantly underestimate uncertainty, especially as noise increases. For the error-uncertainty calibration ratio $R(x)$4, MS-PHS stays near unity across noise and time-grid jitter levels, whereas other baselines rise above unity (Leung et al., 1 Oct 2025).

6. Relation to neighboring methods and subsequent developments

MS-PHS GP sits at the intersection of two research directions: structure-preserving GP dynamics and GP-based port-Hamiltonian learning. Structure-preserving Gaussian Process Dynamics combines GPs with structure-preserving numerical integrators such as Runge-Kutta methods, including implicit schemes, and evaluates the implicitly defined Runge-Kutta transition probability through an implicit layer embedded in variational inference. That framework supports both multistep and one-step structure-preserving schemes, but it is formulated as a recurrent variational GP state-space model with sampled GP transitions (Ensinger et al., 2021). By contrast, MS-PHS GP uses variable-step multistep integrator constraints as finite linear functionals and obtains closed-form Gaussian conditioning without latent states (Leung et al., 1 Oct 2025).

A concise comparison of adjacent constructions is given below.

Method	Core construction	Distinctive feature
GP-PHS (Beckers et al., 2023)	GP prior on $R(x)$5 within PHS dynamics	Passivity and compositionality
Structure-preserving GP dynamics (Ensinger et al., 2021)	GP-learned vector field inside structure-preserving integrator	Implicit layer and variational inference
MS-PHS GP (Leung et al., 1 Oct 2025)	GP prior on Hamiltonian plus variable-step multistep constraints	Closed-form conditioning without latent states
PFEM-GP-dpHs (Courteville et al., 15 Dec 2025)	GP prior on Hamiltonian functional plus PFEM	Late lumping for distributed systems

Later work broadens the same design philosophy. PFEM-GP-dpHs combines the partitioned finite element method with GP learning for infinite-dimensional distributed port-Hamiltonian systems, uses a late lumping approach to decouple discretization of GP hyperparameters from discretization of the distributed dynamics, and models the mean of the GP prior of the Hamiltonian as a quadratic form so that the GP kernel focuses on the nonlinear part of a given distributed system (Courteville et al., 15 Dec 2025). In a different direction, data-driven identification of port-Hamiltonian DAE systems uses multi-task Gaussian processes to estimate nonlinear effort functions for index-one and index-three pH-DAEs, extending GP-based identification beyond ordinary differential equation models (Zaspel et al., 2024).

Control-oriented developments also exploit the uncertainty quantification furnished by Hamiltonian GP models. Gaussian Process Port-Hamiltonian Systems underlie a GP-based control approach that uses passivity-based robust control with interconnection and damping assignment to establish probabilistic stability guarantees (Beckers, 2023). For distributed systems, GP-dPHS learning has been combined with boundary control by interconnection, and the posterior uncertainty has been incorporated into an energy-based robustness analysis yielding probabilistic conditions under which the closed-loop trajectories remain bounded despite model mismatch (Beckers et al., 5 Apr 2026). A plausible implication is that the MS-PHS GP viewpoint—namely, learning Hamiltonian structure while respecting the numerical and physical operators that connect data to dynamics—provides a template for broader structure-preserving Bayesian system identification beyond finite-dimensional, fully observed settings.