Port-Hamiltonian Neural ODE
Source paper: Model reduction of port-Hamiltonian systems via neural networks arXiv:2608.30788 ⓘ · analyzed Sep 1, 2026
AI-generated research hypothesis, automatically tested. Not peer-reviewed.
Idea description
Replace an unconstrained neural ODE vector field with a learned port-Hamiltonian vector field whose energy gradient drives the dynamics, whose interconnection matrix is skew-symmetric, and whose dissipation matrix is positive semidefinite. The resulting model remains expressive through state-dependent neural matrices while guaranteeing non-increasing learned energy in the unforced case.
Formulas
Mathematical statement
The paper models the state x(t) in R^N with the port-Hamiltonian equations dx/dt = (J(x)-R(x))∇_x H(x)+Bu(t) and y(t)=B^T∇_x H(x), where H:R^N→R is the Hamiltonian or stored energy, J(x)∈R^{N×N} is the interconnection matrix, R(x)∈R^{N×N} is the dissipation matrix, B is the input-port matrix, u(t) is the input, and y(t) is the output. The structural assumptions are J(x)^T=-J(x) and R(x)⪰0 for every x. For u=0, the energy derivative is dH(x(t))/dt = ∇H(x)^T(J(x)-R(x))∇H(x) = -∇H(x)^T R(x)∇H(x)≤0, because v^T Jv=0 for every vector v and v^T Rv≥0 when R is positive semidefinite. Implement these constraints by predicting an unconstrained A_theta(z) and L_theta(z), then setting J_theta(z)=A_theta(z)-A_theta(z)^T and R_theta(z)=L_theta(z)L_theta(z)^T+epsilon I with epsilon>0. The latent state z∈R^d replaces x, H_theta(z) is a scalar neural energy, and B_theta may be fixed or learned. The explicit factorization is the transferable guarantee; it avoids penalty-only enforcement of symmetry or stability.
Implementation notes
1. Integration point: use this as the vector field of a latent neural ODE or continuous-time state-space model. Encode an observation sequence into an initial latent state z0∈R^d, integrate z with an ODE solver, and decode z(t) to predictions. At each function evaluation, feed z and optionally u(t) into three MLP heads: a scalar energy head h_theta(z), a matrix head A_theta(z), and a factor head L_theta(z). Define H_theta(z)=softplus(h_theta(z))+0.5 z^T Q z, with Q=Q_raw Q_raw^T+delta I and delta>0, so the learned energy is bounded below and locally coercive. 2. Pseudocode: g = grad_z(H_theta(z)); A = A_net(z); J = A - transpose(A); L = L_net(z); R = L @ transpose(L) + eps*I; dz = (J-R) @ g + B @ u; return dz. Use fixed-step RK4 first, then compare with an adaptive solver. 3. Computed exactly from the mathematics: J is skew-symmetric, R is positive semidefinite, and the unforced energy derivative is nonpositive. Estimate empirically the trajectory loss, energy drift, Jacobian spectral abscissa, solver function evaluations, and long-horizon rollout error. For forced data, optionally add the balance residual dH/dt - (y^T u - g^T R g); do not use symmetry or PSD penalties because the parameterization enforces them exactly. 4. First cheap experiment: fit a 2D or 4D nonlinear mass-spring-damper benchmark and compare against an equally sized unconstrained neural ODE and a standard Hamiltonian neural network. Train on short trajectories from random initial states with random forcing, then test rollouts 10–50 times longer. Success means lower long-horizon state error, no energy growth on unforced rollouts, and comparable one-step error at similar parameter count. Report wall-clock cost and ODE function evaluations as well.
Verification
Mechanism evidence: Not confirmed
Practical benchmark: Not run
Stage 1 — Mechanism check agent confidence 9/10
Built a port-Hamiltonian neural vector field with exact skew-symmetric J=A-A^T, PSD R=LL^T+epsilon I, and an autograd-computed learned energy gradient. The mechanism manifested: skew residual was exactly 0, the minimum R eigenvalue stayed at or above epsilon=0.03, and the continuous energy identity matched numerically to 3.3e-16 with no positive derivatives. The dissipation sweep confirmed the predicted quadratic dependence on factor scale, and the PH model had lower long-horizon RMSE than the unconstrained baseline, though with higher parameter count and training cost.
- Agent confidence
- 9/10
- Baseline
- Unconstrained neural vector field: long-horizon RMSE 0.5380; 1,218 parameters; training time 0.561 s.
- Idea
- Port-Hamiltonian neural ODE: long-horizon RMSE 0.1530; 3,757 parameters; training time 2.724 s; maximum skew residual 0; minimum R eigenvalue 0.030000; energy identity residual 3.33e-16; RK4 energy increases 0.
- Predictions
- J antisymmetry predicted exact and observed maximum residual 0. R eigenvalue predicted >= epsilon=0.03 and observed minimum 0.0300000. Dissipation dH/dt was predicted affine in lambda^2: slope -19.3686 predicted versus -19.4580 observed.
Limitations: Only a small 2D damped nonlinear oscillator was tested, with fixed-step RK4, short training, and unequal parameter counts. Adaptive solvers, forced balance residuals, a standard Hamiltonian baseline, FLOP-matched training, and larger datasets were not tested. Discrete RK4 energy non-increase is empirical rather than an unconditional guarantee.
How to run: python3 experiment.py
Verdict computed by deterministic test code from paired-seed statistics — not by the language model.
Artifacts
- experiment.py 6.3 KB View Raw JSON
- report.md 1.8 KB View
- stage2_bench.py 5.3 KB View Raw JSON