Symplectic Neural Flows

• Symplectic neural flows are neural architectures that enforce canonical symplectic structure, preserving phase-space volume and Hamiltonian dynamics.
• They employ design patterns like Real-NVP layers, symplectic integrators, and adjoint methods to maintain invariant-preserving maps during training.
• Applications span canonical transformations, density modeling, and nonseparable Hamiltonians, proving effective in latent mode discovery and trajectory predictions.

Symplectic neural flows are neural architectures and learning schemes that parameterize maps or time-dependent flows on phase space while enforcing canonical symplectic structure. In canonical coordinates $x=(q,p)\in\mathbb R^{2n}$, this typically means requiring the Jacobian of a map $\Phi$ to satisfy

$\bigl(D\Phi(x)\bigr)^{T}J\bigl(D\Phi(x)\bigr)=J,\qquad J=\begin{pmatrix}0&I_n\-I_n&0\end{pmatrix},$

or, in continuous time, constructing dynamics through Hamilton’s equations $\dot x=J\nabla_x H$. The resulting models are used for canonical transformation, structure-preserving neural integration, density modeling, collective-variable discovery, and learning of Hamiltonian or volume-preserving dynamics. The literature represented here includes neural canonical transformations to latent harmonic oscillators, explicit and implicit symplectic neural integrators, exact symplectomorphism networks, locally-symplectic architectures for divergence-free systems, and augmented schemes for nonseparable Hamiltonians (Li et al., 2019, Maslovskaya et al., 2024, Canizares et al., 2024).

1. Canonical geometry and the defining constraint

The common mathematical substrate is the phase space $(q,p)$, the canonical two-form $\omega=\sum_i dq_i\wedge dp_i$, and Hamilton’s equations

$$\dot q_i=\frac{\partial H}{\partial p_i},\qquad \dot p_i=-\frac{\partial H}{\partial q_i}.$$

A symplectic map preserves $\omega$, equivalently satisfies $D\Phi^\top J D\Phi=J$, and therefore preserves phase-space volume in the Liouville sense. Several works treat this as the central architectural invariant rather than as a soft regularizer (Li et al., 2019, He et al., 2024).

Within this framework, different papers emphasize different targets. In neural canonical transformation, the goal is to map physical variables $(q,p)$ to latent variables $\Phi$0 whose Hamiltonian is a decoupled harmonic oscillator,

$\Phi$1

so that the latent density becomes a factorized Gaussian and the learned frequencies $\Phi$2 order slow and fast modes (Li et al., 2019). In symplectic neural integrators, the target is instead a time-$\Phi$3 approximation to Hamiltonian flow, often by composing elementary exact or symplectic substeps (Maslovskaya et al., 2024, Canizares et al., 2024).

Related work broadens the scope beyond canonical Hamiltonian systems in the strict narrow sense. Locally-symplectic networks begin from the Feng–Shang theorem that a divergence-free vector field can be decomposed into pairwise Hamiltonian subsystems, and then learn the resulting flow by composing local symplectic modules (Bajārs, 2021). This suggests that, in current usage, symplectic neural flows form a geometric modeling family centered on symplecticity but not restricted to one architectural template.

2. Architectural constructions

The principal construction patterns are summarized below.

Framework	Core symplectic mechanism	Target setting
Neural Canonical Transformation	Stacked point transformations with Real-NVP coordinate maps and induced momentum update	Canonical transformation and latent mode discovery
SPRK Nets	Higher-order explicit symplectic Partitioned Runge–Kutta layers	Deep architectures and learned dynamical structure
SpNN	Composition of q-shearing, p-shearing, and symplectic stretching blocks	Exact invertible symplectomorphisms
Taylor-net	Symmetric Taylor-series subnetworks plus fourth-order symplectic integrator	Continuous long-term prediction
NSSNN	Augmented symplectic integrator on extended phase space $\Phi$4	Nonseparable Hamiltonians
LSNet / SLSNet	Local 2D symplectic Up/Low modules and symmetric reversal	Volume-preserving dynamics
SympFlow	Composition of exact position-only and momentum-only Hamiltonian subflows	Symplectic approximation and flow-map discovery

A canonical example is the point-transformation layer of "Neural Canonical Transformation with Symplectic Flows" (Li et al., 2019): $\Phi$5 Here $\Phi$6 is an invertible coordinate map, implemented with Real-NVP coupling blocks, and the full canonical map is obtained by stacking such layers. Because each layer is volume preserving with Jacobian determinant exactly one, no Jacobian determinants need to be tracked.

A different exact-construction strategy appears in "Deep Neural Networks with Symplectic Preservation Properties" (He et al., 2024). Its symplectomorphism neural network is a composition of three elementary maps: q-shearing,

$\Phi$7

p-shearing,

$\Phi$8

and symplectic stretching,

$\Phi$9

Each block has an explicit inverse and satisfies the symplectic condition exactly.

Integrator-based architectures form another major branch. "Symplectic Methods in Deep Learning" constructs symplectic networks from higher-order explicit separable-Hamiltonian SPRK schemes, including second-, third-, and fourth-order methods with Butcher arrays satisfying the discrete symplectic condition $\bigl(D\Phi(x)\bigr)^{T}J\bigl(D\Phi(x)\bigr)=J,\qquad J=\begin{pmatrix}0&I_n\-I_n&0\end{pmatrix},$0 (Maslovskaya et al., 2024). "Symplectic Neural Flows for Modeling and Discovery" and "Hamiltonian Matching for Symplectic Neural Integrators" build a layer as the composition of exact position-only and momentum-only Hamiltonian subflows generated by scalar time-dependent potentials $\bigl(D\Phi(x)\bigr)^{T}J\bigl(D\Phi(x)\bigr)=J,\qquad J=\begin{pmatrix}0&I_n\-I_n&0\end{pmatrix},$1 and $\bigl(D\Phi(x)\bigr)^{T}J\bigl(D\Phi(x)\bigr)=J,\qquad J=\begin{pmatrix}0&I_n\-I_n&0\end{pmatrix},$2, so that the full network is symplectic by construction and admits backward-error analysis through an associated modified Hamiltonian (Canizares et al., 2024, Canizares et al., 2024).

3. Objectives, estimation strategies, and gradient computation

Training objectives vary with what is assumed known. When the physical Hamiltonian $\bigl(D\Phi(x)\bigr)^{T}J\bigl(D\Phi(x)\bigr)=J,\qquad J=\begin{pmatrix}0&I_n\-I_n&0\end{pmatrix},$3 is available analytically, neural canonical transformation minimizes a variational free energy,

$\bigl(D\Phi(x)\bigr)^{T}J\bigl(D\Phi(x)\bigr)=J,\qquad J=\begin{pmatrix}0&I_n\-I_n&0\end{pmatrix},$4

which is equivalent, up to $\bigl(D\Phi(x)\bigr)^{T}J\bigl(D\Phi(x)\bigr)=J,\qquad J=\begin{pmatrix}0&I_n\-I_n&0\end{pmatrix},$5, to the expectation of $\bigl(D\Phi(x)\bigr)^{T}J\bigl(D\Phi(x)\bigr)=J,\qquad J=\begin{pmatrix}0&I_n\-I_n&0\end{pmatrix},$6 under the model density. When samples from $\bigl(D\Phi(x)\bigr)^{T}J\bigl(D\Phi(x)\bigr)=J,\qquad J=\begin{pmatrix}0&I_n\-I_n&0\end{pmatrix},$7 are available instead, the same framework uses maximum likelihood

$\bigl(D\Phi(x)\bigr)^{T}J\bigl(D\Phi(x)\bigr)=J,\qquad J=\begin{pmatrix}0&I_n\-I_n&0\end{pmatrix},$8

and for separable Hamiltonians augments coordinate-only observations by drawing Gaussian momenta on the fly (Li et al., 2019).

For map regression, SpNN uses a standard mean-squared error objective

$\bigl(D\Phi(x)\bigr)^{T}J\bigl(D\Phi(x)\bigr)=J,\qquad J=\begin{pmatrix}0&I_n\-I_n&0\end{pmatrix},$9

with no penalty terms for deviation from symplecticity or invertibility because both properties are built into the architecture. In the reported experiments, optimization used Adamax with learning rate $\dot x=J\nabla_x H$0, decay factor $\dot x=J\nabla_x H$1 every $\dot x=J\nabla_x H$2 epochs, trained for $\dot x=J\nabla_x H$3 epochs (He et al., 2024).

SympFlow introduces a two-term physics-informed objective: an ODE-residual term enforcing Hamilton’s equations along the learned flow and a Hamiltonian-matching term comparing the network’s exact modified Hamiltonian $\dot x=J\nabla_x H$4 to the target Hamiltonian $\dot x=J\nabla_x H$5. In the supervised variant, it regresses directly on noisy trajectory data with

$\dot x=J\nabla_x H$6

The reported setups use Adam with learning rate $\dot x=J\nabla_x H$7, batch size $\dot x=J\nabla_x H$8, and up to $\dot x=J\nabla_x H$9 epochs (Canizares et al., 2024, Canizares et al., 2024).

The most explicit treatment of gradient computation appears in "Symplectic Neural Networks for learning Generalized Hamiltonians" (Choudhary et al., 25 Jun 2026). There, a parameterized Hamiltonian $(q,p)$0 is integrated with an implicit midpoint or Störmer–Verlet scheme, and the discrete adjoint system coincides with the forward symplectic scheme through cotangent-lift symmetry. The resulting gradient

$(q,p)$1

can be computed by a backward pass solving the same implicit steps in reverse time. The paper further reports a predictor–corrector implicit midpoint step with $(q,p)$2 fixed-point iterations, where $(q,p)$3 typically suffices.

4. Extensions beyond separability and exact canonical maps

A persistent theme in the literature is that symplectic neural flows are not confined to separable Hamiltonians. "Nonseparable Symplectic Neural Networks" addresses general $(q,p)$4 by lifting the dynamics to an extended phase space $(q,p)$5 with augmented Hamiltonian

$(q,p)$6

following Tao (2016). A second-order symmetric splitting into $(q,p)$7, $(q,p)$8, and $(q,p)$9 yields an explicit symplectic integrator. The training loss includes both physical and auxiliary variables, and the reported experiments use $\omega=\sum_i dq_i\wedge dp_i$0, Adam with learning rate $\omega=\sum_i dq_i\wedge dp_i$1, decay $\omega=\sum_i dq_i\wedge dp_i$2 every $\omega=\sum_i dq_i\wedge dp_i$3 epochs, and batch size $\omega=\sum_i dq_i\wedge dp_i$4 (Xiong et al., 2020).

Locally-symplectic architectures extend the design logic from Hamiltonian systems to general divergence-free dynamics. LocSympNets decompose the vector field into adjacent-coordinate Hamiltonian subsystems and implement each with paired Up/Low modules whose Jacobians are $\omega=\sum_i dq_i\wedge dp_i$5 symplectic triangular blocks. SymLocSympNets then enforce time reversibility by constructing a symmetric map $\omega=\sum_i dq_i\wedge dp_i$6 with $\omega=\sum_i dq_i\wedge dp_i$7 (Bajārs, 2021).

Taylor-net occupies an intermediate position. It approximates $\omega=\sum_i dq_i\wedge dp_i$8 with two symmetric Taylor-series subnetworks whose Jacobians are constrained to be symmetric, and then embeds a fourth-order Forest–Ruth or Yoshida-type symplectic integrator into a neural ODE framework (Tong et al., 2020). This is still an intrinsic symplectic design, but the structural prior resides in the Taylor expansion and symmetric factorization of the gradient network rather than in exact elementary symplectomorphisms.

5. Empirical regimes and reported applications

The early canonical-transformation formulation demonstrates both mechanistic interpretability and practical modeling. On a two-dimensional ring potential with Hamiltonian $\omega=\sum_i dq_i\wedge dp_i$9, the variationally trained flow finds a polar transformation automatically: the slow latent coordinate $$\dot q_i=\frac{\partial H}{\partial p_i},\qquad \dot p_i=-\frac{\partial H}{\partial q_i}.$$0 correlates with the polar angle $$\dot q_i=\frac{\partial H}{\partial p_i},\qquad \dot p_i=-\frac{\partial H}{\partial q_i}.$$1, and $$\dot q_i=\frac{\partial H}{\partial p_i},\qquad \dot p_i=-\frac{\partial H}{\partial q_i}.$$2 with the radius $$\dot q_i=\frac{\partial H}{\partial p_i},\qquad \dot p_i=-\frac{\partial H}{\partial q_i}.$$3. On a harmonic chain, the learned frequencies $$\dot q_i=\frac{\partial H}{\partial p_i},\qquad \dot p_i=-\frac{\partial H}{\partial q_i}.$$4 match the known normal-mode dispersion $$\dot q_i=\frac{\partial H}{\partial p_i},\qquad \dot p_i=-\frac{\partial H}{\partial q_i}.$$5, and the Jacobian $$\dot q_i=\frac{\partial H}{\partial p_i},\qquad \dot p_i=-\frac{\partial H}{\partial q_i}.$$6 reproduces the analytic sine-mode shapes for the two slowest modes. On alanine dipeptide, trained by maximum likelihood on $$\dot q_i=\frac{\partial H}{\partial p_i},\qquad \dot p_i=-\frac{\partial H}{\partial q_i}.$$7 trajectories at $$\dot q_i=\frac{\partial H}{\partial p_i},\qquad \dot p_i=-\frac{\partial H}{\partial q_i}.$$8, the two slowest latent coordinates show high mutual information with the backbone dihedrals $$\dot q_i=\frac{\partial H}{\partial p_i},\qquad \dot p_i=-\frac{\partial H}{\partial q_i}.$$9 and $\omega$0, computed via the Kraskov estimator. On MNIST, treating each $\omega$1 image as $\omega$2 coordinates and augmenting with $\omega$3 momenta, a simple MLP classifier fed only the $\omega$4 slow latent coordinates achieves $\omega$5 accuracy for $\omega$6, and conceptual compression using the first $\omega$7–$\omega$8 slow modes yields reasonable reconstructions (Li et al., 2019).

Higher-order explicit symplectic integrators have also been evaluated as deep networks. In 2-D point classification with feature-space augmentation to $\omega$9 and $D\Phi^\top J D\Phi=J$0 layers, reported test accuracies are $D\Phi^\top J D\Phi=J$1 for Hamiltonian-Euler, $D\Phi^\top J D\Phi=J$2 for second-order SPRK, $D\Phi^\top J D\Phi=J$3 for third-order SPRK, and $D\Phi^\top J D\Phi=J$4 for fourth-order SPRK. In learning Kepler dynamics from $D\Phi^\top J D\Phi=J$5 trajectories and $D\Phi^\top J D\Phi=J$6 points, the final $D\Phi^\top J D\Phi=J$7 test error is $D\Phi^\top J D\Phi=J$8 for Hamiltonian-Euler, $D\Phi^\top J D\Phi=J$9 for second-order SPRK, $(q,p)$0 for third-order SPRK, and $(q,p)$1 for fourth-order SPRK (Maslovskaya et al., 2024).

Other architectures focus on exact symplectomorphism learning or long-time prediction. SpNN learns the canonical map $(q,p)$2 with final $(q,p)$3, with errors concentrated near singularities such as $(q,p)$4 and $(q,p)$5-boundaries (He et al., 2024). Taylor-net reports pendulum prediction error $(q,p)$6, compared with $(q,p)$7 for HNN and $(q,p)$8 for ODE-net, and states convergence in $(q,p)$9–$\Phi$00 epochs while HNN and plain ODE-net require $\Phi$01–$\Phi$02 more epochs to reach comparable $\Phi$03 losses (Tong et al., 2020). NSSNN reports long-term stability on Tao’s nonseparable example up to $\Phi$04 with negligible phase drift, robustness to up to $\Phi$05 uniform noise, and scaling to $\Phi$06 vortex particles (Xiong et al., 2020). SympFlow on the harmonic oscillator and Hénon–Heiles yields bounded energy behavior qualitatively similar to classical symplectic methods, while unconstrained MLP baselines drift (Canizares et al., 2024, Canizares et al., 2024). For generalized HNN training with implicit symplectic solvers, representative results include Hénon–Heiles energy drift after $\Phi$07 steps of explicit RK4 at $\Phi$08 versus symplectic NN at $\Phi$09, memory scaling from $\Phi$10 for backprop while the adjoint remains at $\Phi$11, and runtime for $\Phi$12 that is $\Phi$13–$\Phi$14 faster per iteration than backprop (Choudhary et al., 25 Jun 2026).

6. Guarantees, misconceptions, and limitations

A central guarantee is exact or built-in symplecticity. In map-based models this appears as $\Phi$15 up to floating-point rounding and, in several cases, explicit inverses such as $\Phi$16, $\Phi$17, and $\Phi$18 (He et al., 2024). In integrator-based models it appears through discrete symplectic conditions on Butcher arrays, through exact subflow composition, or through symplectic implicit midpoint and Störmer–Verlet schemes (Maslovskaya et al., 2024, Choudhary et al., 25 Jun 2026).

These guarantees should not be conflated with exact preservation of every invariant. One limitation stated explicitly for SympFlow is that exact symplecticity does not guarantee other invariants such as momenta; another is that projection for non-conservative systems breaks symplecticity in the final output. High-frequency time dependence can also increase MLP complexity (Canizares et al., 2024). In LocSympNet experiments, many restarts fail beyond $\Phi$19 in the charged-particle example, illustrating the difficulty of quasi-periodic learning even when short-time accuracy is good (Bajārs, 2021).

Several papers also address common numerical concerns. Symplecticity is linked to non-vanishing gradients in deep architectures: for a symplectic map $\Phi$20, the Jacobian norm satisfies $\Phi$21, and one cited result states that all singular values of $\Phi$22 are $\Phi$23, which prevents vanishing Jacobian products in back-propagation (Maslovskaya et al., 2024). At the same time, exact structure preservation may come with additional computational complexity: implicit symplectic integrators are described as more computationally intensive and make backpropagation through the ODE solver non-trivial, motivating predictor–corrector solvers and discrete-adjoint training (Choudhary et al., 25 Jun 2026).

A final misconception is that symplectic neural flows are synonymous with separable-Hamiltonian models. The available literature directly contradicts that view: NSSNN targets nonseparable Hamiltonians with augmented splitting, generalized HNN training uses implicit midpoint for general $\Phi$24, and LocSympNets target divergence-free systems through local Hamiltonian descriptions rather than a single global separable Hamiltonian (Xiong et al., 2020, Choudhary et al., 25 Jun 2026, Bajārs, 2021). Taken together, these results suggest that the defining feature of symplectic neural flows is not a single parameterization, but the hard incorporation of symplectic or locally symplectic geometry into the learned map, integrator, or latent representation.