A global structure-preserving kernel method for the learning of Poisson systems

Jianyu Hu¹ Juan-Pablo Ortega¹ and Daiying Yin¹

Abstract

A structure-preserving kernel ridge regression method is presented that allows the recovery of globally defined, potentially high-dimensional, and nonlinear Hamiltonian functions on Poisson manifolds out of datasets made of noisy observations of Hamiltonian vector fields. The proposed method is based on finding the solution of a non-standard kernel ridge regression where the observed data is generated as the noisy image by a vector bundle map of the differential of the function that one is trying to estimate. Additionally, it is shown how a suitable regularization solves the intrinsic non-identifiability of the learning problem due to the degeneracy of the Poisson tensor and the presence of Casimir functions. A full error analysis is conducted that provides convergence rates using fixed and adaptive regularization parameters. The good performance of the proposed estimator is illustrated with several numerical experiments.

¹¹footnotetext: Jianyu Hu, Juan-Pablo Ortega, and Daiying Yin are with the Division of Mathematical Sciences, School of Physical and Mathematical Sciences, Nanyang Technological University, Singapore. Their email addresses are Jianyu.Hu@ntu.edu.sg, Juan-Pablo.Ortega@ntu.edu.sg, and YIND0004@e.ntu.edu.sg, respectively.

1 Introduction

The study of Poisson systems plays a crucial role in understanding a wide range of physical phenomena, from classical mechanics to fluid dynamics [Mars 99, Orte 04, Holm 09, Holm 11a, Holm 11b]. These systems, characterized by their inherent geometric structures and conservation laws, present unique challenges in modeling and simulation. Traditional methods for learning and approximating Poisson systems often require extensive domain knowledge and computational resources. In recent years, several machine learning techniques have emerged as powerful tools for addressing complex scientific problems that could offer new avenues for efficiently learning the structures underlying Poisson systems.

Recent literature has extensively explored the learning and prediction of some Poisson systems. These studies incorporate varying degrees of domain-specific knowledge into learning algorithms, enhancing their ability to accurately model and predict such systems’ behavior. An initial approach involved directly learning the Hamiltonian [Grey 19, Davi 23, Hu 24] or Lagrangian [Cran 20, Chen 22] functions using neural networks or other machine learning models, such as Gaussian processes [Offe 24, Pfö 24, Beck 22]. It is worth pointing out that proper variational integrators [Mars 01, Leok 12] should be adopted to match the predicted vector field with observed time series data [Zhu 20]. Researchers have also proposed learning the generating function that produces the flow map [Chen 21, Rath 21]. Alternative approaches aim to parameterize a (possibly universal) family of symplectic or Poisson maps to fit the flow map directly [Chen 19, Jin 20, Baja 23, Jin 22, Eldr 24b], with SympNet being a prominent example [Jin 20]. Using SympNet, symplectic autoencoders have been developed to achieve a structure-preserving model reduction of Hamiltonian systems [Bran 23], effectively mapping high-dimensional systems to low-dimensional representations. Additionally, physics-informed neural networks have become popular for solving known physical partial differential equations with high precision [Rais 19] and convergence analysis in certain scenarios [Doum 23, Shin 20], though they enforce the PDEs as loss functions and thus do not preserve the exact structure. In general, to characterize and replicate the flow of (port-)Hamiltonian systems with exact structure preservation remains challenging, with a notable exception in the literature for linear port-Hamiltonian systems [Orte 24].

This paper proposes a non-standard kernel ridge regression method that allows the recovery of globally defined and potentially high-dimensional and nonlinear Hamiltonian functions on Poisson manifolds out of datasets made of noisy observations of Hamiltonian vector fields. We now spell out three key features of the approach proposed in this paper.

First, much of the existing research in this field operates under the assumption that the underlying dynamical system is Hamiltonian or that the flow map is symplectic [Chen 19, Jin 20], allowing for the use of generating function models or universal parameterizations of symplectic maps to directly model the flow. However, this assumption does not extend to general Poisson systems on manifolds for which the Hamiltonian function cannot be globally expressed in canonical coordinates, and the flow map is generally not symplectic. Moreover, due to the potential degeneracy of the Poisson tensor, Poisson manifolds typically possess a family of functions known as Casimirs, which do not influence the Hamiltonian vector fields. Preserving these Casimirs is an intrinsic requirement for any machine learning algorithm applied to such systems. Inspired by the Darboux-Weinstein theorem [Wein 83, Libe 12, Vais 12], recent approaches have proposed learning coordinate transformations to convert coordinates into canonical ones, followed by inverse transformations to retrieve the dynamics [Jin 22, Baja 23]. While the local existence of these transformations is guaranteed, their global behavior remains unclear, particularly in cases where coordinate patches overlap, and numerical errors on the preservation of Casimir values accumulate. The Lie-Poisson Network [Eldr 24b, Eldr 24a] attempts to approximate flows of the special case of Lie-Poisson systems using compositions of flow maps of linear systems, preserving Casimirs exactly through projection but not addressing the expressiveness of the constructed flows. Rather than learning the flow map, we propose to discover an underlying Hamiltonian function of the Poisson system from Hamiltonian vector fields in a structure-preserving manner. The presence of Casimir functions renders the identification of the Hamiltonian function from the observed vector field infeasible, which is a challenge that lies beyond the scope of the work already presented by the authors in [Hu 24] and necessitates the development of new algorithms to address this issue.

Second, while significant research has focused on the forward problem of learning the flow map of dynamical systems using approaches like Physics-Informed neural networks (PINNs), symbolic regression [Brun 16], Gaussian processes, or kernel methods, the inverse problem of uncovering unknown Hamiltonian functions from data remains relatively underexplored. Unlike neural networks or likelihood-based methods that rely on iterative or stochastic gradient descent to optimize objective functions, kernel methods offer a more straightforward, easy-to-train, and data-efficient alternative. By transforming the optimization problem into a convex one, kernel methods enable the derivation of a closed-form solution for the Hamiltonian function estimation problem in a single step, even on manifolds. Another advantage of kernel methods is that by selecting a kernel with inherent symmetry, the learning scheme can be designed to automatically preserve it, which is particularly important when learning mechanical systems. We recall that due to the celebrated Noether’s Theorem, symmetries usually carry in their wake additional conserved quantities when dealing with mechanical systems. The Reproducing Kernel Hilbert Space (RKHS) framework also facilitates rigorous estimation and approximation error bounds, offering both theoretical insights and interpretability. Additionally, it is well-suited for the design of online learning paradigms.

Finally, while the literature on learning physical dynamical systems is extensive and the idea of structure-preserving kernel regression has already been proposed to recover either interacting potentials of particle systems [Feng 23] or Hamiltonian functions on Euclidean spaces [Hu 24], the majority of studies have been confined to the context of Euclidean spaces, with a lack of attention to scenarios where the underlying configuration space is a more general manifold, which is typically the case in many applications like rigid body mechanics, robotics, fluids, or elasticity. This generalization introduces a two-fold challenge in model construction or for the learning of either the flows of dynamical systems or their corresponding Hamiltonian functions, namely, preserving simlutaneously the manifold and the variational structures. When a model is designed to predict Hamiltonian vector fields on a manifold, ensuring that the predicted flow remains on the manifold and preserves symplectic or Poisson structures necessitates structure-preserving integrators. However, even in Euclidean spaces, the design of variational integrators is complex and often requires implicit methods, indicating a heightened difficulty in the manifold setting. Additionally, standard methods such as Galerkin projections are typically approximately but not exactly manifold-preserving [Boro 20] and typically fail to preserve the variational structure. In [Magg 21], an algorithm is proposed to identify interacting potentials on a Riemannian manifold, though the estimated potential remains dependent only on the Euclidean scalar distance. This approach requires a prior selection of a finite set of basis functions and coordinate charts or embeddings. In contrast, our work introduces a kernel method that globally estimates a Hamiltonian function for Poisson systems on the phase-space manifold and yields a globally well-defined and chart-independent solution. This sets it apart from existing approaches in the literature [Feng 23, Magg 21, Lu 19].

Main Results.

This paper proposes a method to learn a Hamiltonian function $H:P\rightarrow\mathbb{R}$ of a globally defined Hamiltonian system on a $d$ -dimensional Poisson manifold $P$ from noisy realizations of the corresponding Hamiltonian vector field $X_{H}\in\mathfrak{X}(P)$ , that is:

\displaystyle\mathbf{X}_{\sigma^{2}}^{(n)}:=X_{H}(\mathbf{Z}^{(n)})+\bm{\varepsilon}^{(n)},\quad n=1,\ldots,N,

where $\mathbf{Z}^{(n)}$ are IID random variables with values on the Poisson manifold $P$ and with the same distribution $\mu_{\mathbf{Z}}$ , and $\bm{\varepsilon}^{(n)}$ are independent random variables on the tangent spaces $T_{\mathbf{z}^{(n)}}P$ ( $\mathbf{z}^{(n)})$ is a particular realization of $\mathbf{Z}^{(n)}$ with mean zero and variance $\sigma^{2}I_{d}$ . In this setup, we consider the inverse problem consisting in recovering the Hamiltonian function $H$ , or a dynamically equivalent version of it, by solving the following optimization problem,

\displaystyle\widehat{h}_{\lambda,N}:=\mathop{\arg\min}\limits_{h\in\mathcal{H}_{K}}\frac{1}{N}\sum_{n=1}^{N}\left\|X_{h}(\mathbf{Z}^{(n)})-\mathbf{X}^{(n)}_{\sigma^{2}}\right\|_{g}^{2}+\lambda\|h\|_{\mathcal{H}_{K}}^{2},

(1.1)

where $X_{h}$ is the Hamiltonian vector field corresponding to the function $h\in{\mathcal{H}}_{K}$ , ${\mathcal{H}}_{K}$ is the reproducing kernel Hilbert space (RKHS) associated to a kernel function $K:P\times P\rightarrow\mathbb{R}$ that will be introduced later on, $g$ is a Riemannian metric on $P$ , and $\lambda\geq 0$ is a Tikhonov regularization parameter. The solution $\widehat{h}_{\lambda,N}$ in the RKHS is sought to minimize the natural regularized empirical risk defined in this setup. We shall call the solution $\widehat{h}_{\lambda,N}$ of the optimization problem (1.1) the structure-preserving kernel estimator of the Hamiltonian function.

We now summarize the outline and the main contributions of the paper.

1.

In Section 2, we review basic concepts in Riemannian geometry, including differentials on Riemannian manifolds and the definitions of random variables and integration in this context. Additionally, we provide the necessary background on Poisson mechanics, covering topics such as Hamiltonian vector fields, Casimir functions, compatible structures, and Lie-Poisson systems, illustrated with two real-world examples: rigid-body dynamics and underwater vehicle dynamics.
2.

In Section 3, we extend the differentiable reproducing property, previously established for the compact [Zhou 08] and Euclidean cases [Hu 24] to Riemannian manifolds. Specifically, we show that the differential of a function $f\in\mathcal{H}_{K}$ in the reproducing kernel Hilbert space (RKHS) $\mathcal{H}_{K}$ associated to a kernel function $K:P\times P\rightarrow\mathbb{R}$ , when paired with a vector $v\in T_{x}P$ , can be represented as the RKHS inner product of $f\in\mathcal{H}_{K}$ with the differential of the kernel section of $K$ evaluated along $v\in T_{x}P$ . This property is pivotal throughout our analysis. Following this, we derive the solution of the optimization problem (1.1) using operator representations and ultimately express it as a closed-form formula made of a linear combination of kernel sections, similar to the standard Representer theorem. Importantly, this formula provides a globally defined expression on the Poisson manifold $P$ , independent of the choice of local coordinates, though we also present the formula in local charts. In Subsection 3.5, we study the case in which the Hamiltonian that needs to be learned is known to be a priori invariant with respect to the canonical action of a given Lie group and show that if the kernel is chosen to be invariant under that group action, then the flow of the estimated Hamiltonian inherently conserves momentum maps, as a direct consequence of Noether’s theorem.
3.

In Section 4, we utilize the operator representation of the estimator to derive bounds for estimation and approximation errors. These bounds demonstrate that, as the number of random samples increases, the estimator converges to a true Hamiltonian function with high probability.
4.

In Section 5, we implement the structure-preserving kernel ridge regression in Python and conduct extensive numerical experiments. These experiments focus on two primary scenarios: Hamiltonian systems on symplectic manifolds (Subsection 5.2) and Lie-Poisson systems (Subsection 5.1). In subsection 5.2, we aim to recover the Hamiltonian function for the classical two-vortex system on products of $2$ -spheres, which exhibits singularities, as well as another well-behaved Hamiltonian function represented as a 3-norm defined on the same manifold. In Subsection 5.1, we address symmetric Poisson systems on cotangent bundles of Lie groups. Through Lie-Poisson reduction, these systems can be modeled as Lie-Poisson systems on the dual of a Lie algebra. Numerical simulations validate the effectiveness, ease of training, and data efficiency of our proposed kernel-based approach.

We emphasize that, although this paper focuses on recovering Hamiltonian functions from observed Hamiltonian vector fields, the proposed framework is more broadly applicable. Specifically, it can be extended to any scenario where the observed quantity (in this paper $J\nabla H$ ) is the image by a vector bundle map of the differential of a function (in this paper $H$ ) on the manifold. This makes the framework versatile for various systems beyond Hamiltonian dynamics; for example, in gradient systems with an underlying control-dependent potential function, the same algorithm can be applied to learn an optimal control force based on sensor-measured external forces.

2 Preliminaries on Poisson mechanics

2.1 Hamiltonian systems on Poisson manifolds

The aim of this paper is to learn Hamiltonian systems on Poisson manifolds. In this subsection, we will provide a brief introduction to Poisson mechanics relevant to the learning problems; for more details, see [Mars 99, Orte 04, Libe 12, Vais 12]. Let $P$ be a manifold and denote by $C^{\infty}(P)$ the space of all smooth functions defined on it.

Definition 2.1.

A Poisson bracket (or a Poisson structure) on a manifold $P$ is a bilinear map $\{\cdot,\cdot\}:C^{\infty}(P)\times C^{\infty}(P)\to C^{\infty}(P)$ such that:

(i)

$\left(C^{\infty}(P),\{\cdot,\cdot\}\right)$ is a Lie algebra.

(ii)

$\{\cdot,\cdot\}$ is a derivation in each factor, that is,

\{FG,H\}=\{F,H\}G+F\{G,H\},

for all $F,G$ , and $H\in C^{\infty}(P)$ .

A manifold $P$ endowed with a Poisson bracket $\{\cdot,\cdot\}$ on $C^{\infty}(P)$ is called a Poisson manifold.

Example 2.2 (Symplectic bracket).

Any symplectic manifold $(P,\omega)$ is a Poisson manifold with the Poisson bracket defined by the symplectic form $\omega$ as follows:

\displaystyle\{F,G\}(z)=\omega(X_{F}(z),X_{G}(z)),\quad\mbox{$F,G\in C^{\infty}(P)$, }

(2.1)

where $X_{F},X_{G}$ are the Hamiltonian vector fields of $F,G\in C^{\infty}(P)$ , respectively, uniquely determined by the equalities $\mathbf{i}_{X_{F}}\omega=\mathbf{d}F$ and $\mathbf{i}_{X_{G}}\omega=\mathbf{d}G$ ( $\mathbf{i}$ denotes the interior product and $\mathbf{d}$ the exterior derivative).

Example 2.3 (Lie-Poisson bracket).

If $\mathfrak{g}$ is a Lie algebra, then its dual algebra $\mathfrak{g}^{*}$ is a Poisson manifold with respect to the Lie-Poisson brackets $\{\cdot,\cdot\}_{+}$ and $\{\cdot,\cdot\}_{-}$ defined by

\displaystyle\{F,G\}_{\pm}(\mu)=\pm\left\langle\mu,\left[\frac{\delta F}{\delta\mu},\frac{\delta G}{\delta\mu}\right]\right\rangle,\quad\mbox{$F,G\in C^{\infty}(\mathfrak{g}^{\ast})$,}

(2.2)

where $[\cdot,\cdot]$ is the Lie bracket on $\mathfrak{g}$ , $\langle\cdot,\cdot\rangle$ denotes the pairing between $\mathfrak{g}^{*}$ and $\mathfrak{g}$ , and $\frac{\delta F}{\delta\mu}\in\mathfrak{g}$ is the functional derivative of $F$ at $\mu\in\mathfrak{g}^{*}$ defined as the unique element that satisfies

\displaystyle\lim_{\varepsilon\to 0}\frac{1}{\varepsilon}\left[F(\mu+\varepsilon\delta\mu)-F(\mu)\right]=\left\langle\delta\mu,\frac{\delta F}{\delta\mu}\right\rangle\quad\mbox{for any $\delta\mu\in\mathfrak{g}^{\ast}$}.

(2.3)

By the derivation property of the Poisson bracket, the value of the bracket $\{F,G\}$ at $z\in P$ (and thus $X_{F}(z)$ as well) depends on $F$ only through the differential $\mathbf{d}F(z)$ . Thus, a Poisson tensor associated with the Poisson bracket can be defined as follows.

Definition 2.4.

Let $(P,\left\{\cdot,\cdot\right\})$ be a Poisson manifold. The Poisson tensor of this Poisson bracket is the contravariant anti-symmetric two-tensor

B:T^{*}P\times T^{*}P\rightarrow\mathbb{R},

defined by

\displaystyle B(z)\left(\alpha_{z},\beta_{z}\right)=\{F,G\}(z),\quad\mbox{for any $\alpha_{z},\beta_{z}\in T_{z}^{\ast}P$, $z\in P$},

(2.4)

and where $F,G\in C^{\infty}(P)$ are any two functions such that $\mathbf{d}F(z)=\alpha_{z}$ and $\mathbf{d}G(z)=\beta_{z}\in T_{z}^{*}P$ .

Hamiltonian vector fields and Hamilton’s equations. Given a Poisson manifold $(P,\left\{\cdot,\cdot\right\})$ and $H\in C^{\infty}(P)$ , there is a unique vector field $X_{H}\in\mathfrak{X}(P)$ such that

\displaystyle X_{H}[F]=\{F,H\},\quad\mbox{for all $F\in C^{\infty}(P)$}.

(2.5)

The unique vector field $X_{H}$ is called the Hamiltonian vector field associated with the Hamiltonian function $H$ . Let $F_{t}:P\rightarrow P$ be the flow of the Hamiltonian vector field $X_{H}$ ; then, for any $G\in C^{\infty}(P)$ ,

\displaystyle\frac{d}{dt}(G\circ F_{t})=\{G,H\}\circ F_{t}=\{G\circ F_{t},H\},\quad\text{or}\quad\dot{G}=\{G,H\}.

(2.6)

We call (2.6) the Hamilton equations associated with the Hamiltonian function $H$ . Let $B^{\sharp}:T^{*}P\to TP$ be the vector bundle map associated to the Poisson tensor $B$ defined in (2.4), that is,

\displaystyle B(z)\left(\alpha_{z},\beta_{z}\right)=\alpha_{z}\cdot B^{\sharp}(z)(\beta_{z}).

(2.7)

Then the Hamiltonian vector field defined in (2.5) can be expressed as $X_{H}=B^{\sharp}(\bm{\mathrm{d}}H)$ . Furthermore, the Poisson bracket is given by

\displaystyle\{F,H\}(z)=B(z)(\mathbf{d}F(z),\mathbf{d}H(z))=\mathbf{d}F(z)\cdot B^{\sharp}(z)(\mathbf{d}H(z)).

Poisson geometry is closely related to symplectic geometry: for instance, every Poisson bracket determines a (generalized) foliation of a Poisson manifold into symplectic (immersed) submanifolds such that the leaf that passes through a given point has the dimension of the rank of the Poisson tensor at that point. We recall that the non-degeneracy hypothesis on symplectic forms implies that the Poisson tensor, in that case, has constant maximal rank. The following paragraph describes a distinctive feature of Poisson manifolds with respect to the symplectic case.

Casimir functions. As we just explained, an important difference between the symplectic and Poisson cases is that the Poisson tensor could be degenerate, which leads to the emergence of Casimir functions. A function $C\in C^{\infty}(P)$ is called a Casimir function if $\{C,F\}=0$ for all $F\in C^{\infty}(P)$ , that is, $C$ is constant along the flow of all Hamiltonian vector fields, or, equivalently, $X_{C}=0$ , that is, $C$ generates trivial dynamics. Casimir functions are conserved quantities along the equations of motion associated with any Hamiltonian function. The Casimirs of a Poisson manifold are the center $\mathcal{C}(P)$ of the Poisson algebra $(C^{\infty}(P),\left\{\cdot,\cdot\right\})$ which is defined by

\displaystyle\mathcal{C}(P)=\{C\in C^{\infty}(P)\mid\{C,F\}=0,\text{ for all }F\in C^{\infty}(P)\}.

If the Poisson tensor is non-degenerate (symplectic case), then $\mathcal{C}(P)$ consists only of all constant functions on $P$ . The following Lie-Poisson systems are two examples that admit non-constant Casimirs.

Example 2.5 (Rigid body Casimirs [Mars 99]).

The motion of a rigid body follows a Lie-Poisson equation on the dual of the Lie algebra $\mathfrak{so}(3)^{*}$ of the special orthogonal group $SO(3)$ determined by the Hamiltonian function

\displaystyle H(\Pi)=\frac{1}{2}\Pi^{\top}\mathbb{I}^{-1}\Pi,

where $\mathbb{I}=\operatorname{diag}\{I_{1},I_{2},I_{3}\}$ is a diagonal matrix. The Poisson bracket on $\mathfrak{so}(3)^{*}\simeq\mathbb{R}^{3}$ is given by

\displaystyle\{F,K\}(\Pi)=-\Pi\cdot(\nabla F\times\nabla K).

(2.8)

Hence, it can be easily checked that $C(\Pi)=\|\Pi\|^{2}$ is a Casimir function of Poisson structure (2.8). Furthermore, for any analytic function $\Phi:\mathbb{R}\to\mathbb{R}$ , the function $\Phi\circ C$ is also a Casimir.

Example 2.6 (Underwater vehicle Casimirs [Leon 97]).

The underwater vehicle dynamics can be described using the Lie-Poisson structure on $\mathfrak{so}(3)^{*}\times{\mathbb{R}^{3}}^{*}\times{\mathbb{R}^{3}}^{*}$ and the Hamiltonian function

\displaystyle H(\Pi,\mathrm{Q},\Gamma)=\frac{1}{2}\left(\Pi^{\top}A\Pi+2\Pi^{\top}B^{\top}\mathrm{Q}+\mathrm{Q}^{\top}C\mathrm{Q}-2mg\left(\Gamma\cdot r_{\mathrm{G}}\right)\right),

where $A,B,C$ are matrices coming from first principles and physical laws, $m$ is the mass of the body, $g$ is the gravity constant, and $r_{G}$ is the vector that links the center of buoyancy to the center of gravity. For the underwater vehicle, $\Pi$ and $\mathrm{Q}$ correspond, respectively, to angular and linear components of the impulse of the system (roughly, the momentum of the system less a term at infinity). The vector $\Gamma$ describes the direction of gravity in body-fixed coordinates. The Poisson bracket on $\mathfrak{so}(3)^{*}\times{\mathbb{R}^{3}}^{*}\times{\mathbb{R}^{3}}^{*}$ is

\{F,K\}(\Pi,\mathrm{Q},\Gamma)=\nabla F^{\top}\Lambda(\Pi,\mathrm{Q},\Gamma)\nabla K,

where $F$ and $K$ are smooth functions on $\mathfrak{so}(3)^{*}\times{\mathbb{R}^{3}}^{*}\times{\mathbb{R}^{3}}^{*}$ and the Poisson tensor $\Lambda$ is given by

\Lambda(\Pi,\mathrm{Q},\Gamma)=\left(\begin{array}[]{ccc}\hat{\Pi}&\hat{\mathrm{Q}}&\hat{\Gamma}\\ \hat{\mathrm{Q}}&0&0\\ \hat{\Gamma}&0&0\end{array}\right).

In this equation, the operation ${}^{\wedge}:\mathbb{R}^{3}\rightarrow\mathfrak{so}(3)$ is the standard isomorphism of $\mathbb{R}^{3}$ with the Lie algebra of the rotation group $SO(3)$ and is defined by $\hat{v}w=v\times w$ for $v,w\in\mathbb{R}^{3}$ . There are six Casimir functions for the underwater vehicle dynamics

	$\displaystyle C_{1}(\Pi,\mathrm{Q},\Gamma)=\mathrm{Q}\cdot\Gamma,\quad C_{2}(\Pi,\mathrm{Q},\Gamma)=\\|\mathrm{Q}\\|^{2},\quad C_{3}(\Pi,\mathrm{Q},\Gamma)=\\|\Gamma\\|^{2},$
	$\displaystyle C_{4}(\Pi,\mathrm{Q},\Gamma)=\Pi\cdot\mathrm{Q},\quad C_{5}(\Pi,\mathrm{Q},\Gamma)=\Pi\cdot\Gamma,\quad C_{6}(\Pi,\mathrm{Q},\Gamma)=\\|\Pi\\|^{2}.$

It can be easily checked that for any analytic function $\Phi:\mathbb{R}^{6}\to\mathbb{R}$ , the function $\Phi\circ(C_{1},\ldots,C_{6})$ is also a Casimir.

The momentum map and conservation laws. Symmetries are particularly important in mechanics since they are, on many occasions, associated with the appearance of conservation laws. This fact is the celebrated Noether’s Theorem [Noet 18]. A modern and elegant way to encode this fact uses the momentum map, an object introduced in [Lie 93, Sour 65, Sour 69, Kost 09, Sour 66]. The definition of the momentum map requires a canonical Lie group action, that is, an action that respects the Poisson bracket. Its existence is guaranteed when the infinitesimal generators of this action are Hamiltonian vector fields. In other words, if the Lie algebra $\mathfrak{g}$ of the group $G$ that acts canonically on the Poisson manifold $(P,\{\cdot,\cdot\})$ then, for each $\xi\in\mathfrak{g}$ with associated infinitesimal generator vector field $\xi_{P}\in\mathfrak{X}(P)$ , we require the existence of a globally defined function $\mathbf{J}^{\xi}\in C^{\infty}(P)$ , such that

\xi_{M}=X_{\mathbf{J}^{\xi}}.

Definition 2.7.

Let $\mathfrak{g}$ be the Lie algebra of a Lie group $G$ acting canonically on the Poisson manifold $(P,\,\{\cdot,\cdot\})$ . Suppose that for any $\xi\in\mathfrak{g}$ , the vector field $\xi_{P}$ is globally Hamiltonian, with Hamiltonian function $\mathbf{J}^{\xi}\in C^{\infty}(P)$ . The map $\mathbf{J}:P\rightarrow\mathfrak{g}^{\ast}$ defined by the relation

\langle\mathbf{J}(z),\,\xi\rangle=\mathbf{J}^{\xi}(z),

for all $\xi\in\mathfrak{g}$ and $z\in M$ , is called a momentum map of the $G$ –action.

Notice that momentum maps are not uniquely determined; indeed, $\mathbf{J}_{1}$ and $\mathbf{J}_{2}$ are momentum maps for the same canonical action if and only if for any $\xi\in\mathfrak{g}$ , $\mathbf{J}_{1}^{\xi}-\mathbf{J}_{2}^{\xi}\in\mathcal{C}(P)$ . Obviously, if $P$ is symplectic and connected, then $\mathbf{J}$ is determined up to a constant in $\mathfrak{g}^{\ast}$ .

Noether’s theorem is formulated in terms of the momentum map by saying that the fibers of the momentum map are preserved by the flow of the Hamiltonian vector field associated with any $G$ -invariant Hamiltonian function.

Example 2.8 (The linear momentum).

Take the phase space of the $N$ –particle system, that is, $T^{\ast}\mathbb{R}^{3N}$ . The additive group $\mathbb{R}^{3}$ acts on it by spatial translation on each factor: $\bm{v}\cdot(\bm{q}_{i},\,\bm{p}^{i})=(\bm{q}_{i}+\bm{v},\,\bm{p}^{i})$ , with $i=1,\ldots,\,N$ . This action is canonical and has an associated momentum map that coincides with the classical linear momentum

\begin{array}[]{cccc}\mathbf{J}:&T^{\ast}\mathbb{R}^{3N}&\longrightarrow&{\rm Lie}(\mathbb{R}^{3})\simeq\mathbb{R}^{3}\\ &(\bm{q}_{i},\,\bm{p}^{i})&\longmapsto&\sum_{i=1}^{N}\bm{p}_{i}.\end{array}

Example 2.9 (The angular momentum).

Let ${\rm SO}(3)$ act on $\mathbb{R}^{3}$ and then, by lift, on $T^{\ast}\mathbb{R}^{3}$ , that is, $A\cdot(\bm{q},\,\bm{p})=(A\bm{q},\,A\bm{p})$ . This action is canonical and has an associated momentum map

\begin{array}[]{cccc}\mathbf{J}:&T^{\ast}\mathbb{R}^{3}&\longrightarrow&\mathfrak{so(3)}^{\ast}\simeq\mathbb{R}^{3}\\ &(\bm{q},\,\bm{p})&\longmapsto&\bm{q}\times\bm{p}.\end{array}

which is the classical angular momentum.

Example 2.10 (Lifted actions on cotangent bundles).

The previous two examples are particular cases of the following situation. Let $G$ be a Lie group acting on the manifold $Q$ and then by lift on its cotangent bundle $T^{\ast}Q$ . Any such lifted action is canonical with respect to the canonical symplectic form on $T^{\ast}Q$ and has an associated momentum map $\mathbf{J}:T^{\ast}Q\rightarrow\mathfrak{g}^{\ast}$ given by

\langle\mathbf{J}(\alpha_{q}),\,\xi\rangle=\langle\alpha_{q},\,\xi_{Q}(q)\rangle,

(2.9)

for any $\alpha_{q}\in T^{\ast}Q$ and any $\xi\in\mathfrak{g}$ .

2.2 High-order differentials and random variables on Riemannian manifolds

We start by reviewing higher-order differentials and the spaces of bounded differentiable functions on general manifolds, which will allow us to examine the regularity of functions in a reproducing kernel Hilbert space (RKHS) later on in Section 3.1. Finally, we introduce some basic facts about statistics on Riemannian manifolds that will be needed later on. For further details, see [Penn 19, Hsu 02, Emer 12] and references therein.

Let $M$ be a smooth manifold with an atlas $\{U_{\beta},\varphi_{\beta},\beta\in\mathcal{I}\}$ , where $\mathcal{I}$ is an index set. A function $f:M\to\mathbb{R}$ is said to be in the $C^{k}$ class if for any chart $\varphi_{\beta},\beta\in\mathcal{I}$ , the function $f\circ\varphi_{\beta}:\varphi_{\beta}(U_{\beta})\to\mathbb{R}$ is a $C^{k}$ function and the $k$ -order derivative of $f\circ\varphi_{\beta}$ is continuous.

Definition 2.11 (Differentials).

If $f:M\to\mathbb{R}$ is in the $C^{1}$ class, we define its differential $Df:M\to T^{*}M$ of $f$ as

\displaystyle Df(x)\cdot v:=\left.\frac{d}{dt}\right|_{t=0}f(c(t)),\quad\mbox{for all $x\in M,v\in T_{x}M$},

where $c:[-1,1]\to M$ is a differentiable curve with $c(0)=x$ and $c^{\prime}(0)=v$ .

In the sequel, we shall also denote $\mathbf{d}f$ as the differential of $f$ . The differential $Df$ can be identified with a function on $TM$ . Then, if $f$ is regular enough, we can continue taking the differential of $Df$ on $TM$ . In this way, we can define high-order differentials as follows.

Definition 2.12 (High-order differentials).

If $f:M\to\mathbb{R}$ is in $C^{k}$ class, we define the $k$ -order differential of $f$ denoted as $D^{k}f:T^{k}M\to\mathbb{R}$ inductively to be the differential of $D^{k-1}f:T^{k-1}M\to\mathbb{R}$ .

Remark 2.13.

For a differentiable function $K:M\times M\to\mathbb{R}$ , we denote $D^{(m,l)}K:T^{m}M\times T^{l}M\to\mathbb{R}$ as the $m$ -order differential with respect to the first variable and $l$ -order differential with respect to the second variable.

Remark 2.14 (Functions with bounded differentials).

Having defined high-order differentials, we can now define the spaces $C_{b}^{s}(M)$ of functions with bounded differentials for all positive integers $s$ . We assume first that $M$ is a second-countable manifold, and we hence can construct a complete Riemannian metric $g$ on it [Nomi 61]. A general result [Tu 08, Proposition 12.3] guarantees that the higher order tangent $T^{k}M$ bundles are also second countable and can hence also be endowed with complete metrics $g_{k}$ Let now $C^{s}_{b}(M)$ be the subset of $C^{s}(M)$ given by

\displaystyle C_{b}^{s}(M):=\left\{f\in C^{s}(M)~{}\mid~{}\|f\|_{C_{b}^{s}}:=\|f\|_{\infty}+\sum_{k=1}^{s}\|D^{k}f\|_{\infty}<\infty\right\},

where $\|f\|_{\infty}:=\sup_{x\in M}|f(x)|$ and

\displaystyle\|D^{k}f\|_{\infty}:=\sup_{y\in T^{k-1}M}\ \sup_{v\in T_{y}(T^{k-1}M),v\neq 0}\frac{\left|D^{k}f(y)\cdot v\right|}{\|v\|_{k-1}}.

Here, $\|v\|_{k-1}$ stands for the norm of $v$ in the tangent space $T_{y}T^{k-1}M$ that is defined using the complete Riemannian metric $g_{k-1}$ on $T^{k-1}M$ .

The metric $g=g_{0}$ is also used to define the gradient $\nabla f\in\mathfrak{X}(M)$ of a function $f:M\to\mathbb{R}$ which is the vector field defined by its differential as

\displaystyle g(x)(\nabla f(x),v)=Df(x)\cdot v,\quad\text{for\ all\ }v\in T_{x}M.

For a differentiable function $f:M\to\mathbb{R}$ , a Fundamental Theorem of Calculus can be formulated as

\displaystyle f(x)-f(y)

\displaystyle=\int_{a}^{b}\frac{d}{dt}f(\gamma(t))\mathrm{d}t=\int_{a}^{b}g(\gamma(t))\left(\nabla f(\gamma(t)),\gamma^{\prime}(t)\right)\mathrm{d}t=\int_{a}^{b}Df(\gamma(t))\cdot\gamma^{\prime}(t)\mathrm{d}t,

(2.10)

where $\gamma:[a,b]\to M$ is a differentiable curve with $\gamma(a)=x$ and $\gamma(b)=y$ . Taking absolute values and using the Cauchy-Schwarz inequality, we obtain

\displaystyle|f(x)-f(y)|\leq\|Df\|_{\infty}\int_{a}^{b}\|\gamma^{\prime}(t)\|_{g}\mathrm{d}t=\|Df\|_{\infty}\ \text{length}(\gamma),

for each $f\in C^{1}_{b}(M)$ . Note that, if we choose $\gamma$ to be the geodesic connecting $x$ and $y$ , then $\text{length}(\gamma)$ is the geodesic distance of $x$ and $y$ .

Random variables on Riemannian manifolds. Let $(M,g)$ be a smooth, connected, and complete $d$ -dimensional Riemannian manifold with Riemannian metric $g$ . Existence and uniqueness theorems from the theory of ordinary differential equations guarantee that there exists a unique geodesic integral curve $\gamma_{x,v}$ going through the point $x\in M$ at time $t=0$ with tangent vector $v\in T_{x}M$ . The Riemannian exponential map at $x\in M$ denoted by $\operatorname{exp}_{x}$ maps each vector $v\in T_{x}M$ to the point of the manifold reached in a unit time, that is,

\displaystyle\operatorname{exp}_{x}(v):=\gamma_{x,v}(1),

where $\gamma_{x,v}$ is the geodesic starting at $x$ along the vector $v$ . Let $t_{0}(x,v)>0$ be the first time when the geodesic $t\mapsto\operatorname{exp}_{x}(tv)$ stops to be length-minimizing. The point $\operatorname{exp}_{x}(t_{0}(x,v)v)$ is called a cut point and the corresponding tangent vector $t_{0}(x,v)v$ is called a tangent cut point. The set of all cut points of all geodesics starting from $x\in M$ is the cut locus $C(x)\in M$ and the set of corresponding vectors the tangential cut locus $\mathcal{C}(x)\in T_{x}M$ . Thus, we have $C(x)=\operatorname{exp}_{x}(\mathcal{C}(x))$ . Define the maximal definition domain for the exponential chart as

\displaystyle\mathcal{D}(x):=\left\{tv\in T_{x}M\mid v\in S^{d-1},0\leq t<t_{0}(x,v)\right\},

where $S^{d-1}$ is the unit sphere in the tangent space with respect to the metric $g(x)$ . Then the exponential map $\operatorname{exp}_{x}$ is a diffeomorphism from $\mathcal{D}(x)$ on its image $\operatorname{exp}_{x}(\mathcal{D}(x))$ . Denote by $\operatorname{log}_{x}$ the inverse of the exponential map. The Riemannian distance $d(x,y)=\|\operatorname{log}_{x}(y)\|_{g(x)}$ . Recall that for connected complete Riemannian manifolds we have that

\displaystyle M=\operatorname{exp}_{x}(\mathcal{D}(x))\cup\mathcal{C}(x).

It is easy to check that $\mathcal{D}(x)$ is an open star-shaped subset of $T_{x}M$ , so $\left(\mathcal{D}(x),\exp_{x}\right)$ is a smooth chart and one has that for any Borel measurable function $f:M\rightarrow\mathbb{R}$ on $M$

\int_{M}fdV_{M}=\int_{\mathcal{D}(x)}\left(f\circ\exp_{x}\right)\sqrt{\operatorname{det}\left(G_{x}\right)}d\lambda^{d},

where $dV_{M}$ is the Riemannian measure on the manifold $(M,g),G_{x}:=\exp_{x}^{*}g$ is the pullback of $g$ by $\exp_{x}$ , and $\lambda^{d}$ denotes the $d$ -dimensional Lebesgue measure on $\mathbb{R}^{d}$ .

Definition 2.15 (Random variables on Riemannian manifold).

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space and $(M,g)$ be a Riemannian manifold. A random variable $\xi$ on the Riemannian manifold $M$ is a Borel measurable function from $\Omega$ to $M$ .

Let $\mathcal{B}(M)$ be the Borel $\sigma$ -algebra of $M$ . The random variable $\xi$ on the Riemannian manifold $(M,g)$ has a probability density function $p:M\rightarrow\mathbb{R}$ (real, positive and integrable function) if for all $A\in\mathcal{B}(M)$ ,

\mathbb{P}(X\in A)=\int_{A}p(x)dV_{M}(x)\quad\text{ and }\quad\mathbb{P}(M)=\int_{M}p(x)dV_{M}(x)=1.

The definitions of the mean and variance of a random variable on a Riemannian manifold can be found in the literature [Penn 06].

Remark 2.16 (Random variables on tangent spaces).

Since the tangent space of a manifold at any given point is a vector space, the definition of a random variable on them is the same as for vector spaces. Let $\eta$ be a random variable on the tangent space $T_{x}M$ with density $p:T_{x}M\rightarrow\mathbb{R}$ , that is, $\eta$ is a Borel measurable function from $\Omega$ to $T_{x}M$ .

(i) (Mean and variance). The mean and variance of a random variable $V$ on $T_{x}M$ are defined as

	$\displaystyle\mathbb{E}[\eta]:$	$\displaystyle=\int_{T_{x}M}vp(v)d\lambda^{d}(v)\in T_{x}M,$
	$\displaystyle\operatorname{Var}[\eta]:$	$\displaystyle=\mathbb{E}\left[\\|\eta-\mathbb{E}[\eta]\\|_{g}^{2}\right]=\int_{T_{x}M}\left\\|v-\mathbb{E}[\eta]\right\\|_{g}^{2}p(v)d\lambda^{d}(v),$

where $\lambda^{d}$ is the $d$ -dimensional Lebesgue measure.

(ii) (Independence). Let $\eta_{x}$ and $\eta_{y}$ are two random variables on the tangent spaces $T_{x}M$ and $T_{y}M$ , respectively. Let $\mathcal{B}(T_{x}M)$ and $\mathcal{B}(T_{y}M)$ be the Borel $\sigma$ -algebras of $T_{x}M$ and $T_{y}M$ , respectively. Then we say $\eta_{x}$ and $\eta_{y}$ are independent if for all $A_{x}\in\mathcal{B}(T_{x}M)$ and $A_{y}\in\mathcal{B}(T_{y}M)$ , we have that

\displaystyle\mathbb{P}(\eta_{x}\in A_{x},\eta_{y}\in A_{y})=\mathbb{P}(\eta_{x}\in A_{x})\mathbb{P}(\eta_{y}\in A_{y}).

2.3 Compatible structures and local coordinate representations

This subsection introduces the concept of compatible structure, a vector bundle map from the tangent space of the Poisson manifold to itself that enables us to express Hamiltonian vector fields as images of gradients by a vector bundle map. Using this tool, we can formulate adapted differential reproducing properties in RKHS setups and establish an operator framework later in Section 3 for them. Finally, we shall apply the Poisson-Darboux-Lie-Weinstein theorem to provide local coordinate representations of the compatible structure.

We first endow the Poisson manifold $(P,\left\{\cdot,\cdot\right\})$ with a Riemannian metric $g$ and define the vector bundle map $J:TP\to TP$ by

\displaystyle J(z)v:=B^{\sharp}(z)\left(g^{\flat}(z)(v)\right),\quad\mbox{for all $z\in P,v\in T_{z}P$,}

(2.11)

where $B^{\sharp}:T^{*}P\to TP$ is the vector bundle map defined in (2.7) and $g^{\flat}:TP\to T^{*}P$ is the vector bundle isomorphism determined by the Riemannian metric $g$ . Notice that the vector bundle map $B^{\sharp}$ may not be an isomorphism due to the possible degeneracy of the Poisson tensor $B$ , but nevertheless, the map $J$ is always well-defined since we do not need to use the inverse of $B^{\sharp}$ at any moment. Note that the Poisson tensor $B$ and the Riemannian metric $g$ are linked by the following relation:

\displaystyle B(z)(\alpha_{z},\beta_{z})=\alpha_{z}\cdot B^{\sharp}(z)(\beta_{z})=g(z)\left(g^{\sharp}(z)(\alpha_{z}),B^{\sharp}(z)(\beta_{z})\right),\quad\mbox{for all $\alpha_{z},\beta_{z}\in T^{*}_{z}P$.}

(2.12)

We call $J:TP\to TP$ the compatible structure associated with the Poisson tensor $B$ and the Riemannian metric $g$ . Using the compatible structure $J$ , the Hamiltonian vector field associated with a Hamiltonian function $h$ can be expressed as

\displaystyle X_{h}=B^{\sharp}(\mathbf{d}h)=B^{\sharp}(g^{\flat}(\nabla h))=J\nabla h.

(2.13)

Remark 2.17.

If the Poisson manifold $P$ is symplectic with a symplectic form $\omega$ , then $B^{\sharp}=\omega^{\sharp}$ , and the relationship (2.12) reduces to

\displaystyle g(z)(u,v)=\omega(z)(u,-Jv),\quad u,v\in T_{z}P,

where $J(z)v:=\omega^{\sharp}(z)(g^{\flat}(z)(v))$ is, in this case, an isomorphism from $TP$ to $TP$ . When the compatible structure $J$ defined in (2.12) satisfies $J^{2}=-I$ , it is called a complex structure, and the corresponding manifold is said to be a Kähler manifold. Even-dimensional Euclidean vector spaces are special cases of Kähler manifolds by choosing $J$ as the canonical (Darboux) symplectic matrix $J_{\mathrm{can}}$ , whose Hamiltonian vector field associated with a Hamiltonian function is expressed as $X_{h}=J_{\mathrm{can}}\nabla h$ .

Using the compatible structure $J$ , we obtain a useful property of Casimir functions which will be needed later on in Section 3 to prove the uniqueness of the solution to the learning problem.

Proposition 2.18.

Let $(P,\left\{\cdot,\cdot\right\})$ be a Poisson manifold equipped with a Riemannian manifold $g$ . Then for any vector field $Y\in\mathfrak{X}(P)$ and any $h\in C^{\infty}(P)$ :

\displaystyle g(z)(J(z)\nabla h(z),Y(z))=g(z)(\nabla h(z),-J(z)Y(z)),\quad\mbox{for all $z\in P$}.

(2.14)

Furthermore, Casimir functions are constant along the flow of the vector field $-JY$ , that is, $(-JY)[h]=0$ for any vector field $Y\in\mathfrak{X}(P)$ and any $h\in\mathcal{C}(P)$ .

Proof.

By the definition of the compatible structure $J$ and the vector bundle map $B^{\sharp}$ , we have that

	$\displaystyle g(z)(J(z)\nabla h(z),Y(z))$	$\displaystyle=g(z)\left(B^{\sharp}(z)\mathbf{d}h(z),Y(z)\right)=B(z)(g^{\flat}(z)(Y)(z),\mathbf{d}h(z))$
		$\displaystyle=-B(z)(\mathbf{d}h(z),g^{\flat}(z)(Y)(z))=-\mathbf{d}h(z)\cdot B^{\sharp}(z)(g^{\flat}(z)(Y))(z)$
		$\displaystyle=-g(z)(\nabla h(z),B^{\sharp}(z)(g^{\flat}(z)(Y))(z))=g(z)(\nabla h(z),-J(z)Y(z)),$

where the second and fifth equalities are due to the equation (2.12) and the symmetry of the Riemannian metric $g$ , and the third equality holds thanks to the anti-symmetry of the Poisson tensor $B$ . Furthermore, let $h\in\mathcal{C}(P)$ . Then equation (2.14) implies that for all $z\in P$ ,

\displaystyle(-JY)[h](z)=\mathbf{d}h(z)\cdot(-JY)(z)=g(z)(\nabla h(z),-J(z)Y(z))=g(z)(X_{h}(z),Y(z))=0.

∎

Local coordinate representations. We start by recalling the local Lie-Darboux-Weinstein coordinates of Poisson manifolds (see the original paper [Wein 83] or [Libe 12, Vais 12] for details).

Theorem 2.19 (Lie-Darboux-Weinstein coordinates).

Let $(P,\{\cdot,\cdot\})$ be a $m$ –dimensional Poisson manifold and $z_{0}\in M$ a point where the rank of the Poisson structure equals $2n$ , $0\leq 2n\leq m$ . There exists a chart $(U,\varphi)$ of $P$ whose domain contains the point $z_{0}$ and such that the associated local coordinates, denoted by $(q^{1},\ldots,q^{n},p_{1},\ldots,p_{n},y_{1},\ldots,y_{m-2n})$ , satisfy

\{q^{i},q^{j}\}=\{p_{i},p_{j}\}=\{q^{i},y_{k}\}=\{p_{i},y_{k}\}=0,\quad\text{ and }\quad\{q^{i},p_{j}\}=\delta^{i}_{j},

for all $i,j,k$ , $1\leq i,j\leq n$ , $1\leq k\leq m-2n$ . For all $k,l$ , $1\leq k,l\leq m-2n$ , the Poisson bracket $\{y_{k},y_{l}\}$ is a function of the local coordinates $y^{1},\ldots,y^{m-2n}$ exclusively, and vanishes at $z_{0}$ . Hence, the restriction of the bracket $\{\cdot,\cdot\}$ to the coordinates $y^{1},\ldots,y^{m-2n}$ induces a Poisson structure that is usually referred to as the transverse Poisson structure of $(P,\{\cdot,\cdot\})$ at $m$ . This structure is unique up to isomorphisms.

If the rank of $(P,\{\cdot,\cdot\})$ is constant and equal to $2n$ on a neighborhood $W$ of $z_{0}$ then, by choosing the domain $U$ of the chart contained in $W$ , the coordinates $y$ satisfy

\{y_{k},y_{l}\}=0,

for all $k,l$ , $1\leq k,l\leq m-2n$ .

This theorem proves that if the rank of the Poisson structure around a given point $z_{0}\in P$ is locally constant, then there exist local coordinates $\{z^{1},\ldots,z^{m}\}$ in which the Poisson tensor $B$ can be represented by a constant matrix $\left(B^{ij}\right)_{i,j=1}^{m}$ as $B=B^{ij}\frac{\partial}{\partial z^{i}}\otimes\frac{\partial}{\partial z^{j}}$ . Hence, the Poisson bracket can be locally expressed as (using the Einstein convention in the summation of the indices)

\displaystyle\{F,G\}=B^{ij}\frac{\partial F}{\partial z^{i}}\frac{\partial G}{\partial z^{j}}.

Furthermore, for any one-form $\alpha=\alpha_{j}dz^{j}$ , the vector bundle map $B^{\sharp}$ can be locally expressed as

\displaystyle B^{\sharp}(\alpha)=B^{ij}\alpha_{i}\frac{\partial}{\partial z^{j}}.

(2.15)

Note that $g^{\flat}:TP\to T^{*}P$ is the bundle isomorphism with respect to the Riemannian metric $g:=g_{ij}dz^{i}\otimes dz^{j}$ . Thus, for any $X:=X^{j}\frac{\partial}{\partial z^{j}}$ , the bundle isomorphism $g^{\flat}$ can be expressed as

\displaystyle g^{\flat}(X)=g_{ij}X^{i}dz^{j}.

(2.16)

Therefore, combining equations (2.15) and (2.16), the compatible structure $J$ defined in (2.11) has the local representation

\displaystyle J\left(X^{j}\frac{\partial}{\partial z^{j}}\right)=B^{\sharp}\left(g^{\flat}\left(X^{j}\frac{\partial}{\partial z^{j}}\right)\right)=B^{\sharp}\left(g_{ij}X^{i}dz^{j})\right)=B^{ij}g_{ik}X^{k}\frac{\partial}{\partial z^{j}},

which proves that the matrix representation of the bundle map is given by $J=-Bg$ , that is,

\displaystyle J^{ij}=-B^{ik}g_{kj}.

(2.17)

3 Structure-preserving kernel regressions on Poisson manifolds

This section proposes a structure-preserving kernel regression method to recover Hamiltonian functions on Poisson manifolds, which is a significant generalization of the same learning problem on even-dimensional symplectic vector spaces [Hu 24]. First, in Section 3.1, we establish a differential reproducing property, which not only ensures the existence and boundedness of function gradients in the RKHS but also provides an operator framework to represent the minimizers of our structure-preserving kernel regression problems. In Section 3.3, we present an operator framework to represent the estimator for the structure-preserving regression problem. In Section 3.4, we develop a differential kernel representation of the estimator, providing a closed-form solution that is computationally feasible in applications. The error analysis of this structure-preserving estimator is conducted later in Section 4, where we establish convergence rates with both fixed and adaptive regularization parameters.

3.1 RKHS on Riemannian manifolds

The structure-preserving learning of Hamiltonian systems on symplectic vector spaces has been tackled by proving a differential reproducing property in [Hu 24]. We shall show in this section that this differential reproducing property still holds on a generic Riemannian manifold, and we will give a global version of it, where the term global means that the result that we provide is not formulated in terms of locally defined partial derivatives but of globally defined differentials as in (2.11) that do not depend on specific choices of local coordinates.

We first quickly recall Mercer kernels and the associated reproducing kernel Hilbert spaces (RKHS). A Mercer kernel on a non-empty set $\mathcal{X}$ is a positive semidefinite symmetric function $K:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ , where positive semidefinite means that it satisfies

\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}c_{j}K(x_{i},x_{j})\geq 0,

(3.1)

for any $x_{1},\ldots,x_{n}\in\mathcal{X}$ , $c_{1},\ldots,c_{n}\in\mathbb{R}$ , and any $n\in\mathbb{N}$ . Property (3.1) is equivalent to requiring that the Gram matrices $[K(x_{i},x_{j})]_{i,j=1}^{n}$ are positive semidefinite for any $x_{1},\ldots,x_{n}\in\mathcal{X}$ and any given $n\in\mathbb{N}$ . We emphasize that the positive semidefiniteness of kernels in non-Euclidean spaces requires care in the sense that the natural non-Euclidean generalization of positive definite kernels in Euclidean spaces (like the Gaussian kernel) may not be positive definite [Fera 15, Da C 23a, Da C 23b, Li 23]. However, as the next proposition shows, it is easy to construct Mercer kernels in any non-Euclidean space.

Proposition 3.1.

[Van 12, page 69] For any real-valued function $f$ on an non-empty set $\mathcal{X}$ , the function $K(x,y):=f(x)f(y)$ is a Mercer kernel on $\mathcal{X}$ . Moreover, if $K_{1},K_{2}$ are Mercer kernels on $\mathcal{X}$ , so is $K_{1}K_{2}$ .

A Mercer kernel is the key element to define a reproducing kernel Hilbert space (RKHS) as follows.

Definition 3.2 (RKHS).

Let $K:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ be a Mercer kernel on a nonempty set $\mathcal{X}$ . A Hilbert space $\mathcal{H}_{K}$ of real-valued functions on ${\cal X}$ endowed with the pointwise sum and pointwise scalar multiplication, and with inner product $\langle\cdot,\cdot\rangle_{\mathcal{H}_{K}}$ is called a reproducing kernel Hilbert space (RKHS) associated to $K$ if the following properties hold:

(i)

For all $x\in\mathcal{X}$ , we have that the function $K(x,\cdot)\in\mathcal{H}_{K}$ .

(ii)

For all $x\in\mathcal{X}$ and for all $f\in\mathcal{H}_{K}$ , the following reproducing property holds

f(x)=\langle f,K(x,\cdot)\rangle_{\mathcal{H}_{K}}.

The Moore-Aronszajn Theorem [Aron 50] establishes that given a Mercer kernel $K$ on a set ${\cal X}$ , there is a unique Hilbert space of real-valued functions $\mathcal{H}_{K}$ on ${\cal X}$ for which $K$ is a reproducing kernel, which allows us to talk about the RKHS $\mathcal{H}_{K}$ associated to $K$ . Conversely, an RKHS naturally induces a kernel function that is a Mercer kernel. Therefore, there is a bijection between RKHSs and Mercer kernels.

We now state a differential reproducing property on Riemannian manifolds whose proof can be found in Appendix A.1. The following theorem is a generalization of similar results proved for compact [Nova 18, Zhou 08], bounded [Ferr 12], and unbounded [Hu 24] subsets of Euclidean spaces.

Theorem 3.3 (Differential reproducing property on Riemannian manifolds).

Let $M$ be a complete Riemannian manifold. Let $K:M\times M\rightarrow\mathbb{R}$ be a Mercer kernel such that $K\in C_{b}^{2s+1}(M\times M)$ , for some $s\in\mathbb{N}$ (see Remark 2.14 for the definition of this space). Then the following statements hold:

(i)

For all $1\leq k\leq s$ , we have that the function $D^{(k,0)}K(y,\cdot)\cdot v\in\mathcal{H}_{K}$ for all $y\in T^{k-1}M$ and $v\in T_{y}(T^{k-1}M)$ .

(ii)

A differential reproducing property holds true for all $f\in\mathcal{H}_{K}$ :

D^{k}f(y)\cdot v=\langle D^{(k,0)}K(y,\cdot)\cdot v,f\rangle_{\mathcal{H}_{K}},\quad\mbox{for all $y\in T^{k-1}M$, $v\in T_{y}(T^{k-1}M)$, and $1\leq k\leq s$.}

(3.2)

(iii)

Denote $\kappa^{2}=(s+1)\|K\|_{C_{b}^{s}}$ . The inclusion $J:\mathcal{H}_{K}\hookrightarrow C_{b}^{s}(M)$ is well-defined and bounded:

\displaystyle\|f\|_{C_{b}^{s}}\leqslant\kappa\|f\|_{\mathcal{H}_{K}},\quad\forall~{}f\in\mathcal{H}_{K}.

(3.3)

For convenience, we use the following notation in what follows: given $x\in M$ , we denote by $K_{x}:M\rightarrow\mathbb{R}$ the function on $M$ given by $K_{x}(y):=K(x,y)$ and we call it the kernel section of $K$ at $x\in M$ . Therefore, $K_{\cdot}$ can be regarded as a function on $M\times M$ such that $K_{x}(y):=K(y,x)$ . For any function $f\in\mathcal{H}_{K}$ , we denote by $\langle f,K_{\cdot}\rangle_{\mathcal{H}_{K}}$ the function in $\mathcal{H}_{K}$ given by $\langle f,K_{y}\rangle_{\mathcal{H}_{K}}=f(y)$ (by the reproducing property). Furthermore, for a Poisson manifold $P$ equipped with a Riemannian metric $g$ , we shall also denote by $X_{K_{\cdot}}$ the family of Hamiltonian vector fields parametrized by the elements in $P$ , that is, given any $z\in P$ , $X_{K_{z}}$ is the Hamiltonian vector field of the function $K_{z}$ . In particular, for any $z\in P$ and any vector $v\in T_{z}P$ , $g(z)\left(X_{K_{\cdot}}(z),v\right)$ is a function on the manifold $P$ .

Remark 3.4.

(i) By Theorem 3.3 (i), $D^{(1,0)}K(x,\cdot)\cdot u\in\mathcal{H}_{K}$ for any $u\in T_{x}M$ . Thus, the differential reproducing property (3.2) with $f=D^{(1,0)}K(x,\cdot)\cdot u$ implies that

\displaystyle D^{(1,1)}K(x,y)\cdot(u,v)=\left\langle D^{(1,0)}K(x,\cdot)\cdot u,D^{(1,0)}K(y,\cdot)\cdot v\right\rangle_{\mathcal{H}_{K}},\quad\mbox{for all $y\in M,v\in T_{y}M$.}

Furthermore, in the presence of a Poisson structure and combining with Proposition 2.18, we obtain that for any vector field $Y\in\mathfrak{X}(M)$ ,

\displaystyle g(z)\left(X_{K_{\cdot}}(z),Y(z)\right)=D^{(1,0)}K(z,\cdot)\cdot(-J(z)Y(z))\in\mathcal{H}_{K}.

(ii) Note that unlike the differential reproducing property in [Hu 24] that was formulated in terms of partial derivatives, (3.2) is a coordinate-free expression that has been written down in terms of global differentials.

Corollary 3.5.

Let $M$ be a Riemannian manifold. Denote by $\mathcal{H}_{K}$ the RKHS associated with the kernel $K:M\times M\to\mathbb{R}$ . Then for each vector field $X\in\mathfrak{X}(M)$ and each $f\in\mathcal{H}_{K}$ , we have that

\displaystyle X[f](x)=X[\langle f,K_{\cdot}\rangle_{\mathcal{H}_{K}}](x),\quad\forall\ x\in M.

Proof.

By the reproducing property, we have that for each $f\in\mathcal{H}_{K}$ , $f(x)=\langle f,K_{x}\rangle_{\mathcal{H}_{K}}$ , for all $x\in M$ . Let $F_{t}$ be the flow of the vector field $X$ . Then,

X[\langle f,K_{\cdot}\rangle_{\mathcal{H}_{K}}](x):=\left.\frac{d}{dt}\right|_{t=0}[\langle f,K_{F_{t}(x)}\rangle_{\mathcal{H}_{K}}]=\left.\frac{d}{dt}\right|_{t=0}f(F_{t}(x))=X[f](x).

∎

3.2 The structure-preserving learning problem

As we already explained in the introduction, the data that we are given for the learning problem are noisy realizations of the Hamiltonian vector field $X_{H}\in\mathfrak{X}(P)$ whose Hamiltonian function $H\in C^{\infty}(P)$ we are trying to estimate, that is,

\displaystyle\mathbf{X}_{\sigma^{2}}^{(n)}:=X_{H}(\mathbf{Z}^{(n)})+\bm{\varepsilon}^{(n)},\quad n=1,\ldots,N,

where $\mathbf{Z}^{(n)}$ are IID random variables with values in the Poisson manifold $P$ with the same distribution $\mu_{\mathbf{Z}}$ and $\bm{\varepsilon}^{(n)}$ are independent random variables on the tangent spaces $T_{\mathbf{z}^{(n)}}P$ with mean zero and variance $\sigma^{2}I_{d}$ (see Remark 2.16). We notice that unlike the Euclidean case studied in [Hu 24] where it was supposed that $\bm{\varepsilon}^{(n)}$ have the same distribution in $\mathbb{R}^{d}$ , in the present situation, the tangent spaces $T_{\mathbf{z}^{(n)}}P$ are different vector spaces even though they are isomorphic and hence, we can not assume that the random variables $\bm{\varepsilon}^{(n)}$ in different spaces have the same distribution. However, the independence hypothesis and the assumptions on the moments are sufficient for our structure-preserving kernel regression.

Notation. We shall write the random samples as:

\displaystyle\begin{split}\mathbf{Z}_{N}&:=\mathrm{Vec}\left(\mathbf{Z}^{(1)}|\cdots|\mathbf{Z}^{(N)}\right)\in\Pi_{i=1}^{N}P,\\ \mathbf{X}_{\sigma^{2},N}&:=\mathrm{Vec}\left(\mathbf{X}_{\sigma^{2}}^{(1)}|\cdots|\mathbf{X}_{\sigma^{2}}^{(N)}\right)\in\Pi_{i=1}^{N}T_{\mathbf{Z}^{(i)}}P,\end{split}

where the symbol ‘ $\mathrm{Vec}$ ’ stands for the vectorization of the corresponding matrices. We shall denote by $\mathbf{z}^{(n)}$ and $\mathbf{x}_{\sigma^{2}}^{(n)}$ the realizations of the random variables $\mathbf{Z}^{(n)}$ and $\mathbf{X}_{\sigma^{2}}^{(n)}$ , respectively. The collection of realizations is denoted by

\displaystyle\begin{split}\mathbf{z}_{N}&:=\mathrm{Vec}\left(\mathbf{z}^{(1)}|\cdots|\mathbf{z}^{(N)}\right)\in\Pi_{i=1}^{N}P,\\ \mathbf{x}_{\sigma^{2},N}&:=\mathrm{Vec}\left(\mathbf{x}_{\sigma^{2}}^{(1)}|\cdots|\mathbf{x}_{\sigma^{2}}^{(N)}\right)\in\Pi_{i=1}^{N}T_{\mathbf{z}^{(i)}}P.\end{split}

In the sequel, if $f:P\to\mathbb{R}^{s}$ is a function, we then shall denote the value $\mathrm{Vec}\left(f(\mathbf{Z}^{(1)})|\cdots|f(\mathbf{Z}^{(N)})\right)\in\mathbb{R}^{sN}$ by $f(\mathbf{Z}_{N})$ .

Structure-preserving kernel ridge regression. The strategy behind structure-preservation in the context of kernel regression is that we search the vector field $X\in\mathfrak{X}(P)$ that minimizes a risk functional among those that have the form $X:=X_{h}$ , where $h:P\rightarrow\mathbb{R}$ belongs to the RKHS $\mathcal{H}_{K}$ associated with a kernel $K$ defined on the Poisson manifold $P$ . This approach obviously guarantees that the learned vector field is Hamiltonian. To make the method explicit, we shall be solving the following optimization problem

	$\displaystyle\widehat{h}_{\lambda,N}$	$\displaystyle:=\mathop{\arg\min}\limits_{h\in\mathcal{H}_{K}}\ \widehat{R}_{\lambda,N}(h),$		(3.4)
	$\displaystyle\widehat{R}_{\lambda,N}(h)$	$\displaystyle:=\frac{1}{N}\sum_{n=1}^{N}\\|X_{h}(\mathbf{Z}^{(n)})-\mathbf{X}^{(n)}_{\sigma^{2}}\\|_{g}^{2}+\lambda\\|h\\|_{\mathcal{H}_{K}}^{2},$		(3.5)

where $X_{h}$ is the Hamiltonian vector field of $h$ and $\lambda\geq 0$ is the Tikhonov regularization parameter. We shall refer to the minimizer $\widehat{h}_{\lambda,N}$ as the structure-preserving kernel estimator of the Hamiltonian function $H$ . The functional $\widehat{R}_{\lambda,N}$ is called the regularized empirical risk.

Remark 3.6.

(i) If $K\in C_{b}^{3}(P\times P)$ then the differential reproducing property in Theorem 3.3 implies that the functions in ${\mathcal{H}}_{K}$ are all differentiable because and hence they always have a Hamiltonian vector field associated. This implies that the regularized empirical risk $\widehat{R}_{\lambda,N}$ defined in (3.5) and the associated optimization problem (3.4) is always well-defined in that situation.

(ii) The minimization problem (3.4) might be ill-posed since for any Casimir function $C\in\mathcal{C}\cap\mathcal{H}_{K}$ and $h\in\mathcal{H}_{K}$ , the function $h+C$ shares the same Hamiltonian vector field as $h$ . This could result in the lack of unique solutions for the minimization problem. Later on in Remark 3.14, we show how ridge regularization solves this problem.

Given the probability measure $\mu_{\mathbf{Z}}$ and the Riemannian metric $g$ on the manifold $P$ , a vector field $Y\in\mathfrak{X}(P)$ is said $L^{2}(\mu_{\mathbf{Z}})$ -integrable if the following $L^{2}(\mu_{\mathbf{Z}})$ norm of $Y$ is finite, that is,

\displaystyle\|Y\|_{L^{2}(\mu_{\mathbf{Z}})}:=\int_{P}g(z)(Y(z),Y(z))\mathrm{d}\mu(z)<\infty.

(3.6)

Denote by $L^{2}(\mathfrak{X}(P),g,\mu_{\mathbf{Z}})$ the space consisting of all $L^{2}(\mu_{\mathbf{Z}})$ -integrable vector fields. The measure-theoretic analogue $R_{\lambda}$ of (3.5) is referred to as regularized statistical risk and is defined by

\displaystyle R_{\lambda}(h)

\displaystyle:=\|X_{h}-X_{H}\|_{L^{2}(\mu_{\mathbf{Z}})}^{2}+\lambda\|h\|_{\mathcal{H}_{K}}^{2}+\sigma^{2},

(3.7)

where the $L^{2}(\mu_{\mathbf{Z}})$ norm is defined in (3.6) and $X_{h},X_{H}$ are the Hamiltonian vector fields of $h$ and $H$ , respectively. We denote by $h^{*}_{\lambda}\in\mathcal{H}_{K}$ the best-in-class functional with the minimal regularized statistical risk, that is,

\displaystyle h^{*}_{\lambda}:=\mathop{\arg\min}\limits_{h\in\mathcal{H}_{K}}R_{\lambda}(h).

(3.8)

We note that by the strong law of large numbers, the regularized empirical and statistical risks are consistent within the RKHS, that is, for every $h\in\mathcal{H}_{K}$ , we have that

\displaystyle\lim\limits_{N\to\infty}\mathbb{E}_{\bm{\varepsilon}}\left[\widehat{R}_{\lambda,N}(h)\right]=R_{\lambda}(h),\quad\mbox{almost surely},

where $\mathbb{E}_{\bm{\varepsilon}}$ means taking the conditional expectation for all random variables $\bm{\varepsilon}^{(n)}$ given all $\mathbf{Z}^{(n)}$ .

3.3 Operator representations of the kernel estimator

In order to find and study the solutions of the inverse learning problems (3.4)-(3.5) and (3.7)-(3.8), we introduce the operators $A:\mathcal{H}_{K}\to\mathfrak{X}_{\mathrm{Ham}}(P)$ and $A_{N}:\mathcal{H}_{K}\to T_{\mathbf{Z}_{N}}P:=\Pi_{i=1}^{N}T_{\mathbf{Z}^{(i)}}P$ as

\displaystyle Ah:=X_{h}=J\nabla h,\text{ and }A_{N}h:=X_{h}(\mathbf{Z}_{N})=\frac{1}{\sqrt{N}}\mathbb{J}_{N}\nabla h(\mathbf{Z}_{N}),\quad\text{ for all }h\in\mathcal{H}_{K},

(3.9)

where $\mathbb{J}_{N}=\operatorname{diag}\{J(\mathbf{Z}^{(1)}),\ldots,J(\mathbf{Z}^{(N)})\}$ and $J$ is the compatible structure defined in (2.11). By Theorem 3.3, if the Mercer kernel $K\in C_{b}^{3}(P\times P)$ , then the functions in $\mathcal{H}_{K}$ are differentiable, and hence the operators $A$ and $A_{N}$ are well-defined. We aim to show that the image of the operator $A$ belongs to the space $L^{2}(\mathfrak{X}(P),g,\mu_{\mathbf{Z}})$ . To achieve this, we must make the following boundedness assumption on the compatible structure $J$ .

Assumption 3.7.

The compatible structure $J$ defined in (2.11) satisfies

\displaystyle g(z)(J(z)v,J(z)v)\leq\gamma(z)g(z)(v,v),\quad\text{for all }z\in P,\text{ and }v\in T_{z}P,

(3.10)

where $\gamma$ is a positive function on the Poisson manifold $P$ bounded above by some constant $C^{2}$ ( $C>0$ ).

Using this boundedness condition on the compatible structure $J$ , we can obtain the boundedness of the operators $A$ and $A_{N}$ defined in (3.9) as maps from $\mathcal{H}_{K}$ to $L^{2}(\mathfrak{X}(P),g,\mu_{\mathbf{Z}})$ and from $\mathcal{H}_{K}$ to $T_{\mathbf{Z}_{N}}P$ , respectively. A detailed proof of the following proposition can be found in the Appendix A.2

Proposition 3.8.

Let $P$ be a Poisson manifold and let $K\in C_{b}^{3}(P\times P)$ be a Mercer kernel. Suppose that Assumption 3.7 holds. Then, the operator $A$ defined in (3.9) is a bounded linear operator that maps $\mathcal{H}_{K}$ into $L^{2}(\mathfrak{X}(P),g,\mu_{\mathbf{Z}})$ with an operator norm $\|A\|$ that satisfies $\|A\|\leq\kappa C$ , where $\kappa^{2}=2\|K\|_{C_{b}^{2}}$ . The adjoint operator $A^{*}:L^{2}(\mathfrak{X}(P),g,\mu_{\mathbf{Z}})\longrightarrow{\mathcal{H}}_{K}$ of $A$ is given by

\displaystyle A^{*}Y=\int_{P}g(z)(X_{K_{\cdot}}(z),Y(z))\mathrm{d}\mu_{\mathbf{Z}}(z),\quad\mbox{for all $Y\in L^{2}(\mathfrak{X}(P),g,\mu_{\mathbf{Z}})$.}

As a consequence, the bounded linear operator $Q:\mathcal{H}_{K}\longrightarrow\mathcal{H}_{K}$ , defined by

\displaystyle Qh:=A^{*}Ah=\int_{P}g(z)(X_{K_{\cdot}}(z),X_{h}(z))\mathrm{~{}d}\mu_{\mathbf{Z}}(z),

(3.11)

is a positive trace class operator that satisfies $\operatorname{Tr}(Q)\leq\frac{1}{2}d\kappa^{2}C^{2}$ .

Remark 3.9.

(i) The proof of Proposition 3.8 shows that the boundedness condition of the compatible structure $J$ in Assumption 3.7 can be weakened by requiring that the function $\gamma$ is positive $L^{1}$ -integrable with respect to the probability measure $\mu_{\mathbf{Z}}$ . This condition obviously holds if $\mu_{\mathbf{Z}}$ is compactly supported or it is a Gaussian probability measure, which shows this boundedness condition holds for a large class of compatible structures.

(ii) When the Poisson manifold $P$ is symplectic and the compatible structure $J$ satisfies $J^{2}=-I$ then it is a complex structure and $P$ becomes a Kähler manifold. Even-dimensional vector spaces are special cases of Kähler manifolds by picking as the complex structure the one naturally associated to the canonical (Darboux) symplectic form. In these cases, the boundedness condition (3.10) naturally holds by choosing $\gamma\equiv 1$ . This is why this condition never appears in [Hu 24].

We now deal with the empirical version $A_{N}$ of the operator $A$ defined in (3.9). For convenience, we define the inner product $g_{N}(\cdot,\cdot)$ on the product space $T_{\mathbf{Z}_{N}}P$ as follows:

\displaystyle g_{N}(u,v):=\sum_{i=1}^{N}g\left(\mathbf{Z}^{(i)}\right)(u_{i},v_{i}),\quad u=(u_{1},\ldots,u_{N}),v=(v_{1},\ldots,v_{N})\in T_{\mathbf{Z}_{N}}P,

(3.12)

where $u_{i},v_{i}\in T_{\mathbf{Z}^{(i)}}P$ for $i=1,\ldots,N$ . Furthermore, we denote by $\|\cdot\|_{g_{N}}$ the corresponding norm.

Proposition 3.10.

Suppose that Assumption 3.7 holds. Then, the operator $A_{N}:\mathcal{H}_{K}\rightarrow T_{\mathbf{Z}_{N}}P$ defined in (3.9) is a bounded linear operator whose operator norm $\|A_{N}\|\leq\kappa C$ . The adjoint operator $A_{N}^{*}:T_{\mathbf{Z}_{N}}P\rightarrow\mathcal{H}_{K}$ of $A_{N}$ is a finite rank operator given by

\displaystyle A^{*}_{N}W=\frac{1}{\sqrt{N}}g_{N}(W,X_{K_{\cdot}}(\mathbf{Z}_{N})),\quad W\in T_{\mathbf{Z}_{N}}P.

Moreover, the operator $Q_{N}$ defined by

\displaystyle Q_{N}h:=A^{*}_{N}A_{N}h=\frac{1}{N}g_{N}(X_{h}(\mathbf{Z}_{N}),X_{K_{\cdot}}(\mathbf{Z}_{N})),\quad h\in\mathcal{H}_{K},

(3.13)

is a positive-semidefinite compact operator.

Proof.

The explicit expressions of $A_{N}^{\ast}$ and $Q_{N}$ follow from a direct calculation. Under Assumption 3.7 it follows from the equation (A.4) that

\|A_{N}h\|^{2}=\frac{1}{N}\|\mathbb{J}_{N}\nabla h(\mathbf{Z}_{N})\|_{g_{N}}^{2}\leq C^{2}\kappa^{2}\|h\|^{2}_{\mathcal{H}_{K}},

which implies that $A_{N}$ is bounded and that $\|A_{N}\|\leq\kappa C$ . The compactness of the operator $Q_{N}$ follows from a proof similar to that of Proposition 3.8 for $Q$ . ∎

Having defined the operators $A$ and $A_{N}$ , we obtain an operator representation of the minimizers of the learning problems (3.4)-(3.5) and (3.7)-(3.8) that is similar to the one that was already stated in in [Hu 24].

Proposition 3.11.

Let $\widehat{h}_{\lambda,N}$ and $h_{\lambda}^{*}$ be the minimizers of (3.5) and (3.7), respectively. Then, for all $\lambda>0$ , the minimizers $h_{\lambda}^{*}$ and $\widehat{h}_{\lambda,N}$ are unique and they are given by

	$\displaystyle h_{\lambda}^{*}:$	$\displaystyle=(Q+\lambda I)^{-1}A^{*}X_{H}.$
	$\displaystyle\widehat{h}_{\lambda,N}:$	$\displaystyle=\frac{1}{\sqrt{N}}(Q_{N}+\lambda I)^{-1}A_{N}^{}\mathbf{X}_{\sigma^{2},N},\mbox{ with }Q_{N}=A_{N}^{}A_{N}.$		(3.14)

3.4 The global solution of the learning problem

We now derive a kernel representation for the solution of the learning problem (3.4)-(3.5). More precisely, we derive a closed-form expression for the estimator of a Hamiltonian function of the data-generating process by defining a generalized differential Gram matrix on manifolds. We shall see that the symmetric and positive semidefinite properties of this generalized differential Gram matrix guarantee the uniqueness of the estimator. The main expression that will be obtained is stated below in Theorem 3.13 and amounts to a differential version of the Representer Theorem for Poisson manifolds and that is how we shall refer to it. In order to derive it, we first define a generalized differential Gram matrix $G_{N}:T_{\mathbf{Z}_{N}}P\rightarrow T_{\mathbf{Z}_{N}}P$ as

\displaystyle G_{N}\mathbf{c}:=X_{g_{N}(\mathbf{c},X_{K_{\cdot}}(\mathbf{Z}_{N}))}(\mathbf{Z}_{N}),\quad\mathbf{c}\in T_{\mathbf{Z}_{N}}P.

(3.15)

As it was studied in [Hu 24], in symplectic vector spaces endowed with the canonical symplectic form $\omega_{\mathrm{can}}$ , the generalized differential Gram matrix $G_{N}$ defined in (3.15) reduces to the usual differential Gram matrix $\mathbb{J}_{\mathrm{can}}\nabla_{1,2}K(\mathbf{Z}_{N},\mathbf{Z}_{N})\mathbb{J}_{\mathrm{can}}^{\top}$ , where $\nabla_{1,2}K$ represents the matrix of partial derivatives of $K$ with respect to the first and second arguments.

The general differential Gram matrix $G_{N}$ defined in (3.15) corresponds, roughly speaking, to the operator $Q_{N}$ defined in (3.13) (see also [Hu 24]). In the following proposition that is proved in the appendix, we see that the generalized differential Gram matrix $G_{N}$ is symmetric and positive semidefinite.

Proposition 3.12.

Given a Mercer kernel $K\in C_{b}^{3}(P\times P)$ on the $d$ -dimensional Poisson manifold $(P,\left\{\cdot,\cdot\right\})$ and endowed with a Riemannian metric $g$ , the generalized differential Gram matrix $G_{N}$ defined in (3.15) is symmetric and positive semidefinite with respect to the metric $g_{N}$ in (3.12), that is,

g_{N}(\mathbf{c},G_{N}\mathbf{c}^{\prime})=g_{N}(\mathbf{c}^{\prime},G_{N}\mathbf{c})~{}\text{ and }~{}g_{N}(\mathbf{c},G_{N}\mathbf{c})\geq 0~{}\text{ for all }\mathbf{c},\mathbf{c}^{\prime}\in T_{\mathbf{Z}_{N}}P.

Theorem 3.13 (Differential Representer Theorem on Poisson manifolds).

Let $K\in C_{b}^{3}(P\times P)$ be a Mercer kernel on the Poisson manifold $(P,\left\{\cdot,\cdot\right\})$ and endowed with a Riemannian metric $g$ . Suppose that Assumption 3.7 holds. Then for every $\lambda>0$ , the estimator $\widehat{h}_{\lambda,N}$ introduced in (3.4) as the minimizer of the regularized empirical risk functional $\widehat{R}_{\lambda,N}$ in (3.5) can be written as

\widehat{h}_{\lambda,N}=g_{N}(\widehat{\bf c},X_{K_{\cdot}}(\mathbf{Z}_{N})),

(3.16)

where $\widehat{\mathbf{c}}\in T_{\mathbf{Z}_{N}}P$ is given by

\displaystyle\widehat{\mathbf{c}}=(G_{N}+\lambda NI)^{-1}\mathbf{X}_{\sigma^{2},N}.

Proof.

Denote the space $\mathcal{H}_{K}^{N}$ by

\mathcal{H}_{K}^{N}:=\left\{g_{N}(X_{K_{\cdot}}(\mathbf{Z}_{N}),\mathbf{c}),\ \mathbf{c}\in T_{\mathbf{Z}_{N}}P\right\},

By part (i) of Theorem 3.3, $\mathcal{H}_{K}^{N}$ is a subspace of $\mathcal{H}_{K}$ . Then by the representation of the operator $Q_{N}$ in the equation (3.13), we know that $Q_{N}(\mathcal{H}_{K}^{N})\subseteq\mathcal{H}_{K}^{N}$ , that is, $\mathcal{H}_{K}^{N}$ is an invariant space for the operator $Q_{N}$ . This implies that, for any $\lambda>0$ , $(Q_{N}+\lambda I)(\mathcal{H}_{K}^{N})\subseteq\mathcal{H}_{K}^{N}$ . By Proposition 3.12, the operator $Q_{N}$ is positive semi-definite, we can conclude that the restriction $(Q_{N}+\lambda I)|_{\mathcal{H}_{K}^{N}}$ is invertible and since the space $\mathcal{H}_{K}^{N}$ is finite-dimensional then it is also an invariant subspace of $(Q_{N}+\lambda I)|_{\mathcal{H}_{K}^{N}}^{-1}$ , that is

(Q_{N}+\lambda I)|_{\mathcal{H}_{K}^{N}}^{-1}\left(\mathcal{H}_{K}^{N}\right)\subset\mathcal{H}_{K}^{N}.

Then, there exists a vector $\widehat{\mathbf{c}}\in T_{\mathbf{Z}_{N}}P$ such that

\widehat{h}_{\lambda,N}=g_{N}(\widehat{\mathbf{c}},X_{K_{\cdot}}(\mathbf{Z}_{N})).

(3.17)

Then, applying $(Q_{N}+\lambda I)$ on the left of both sides on the estimator $\widehat{h}_{\lambda,N}$ in (3.14) and plugging (3.17) into the identity, we obtain that

\displaystyle\frac{1}{N}g_{N}\left(X_{K_{\cdot}}(\mathbf{Z}_{N}),\mathbf{X}_{\sigma^{2},N}\right)=\frac{1}{N}g_{N}\left(X_{\widehat{h}_{\lambda,N}}(\mathbf{Z}_{N}),X_{K_{\cdot}}(\mathbf{Z}_{N})\right)+\lambda g_{N}\left(\widehat{\mathbf{c}},X_{K_{\cdot}}(\mathbf{Z}_{N})\right).

(3.18)

Therefore, we obtain that

\displaystyle G_{N}\widehat{\mathbf{c}}+\lambda N\widehat{\mathbf{c}}=\mathbf{X}_{\sigma^{2},N}.

Since the matrix $G_{N}+\lambda NI$ is invertible due to the positive semi-definiteness of $G_{N}$ that we proved in Proposition 3.12, we can write that

\displaystyle\widehat{\mathbf{c}}=(G_{N}+\lambda NI)^{-1}\mathbf{X}_{\sigma^{2},N}.

(3.19)

This shows that the function $\widehat{h}_{\lambda,N}$ in (3.17) with $\widehat{\mathbf{c}}$ determined by (3.19) coincides with minimizer (3.14) of the regularized empirical risk functional $\widehat{R}_{\lambda,N}$ in (3.5). Since by Proposition 3.11, this minimizer is unique, the result follows. ∎

Remark 3.14 (On the uniqueness of the estimated Hamiltonian).

Define the kernel of the operator $A$ as follows:

\displaystyle\mathcal{H}_{\mathrm{null}}:=\{h\in\mathcal{H}_{K}\mid Ah=J\nabla h=0\}.

It can be checked that $\mathcal{H}_{\mathrm{null}}$ is a closed vector subspace of $\mathcal{H}_{K}$ . Indeed, let $\left\{f_{n}\right\}_{n\in\mathbb{N}}$ be a convergent sequence in $\mathcal{H}_{\mathrm{null}}$ such that $\left\|f_{n}-f\right\|_{{\mathcal{H}}_{K}}\rightarrow 0$ as $n\rightarrow\infty$ , for some $f\in{\mathcal{H}}_{K}$ . Then, by (A.4) in the proof of Proposition 3.8 and Assumption 3.7, it holds that

\sup_{z\in P}\|X_{f}(z)\|_{g}^{2}=\sup_{z\in P}\|X_{\left(f_{n}-f\right)}(z)\|_{g}^{2}\leq\kappa^{2}C^{2}\|f_{n}-f\|_{\mathcal{H}_{K}}^{2}\rightarrow 0,\quad\mbox{as $n\rightarrow\infty$,}

which implies that $Af=0$ and hence $\mathcal{H}_{\mathrm{null}}$ is a closed subspace of ${\mathcal{H}}_{K}$ . Therefore, $\mathcal{H}_{K}$ can be decomposed as

\displaystyle\mathcal{H}_{K}=\mathcal{H}_{\mathrm{null}}\oplus\mathcal{H}_{\mathrm{null}}^{\bot},

where $\mathcal{H}_{\mathrm{null}}^{\bot}$ denotes the orthogonal complement of $\mathcal{H}_{\mathrm{null}}$ with respect to the RKHS inner product $\left\langle\cdot,\cdot\right\rangle_{\mathcal{H}_{K}}$ . Using this decomposition and the expression of the kernel estimator (3.16), we conclude that

\displaystyle\widehat{h}_{\lambda,N}\in\mathcal{H}_{\mathrm{null}}^{\bot},

(3.20)

which is due to the fact that for any $h\in\mathcal{H}_{\mathrm{null}}$ ,

	$\displaystyle\langle\widehat{h}_{\lambda,N},h\rangle_{\mathcal{H}_{K}}$	$\displaystyle=\langle g_{N}(X_{K_{\cdot}}(\mathbf{Z}_{N}),\widehat{\mathbf{c}}),h\rangle_{\mathcal{H}_{K}}=\langle D^{(1,0)}K(\mathbf{Z}_{N},\cdot)\cdot(-\mathbb{J}_{N}\widehat{\mathbf{c}}),h\rangle_{\mathcal{H}_{K}}$
		$\displaystyle=Dh(\mathbf{Z}_{N})\cdot(-\mathbb{J}_{N}\widehat{\mathbf{c}})=g_{N}(X_{h}(\mathbf{Z}_{N}),\widehat{\mathbf{c}})=0,$

where the second and fourth equalities follow by Proposition 2.18, and the third equality is due to the differential reproducing property (3.2) in Theorem 3.3.

In general, the space $\mathcal{H}_{\mathrm{null}}$ could contain Casimir functions in $\mathcal{C}(P)$ and this could, in principle, cause ill-posedness of the learning problem since adding Casimir functions in $\mathcal{H}_{\mathrm{null}}$ to $\widehat{h}_{\lambda,N}$ would not change the corresponding Hamiltonian vector field. Therefore, one may wonder why, according to Theorem 3.13, the estimator $\widehat{h}_{\lambda,N}$ is unique but not up to Casimir functions. The explanation is in the regularization term. Indeed, let $\widehat{h}_{\lambda,N}$ be the minimizer in (3.16) and let $h\in\mathcal{H}_{\mathrm{null}}$ . Even though $\widehat{h}_{\lambda,N}$ and $\widehat{h}_{\lambda,N}+h$ have the same Hamiltonian vector field associated, it is easy to show that $\widehat{h}_{\lambda,N}+h$ is a minimizer of (3.4) if and only if $h\equiv 0$ . This is because of (3.20) and

	$\displaystyle\widehat{R}_{\lambda,N}(\widehat{h}_{\lambda,N}+h)$	$\displaystyle:=\frac{1}{N}\sum_{n=1}^{N}\\|X_{\widehat{h}_{\lambda,N}+h}(\mathbf{Z}^{(n)})-\mathbf{X}^{(n)}_{\sigma^{2}}\\|_{g}^{2}+\lambda\\|\widehat{h}_{\lambda,N}+h\\|_{\mathcal{H}_{K}}^{2}$
		$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\\|X_{\widehat{h}_{\lambda,N}}(\mathbf{Z}^{(n)})-\mathbf{X}^{(n)}_{\sigma^{2}}\\|_{g}^{2}+\lambda\left(\\|\widehat{h}_{\lambda,N}\\|_{\mathcal{H}_{K}}^{2}+\\|h\\|^{2}_{\mathcal{H}_{K}}+2\langle\widehat{h}_{\lambda,N},h\rangle_{\mathcal{H}_{K}}\right)$
		$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\\|X_{\widehat{h}_{\lambda,N}}(\mathbf{Z}^{(n)})-\mathbf{X}^{(n)}_{\sigma^{2}}\\|^{2}+\lambda\left(\\|\widehat{h}_{\lambda,N}\\|_{\mathcal{H}_{K}}^{2}+\\|h\\|^{2}_{\mathcal{H}_{K}}\right).$

Local coordinate representation of the kernel estimator. In applications, it is important to have local coordinate expressions for the kernel estimator $\widehat{h}_{\lambda,N}$ . Recall that by Theorem 3.13, the kernel estimator is given by

\widehat{h}_{\lambda,N}=g_{N}(\widehat{\bf c},X_{K_{\cdot}}(\mathbf{Z}_{N}))=g_{N}(\widehat{\bf c},\mathbb{J}_{N}\nabla K(\mathbf{Z}_{N},\cdot)),

where $\widehat{\mathbf{c}}\in T_{\mathbf{Z}_{N}}P$ is given by

\displaystyle\widehat{\mathbf{c}}=(G_{N}+\lambda NI)^{-1}\mathbf{X}_{\sigma^{2},N}.

Since $J$ can be locally represented as (2.17), it suffices to write down the Gram matrix $G_{N}$ locally. Let $\mathbf{z}^{(1)},\ldots,\mathbf{z}^{(N)}$ be the available data points and let $\{z_{n}^{1},\ldots,z_{n}^{d}\}$ , $n=1,\ldots,N$ , be local coordinates on open sets containing the points $\mathbf{z}^{(n)}$ , for all $n=1,\ldots,N$ . Let $g^{(n)}$ be the matrix associated to the Riemannian metric $g$ in the local coordinates $\{z_{n}^{1},\ldots,z_{n}^{d}\}$ . We first compute the following function for each $\mathbf{c}\in T_{\mathbf{z}_{N}}P$ ,

\displaystyle g_{N}(\mathbf{c},X_{K_{\cdot}}(\mathbf{z}_{N}))=g_{N}(\mathbf{c},\mathbb{J}_{N}\nabla K(\mathbf{z}_{N},\cdot))=\sum_{i=1}^{N}(\partial_{1}K)^{\top}(\mathbf{z}^{(i)},\cdot){g^{(i)}}^{-1}J^{\top}(\mathbf{z}^{(i)})g^{(i)}\mathbf{c}^{i}.

By the definition of the Gram matrix in (3.15) and combining it with the local representation of $J$ in (2.17), we obtain that

	$\displaystyle(G_{N}\mathbf{c})^{i}$	$\displaystyle=X_{g_{N}(\mathbf{c},X_{K_{\cdot}}(\mathbf{z}_{N}))}(\mathbf{z}^{(i)})=\sum_{j=1}^{N}J(\mathbf{z}^{(i)}){g^{(i)}}^{-1}\partial_{1,2}K(\mathbf{z}^{(i)},\mathbf{z}^{(i)}){g^{(j)}}^{-1}J^{\top}(\mathbf{z}^{(j)})g^{(j)}\mathbf{c}^{j}$
		$\displaystyle=\sum_{j=1}^{N}(-B(\mathbf{z}^{(i)})g^{(i)}){g^{(i)}}^{-1}\partial_{1,2}K(\mathbf{z}^{(j)},\mathbf{z}^{(i)}){g^{(j)}}^{-1}(-{g^{(j)}}^{\top}B^{\top}(\mathbf{z}^{(j)}))g^{(j)}\mathbf{c}^{j}$
		$\displaystyle=\sum_{j=1}^{N}B(\mathbf{z}^{(i)})\partial_{1,2}K(\mathbf{z}^{(j)},\mathbf{z}^{(i)})B^{\top}(\mathbf{z}^{(j)})g^{(j)}\mathbf{c}^{j},$

which shows that the $(i,j)$ -component matrix $G_{N}^{i,j}$ of the general Gram matrix is given by

\displaystyle G_{N}^{i,j}=B(\mathbf{z}^{(i)})\partial_{1,2}K(\mathbf{z}^{(i)},\mathbf{z}^{(j)})B^{\top}(\mathbf{z}^{(j)})g^{(j)}.

(3.21)

3.5 Symmetries and momentum conservation by kernel estimated Hamiltonians

Studies have been carried out to bridge invariant feature learning with kernel methods [Mrou 15]. In this subsection, we will see that if the unknown Hamiltonian function is invariant under a canonical group action, then the structure-preserving kernel estimator can also be made invariant under the same action by imposing certain symmetry constraints on the kernel $K$ . As a result, the fibers of any momentum map associated with that group action that are preserved by the flow of the ground-truth Hamiltonian $H$ will also be preserved by the flow of the estimator $\widehat{h}_{\lambda,N}$ .

Definition 3.15.

Left $G$ be a group acting on the manifold $P$ . A function $h:P\rightarrow\mathbb{R}$ on the Poisson manifold $P$ is called $G$ -invariant if

\displaystyle h(g\cdot x)=h(x),\quad\text{ for all $g,g^{\prime}\in G$ and all $x\in P$.}

A kernel $K:P\times P\rightarrow\mathbb{R}$ is called argumentwise invariant by the group action if

\displaystyle K(g\cdot x,g^{\prime}\cdot x^{\prime})=K(x,x^{\prime}),\quad\text{ for all $g,g^{\prime}\in G$ and all $x,x^{\prime}\in P$.}

Argumentwise invariant kernels can be constructed by double averaging (using the Haar measure if the group is compact) or, alternatively, if $K$ is an arbitrary Mercer kernel and $f$ is a $G$ -invariant function. Then the function $K^{\prime}$ defined as $K^{\prime}(x,y)=K(f(x),f(y))$ is an argumentwise invariant kernel.

Proposition 3.16 (Noether’s Theorem for the structure-preserving kernel estimator).

Let $G$ be a Lie group acting canonically on the Poisson manifold $(P,\left\{\cdot,\cdot\right\})$ . If the Mercer kernel $K\in C_{b}^{3}(P\times P)$ is argumentwise invariant under the group $G$ then the kernel estimator $\widehat{h}_{\lambda,N}$ in Theorem 3.13 is $G$ -invariant. If additionally, the group action has a momentum map $\mathbf{J}:P\longrightarrow\mathfrak{g}^{\ast}$ associated, then its fibers are preserved by the flow $F_{t}$ of the corresponding Hamiltonian vector field $X_{\widehat{h}_{\lambda,N}}$ , that is,

\mathbf{J}\circ F_{t}=\mathbf{J}.

Proof.

By [Gins 12, Property 3.13], if the kernel is argumentwise invariant under $G$ -action, then the corresponding RKHS $\mathcal{H}_{K}$ consists exclusively of $G$ -invariant functions. This implies that the kernel estimator is $G$ -invariant since $\widehat{h}_{\lambda,N}\in\mathcal{H}_{null}^{\bot}\subset\mathcal{H}_{K}$ . The statement on the momentum conservation is a straightforward consequence of Noether’s Theorem (see, for instance, [Mars 99, Theorem 11.4.1]). ∎

4 Estimation and approximation error analysis

We now analyze the extent to which the structure-preserving kernel estimator $\widehat{h}_{\lambda,N}$ can recover the unknown Hamiltonian $H$ on a Poisson manifold $P$ . A standard approach is to decompose the reconstruction error $\widehat{h}_{\lambda,N}-H$ into the sum of what we shall call estimation and approximation errors, and further decompose the estimation error $\widehat{h}_{\lambda,N}-h_{\lambda}^{*}$ into what we shall call sampling error and noisy sampling error. More specifically,

	$\displaystyle\widehat{h}_{\lambda,N}-H$	$\displaystyle=\underbrace{\widehat{h}_{\lambda,N}-h_{\lambda}^{}}_{\text{Estimation error}}\quad+\underbrace{h_{\lambda}^{}-H}_{\text{Approximation error}},$		(4.1)
		$\displaystyle=\underbrace{\widetilde{h}_{\lambda,N}-h_{\lambda}^{}}_{\text{Sampling error}}\quad+\quad\underbrace{\widehat{h}_{\lambda,N}-\widetilde{h}_{\lambda,N}}_{\text{Noisy sampling error}}\quad+\quad\underbrace{h_{\lambda}^{}-H}_{\text{Approximation error}},$		(4.1)

where $h^{*}_{\lambda}\in\mathcal{H}_{K}$ is the best-in-class function introduced in (3.8) that minimizes the regularized statistical risk (3.7) and $\widetilde{h}_{\lambda,N}:=(Q_{N}+\lambda)^{-1}Q_{N}H$ is the noise-free part of the estimator $\widehat{h}_{\lambda,N}$ . Indeed, the noisy sampling error can be computed as

\displaystyle\widehat{h}_{\lambda,N}-\widetilde{h}_{\lambda,N}=\frac{1}{\sqrt{N}}(Q_{N}+\lambda)^{-1}A_{N}^{*}\mathbf{X}_{\sigma^{2},N}-(Q_{N}+\lambda)^{-1}Q_{N}H=\frac{1}{\sqrt{N}}(Q_{N}+\lambda)^{-1}A_{N}^{*}\mathbf{E}_{N}

where the noise vector $\mathbf{E}_{N}$ is

\mathbf{E}_{N}=\mathrm{Vec}\left(\bm{\varepsilon}^{(1)}|\cdots|\bm{\varepsilon}^{(N)}\right)\in T_{\mathbf{Z}_{N}}P,

which, by hypothesis, follows a multivariate distribution with zero mean and variance $\sigma^{2}I_{dN}$ .

Following the approach introduced in [Feng 23, Hu 24], we shall separately analyze these three errors. The estimation of the total error relies on both the operator representations in Section 3.3 and the differential kernel representations in Section 3.4 of the estimator. Since the operator representations of the estimator are similar to those in the Euclidean case, the approximation error and sampling error can be analyzed following an approach very close to the one in [Hu 24]. Therefore, to obtain the total error in (4.1), it is sufficient to analyze the noisy sampling error. As in the previous sections, we assume that $K\in C_{b}^{3}(P\times P)$ is a Mercer kernel. The proofs of the following results can be found in the appendix.

Lemma 4.1.

Let $K\in C_{b}^{3}(P\times P)$ be a Mercer kernel. For any function $h\in\mathcal{H}_{K}$ and $0<\delta<1$ , with probability at least $1-\delta$ , there holds

\displaystyle\|Q_{N}h-Qh\|_{\mathcal{H}_{K}}\leq\left(\sqrt{\frac{8\log(2/\delta)}{N}}+1\right)\sqrt{\frac{\log(2/\delta)}{N}}C^{2}\kappa^{2}\|h\|_{\mathcal{H}_{K}}.

Lemma 4.2 (Error bounds of noisy sampling error).

For any $\delta>0$ , with a probability of at least $1-\delta/2$ , it holds that

\displaystyle\left\|\widetilde{h}_{\lambda,N}-\widehat{h}_{\lambda,N}\right\|_{\mathcal{H}_{K}}\leq\frac{\sigma\kappa}{\lambda}\sqrt{\frac{1}{N}}\left(1+\sqrt{\frac{1}{c}\log(4/\delta)}\right),

where the constant $c$ is related to the Hanson-Wright inequality [Rude 13] (see the proof in Appendix A.5).

Assumption 4.3 (Source condition).

Assume that the unknown Hamiltonian function in the Poisson system (2.6) lies in so-called source space, that is,

\displaystyle H\in\Omega_{S}^{\gamma}:=\{h\in\mathcal{H}_{K}\mid h=Q^{\gamma}\psi,\psi\in\mathcal{H}_{K},\|\psi\|_{\mathcal{H}_{K}}<S\},

(4.2)

for some $\gamma\in(0,1)$ and $S>0$ .

Proposition 4.4.

If a Hamiltonian function $H$ satisfies the source condition (4.2), then $H\in\mathcal{H}_{\mathrm{null}}^{\bot}$ . In other words, $\Omega_{S}^{\gamma}\subset\mathcal{H}_{\mathrm{null}}^{\bot}$ .

Proof.

Recall first that by Proposition 3.8, the operator $Q=A^{*}A$ is a positive compact operator. Let $Q=\sum_{n=1}^{L}\lambda_{n}\langle\cdot,e_{n}\rangle e_{n}$ (possibly $L=\infty$ ) be the spectral decomposition of $Q$ with $0<\lambda_{n+1}<\lambda_{n}$ and $\{e_{n}\}_{n=1}^{L}$ be an orthonormal basis of $\mathcal{H}_{K}$ . Hence for any $\gamma\in(0,1)$ , we have $Q^{\gamma}=\sum_{n=1}^{L}\lambda_{n}^{\gamma}\langle\cdot,e_{n}\rangle_{\mathcal{H}_{K}}e_{n}$ , which implies that $Q^{\gamma}$ is adjoint. Notice that by the representation of the operator $Q$ given in (3.11), for an arbitrary Casimir function $h\in\mathcal{H}_{\mathrm{null}}$ (if it is not an empty set), we have that $Qh=0$ , which implies that $\lambda_{n}\langle h,e_{n}\rangle_{\mathcal{H}_{K}}=0$ for all $n=1,\ldots,L$ . Hence for any $\gamma\in(0,1)$ , we have $\lambda_{n}^{\gamma}\langle h,e_{n}\rangle_{\mathcal{H}_{K}}=0$ for all $n=1,\ldots,L$ . Then we obtain that $Q^{\gamma}h=0$ . Finally for arbitrary $\psi\in\mathcal{H}_{K}$ , we compute the inner product $\langle Q^{\gamma}\psi,h\rangle_{\mathcal{H}_{K}}=\langle\psi,Q^{\gamma}h\rangle_{\mathcal{H}_{K}}=0$ , which yields that $Q^{\gamma}\psi\in\mathcal{H}_{\mathrm{null}}^{\bot}$ for any $\psi\in\mathcal{H}_{K}$ . Therefore, we conclude that $\Omega_{S}^{\gamma}\subset\mathcal{H}_{\mathrm{null}}^{\bot}$ . ∎

Proposition 4.5.

Let $Q=\sum_{n=1}^{L}\lambda_{n}\langle\cdot,e_{n}\rangle e_{n}$ (possibly $L=\infty$ ) be a spectral decomposition of $Q$ , with an orthonormal basis $\{e_{n}\}_{n=1}^{L}$ with $0<\lambda_{n+1}\leq\lambda_{n}$ for all $n=1,\ldots,L-1$ . Then, $\Omega_{S}^{\gamma}$ has the following characterization

\displaystyle\Omega_{S}^{\gamma}=\mathcal{H}_{S}^{\gamma}:=\left\{h\in\mathcal{H}_{\mathrm{null}}^{\bot}\mid\sum_{n=1}^{L}\frac{1}{\lambda_{n}^{2\gamma}}|\langle h,e_{n}\rangle_{\mathcal{H}_{K}}|^{2}<S^{2}\right\}.

Proof.

First, we show that $\mathcal{H}_{S}^{\gamma}\subset\Omega_{S}^{\gamma}$ . For an arbitrary $h\in\mathcal{H}_{S}^{\gamma}$ , we can define a function $\psi\in\operatorname{span}\{e_{n},n\in\mathbb{N}_{+}\}$ such that $\langle\psi,e_{n}\rangle_{\mathcal{H}_{K}}=\frac{1}{\lambda_{n}^{\gamma}}\langle h,e_{n}\rangle_{\mathcal{H}_{K}}$ for all $n\in\mathbb{N}^{+}$ . One computes

\displaystyle\|\psi\|_{\mathcal{H}_{K}}^{2}=\sum_{n\in\mathbb{N}^{+}}\frac{1}{\lambda_{n}^{2\gamma}}|\langle h,e_{n}\rangle_{\mathcal{H}_{K}}|^{2}<S^{2}.

Furthermore, since $h\in\mathcal{H}_{\mathrm{null}}^{\bot}=\operatorname{span}\{e_{n},n\in\mathbb{N}_{+}\}$ , we have that

\displaystyle B^{\gamma}\psi=\sum_{n\in\mathbb{N}^{+}}\lambda_{n}^{\gamma}\langle\psi,e_{n}\rangle_{\mathcal{H}_{K}}e_{n}=\sum_{n\in\mathbb{N}^{+}}\lambda_{n}^{\gamma}\frac{1}{\lambda_{n}^{\gamma}}\langle h,e_{n}\rangle_{\mathcal{H}_{K}}e_{n}=h,

which yields that $h\in\Omega_{S}^{\gamma}$ . Now, we show that $\Omega_{S}^{\gamma}\subset\mathcal{H}_{S}^{\gamma}$ . For an arbitrary $h\in\Omega_{S}^{\gamma}$ , there exists a function $\psi\in\mathcal{H}_{K}$ with norm $\|\psi\|_{\mathcal{H}_{K}}<S$ , such that $B^{\gamma}\psi=h$ . Let $\psi^{\prime}$ be the projection of $\psi$ onto $\mathcal{H}_{\mathrm{null}}^{\bot}$ , then $h=B^{\gamma}\psi^{\prime}$ . Hence, we obtain that

\displaystyle\sum_{n\in\mathbb{N}^{+}}\lambda_{n}^{\gamma}\langle\psi^{\prime},e_{n}\rangle_{\mathcal{H}_{K}}e_{n}=\sum_{n\in\mathbb{N}^{+}}\langle h,e_{n}\rangle_{\mathcal{H}_{K}}e_{n},

which implies that $\langle\psi^{\prime},e_{n}\rangle_{\mathcal{H}_{K}}=\frac{1}{\lambda_{n}^{\gamma}}\langle h,e_{n}\rangle_{\mathcal{H}_{K}}$ for all $n\in\mathbb{N}^{+}$ since $\lambda_{n}>0$ . Therefore,

\displaystyle\|\psi^{\prime}\|_{\mathcal{H}_{K}}^{2}=\sum_{n\in\mathbb{N}^{+}}\frac{1}{\lambda_{n}^{2\gamma}}|\langle h,e_{n}\rangle_{\mathcal{H}_{K}}|^{2}<S^{2}.

Combing Proposition 4.4, we have that $h\in\mathcal{H}_{S}^{\gamma}$ . Therefore, we can conclude that $\Omega_{S}^{\gamma}=\mathcal{H}_{S}^{\gamma}$ . ∎

In Section 3.4, we showed that the kernel estimator $\widehat{h}_{\lambda,N}\in\mathcal{H}_{\mathrm{null}}^{\bot}$ . In this section, we will prove that for an arbitrary Hamiltonian function $H$ in the source space $\Omega_{S}^{\gamma}$ , the kernel estimator $\widehat{h}_{\lambda,N}$ converges to the Hamiltonian $H$ as $N$ tends to infinity.

Notice that the results in Lemma 4.2 are very similar to the error bounds of the noisy sampling error in [Hu 24]. Hence, combining the approximation error and the noisy sampling error bounds obtained in [Hu 24], we immediately formulate probably approximately correct (PAC) bounds for the reconstruction error in which the regularization constant $\lambda$ remains fixed and the size of the estimation sample is allowed to vary independently from it.

Theorem 4.6 (PAC bounds of the total reconstruction error).

Let $\widehat{h}_{\lambda,N}$ be the unique minimizer of the optimization problem (3.4). Suppose that Assumption 4.3 holds, that is, $H\in\Omega_{S}^{\gamma}$ as defined in (4.2). Then, for any $\varepsilon,\delta>0$ , there exist $\lambda>0$ and $n\in\mathbb{N}_{+}$ such that for all $N>n$ , it holds that

\displaystyle\mathbb{P}\left(\left\|\widehat{h}_{\lambda,N}-H\right\|_{\mathcal{H}_{K}}>\varepsilon\right)<\delta.

Notice that by Theorem 4.6, we can make the total reconstruction error arbitrarily small with the probability close to one with a fixed regularization parameter. However, this theorem provides no convergence rates. Similar to the approach followed in [Hu 24], we now consider an adaptive regularization parameter to obtain convergence rates. To this end, we shall assume that

\displaystyle\lambda\propto N^{-\alpha},\quad\alpha>0,

(4.3)

where the symbol $\propto$ means that $\lambda$ has the order $N^{-\alpha}$ as $N\to\infty$ . Then, combining Lemma 4.2 with [Hu 24, Theorem 4.7], we obtain convergence rates for the total reconstruction error when $\alpha\in(0,\frac{1}{3})$ . If we further suppose that the coercivity condition in the following definition holds, we can then improve the convergence integral to $\alpha\in(0,\frac{1}{2})$ .

Definition 4.7 (Coercivity condition).

In the same setup as in Theorem 3.13 we say that the coercivity condition is satisfied, if for all $h\in\mathcal{H}_{K}$ there exists $c_{\mathcal{H}_{K}}>0$ such that

\displaystyle\|Ah\|^{2}_{L^{2}(\mu_{\mathbf{Z}})}=\|X_{h}\|^{2}_{L^{2}(\mu_{\mathbf{Z}})}\geq c_{\mathcal{H}_{K}}\|h\|^{2}_{\mathcal{H}_{K}}.

(4.4)

The largest constant $c_{\mathcal{H}_{K}}$ that satisfies (4.4) is called the coercivity constant.

Theorem 4.8 (Convergence upper rates for the total reconstruction error).

Let $\widehat{h}_{\lambda,N}$ be the unique minimizer of the optimization problem (3.4). Suppose that $\lambda$ satisfies (4.3) and that $H$ satisfies the source condition (4.2), that is, $H\in\Omega_{S}^{\gamma}$ . Then for all $\alpha\in(0,\frac{1}{3})$ , and for any $0<\delta<1$ , with probability as least $1-\delta$ , it holds that

\displaystyle\left\|\widehat{h}_{\lambda,N}-H\right\|_{\mathcal{H}_{K}}\leq C(\gamma,\delta,\kappa)~{}N^{-\min\{\alpha\gamma,\frac{1}{2}(1-3\alpha)\}},

where

C(\gamma,\delta,\kappa)=\max\left\{\|B^{-\gamma}H\|_{\mathcal{H}_{K}},8\sqrt{4\log(8/\delta)}d^{\frac{3}{2}}\kappa^{3}\|H\|_{\mathcal{H}_{K}}\right\}.

Moreover, if the coercivity condition (4.4) holds, then for all $\alpha\in(0,\frac{1}{2})$ , and for any $0<\delta<1$ , and with probability as least $1-\delta$ , it holds that

\displaystyle\|\widehat{h}_{\lambda,N}-H\|_{\mathcal{H}_{K}}\leq C(\gamma,\delta,\sigma,\kappa,c,c_{\mathcal{H}_{K}})~{}N^{-\min\{\alpha\gamma,\frac{1}{2}(1-2\alpha)\}},

where

C(\gamma,\delta,\sigma,\kappa,c,c_{\mathcal{H}_{K}})\\ =\max\left\{\|B^{-\gamma}H\|_{\mathcal{H}_{K}},\sigma\kappa\sqrt{d}\left(1+\sqrt{\frac{1}{c}\log(4/\delta)}\right),4\sqrt{2\log(8/\delta)}d\kappa^{2}\left(2+\kappa\sqrt{\frac{d}{c_{\mathcal{H}_{K}}}}\right)\|H\|_{\mathcal{H}_{K}}\right\}.

Theorem 4.8 guarantees that the structure-preserving kernel estimator provides a function that is close to the data-generating Hamiltonian function with respect to the RKHS norm. As a consequence, we now prove that the flow of the learned Hamiltonian system will uniformly approximate that of the data-generating one.

In the following lemma and proposition, we take the Sasaki metric on the tangent buddle $TP$ induced by the Riemannian metric $g$ on the Poisson manifold $P$ . To define Sasaki metric on $TP$ , let $(x,v)\in TP$ and $V,W$ be tangent vectors to $TP$ at $(x,v)$ . Choose curves in $TP$

\displaystyle\alpha:t\mapsto(x(t),v(t)),\quad\beta:s\mapsto(y(t),w(t)),

with $x(0)=y(0)=x$ , $v(0)=w(0)=v$ and $V=\alpha^{\prime}(0)$ , $W=\beta^{\prime}(0)$ . Then the Sasaki metric is defined as

\displaystyle g_{1}(x,v)(V,W)=g(x)(d\pi(V),d\pi(W))+g(x)\left(\frac{D}{\mathrm{d}t}\bigg{|}_{t=0}v(t),\frac{D}{\mathrm{d}t}\bigg{|}_{t=0}w(t)\right),

(4.5)

where $\pi:TP\to P$ is the canonical projection, $\frac{D}{\mathrm{d}t}$ stands for the covariant derivative. As stated in Remark 2.14, we shall assume that this Sasaki metric (4.5) is complete. Then we obtain the following lemma.

Lemma 4.9.

Let $h\in C_{b}^{2}(P)$ . Then for any $f\in C^{1}(P)$ and $z_{1},z_{2}\in P$ , we have

\displaystyle X_{h}[f](z_{1})-X_{h}[f](z_{2})\leq\|h\|_{C_{b}^{2}}d(z_{1},z_{2}),

where $d$ is the geodesic distance on the Riemannian manifold $(P,g)$ .

Proof.

Since $h\in C_{b}^{2}(P)$ , we have that $Dh\in C_{b}^{1}(TP)$ . Then by the Fundamental Theorem of Calculus (2.10), we obtain that

\displaystyle Dh(y_{1})-Dh(y_{2})=\int_{0}^{1}D^{2}h(\gamma(t))\cdot\gamma^{\prime}(t)\mathrm{d}t\leq\|D^{2}h\|_{\infty}\int_{0}^{1}\|\gamma^{\prime}(t)\|_{g_{1}}\mathrm{d}t,

where $\gamma:[0,1]\to TP$ is a differentiable curve connecting $\gamma(0)=y_{1}$ and $\gamma(1)=y_{2}$ in the Riemannian manifold $(TP,g_{1})$ with Sasaki metric $g_{1}$ induced by the Riemannian metric $g$ on $P$ . Denote $\gamma(t)=(z(t),v(t))$ with $z(t)\in P$ and $v(t)\in T_{z(t)}P$ for all $t\in[0,1]$ . For each $t\in[0,1]$ , let $\alpha_{t}:[0,1]\to TP$ be a smooth curve such that $\alpha_{t}(0)=(z(t),v(t))$ and $\alpha^{\prime}_{t}(0)=(z^{\prime}(t),v^{\prime}(t))$ . Denote $\alpha_{t}=(a_{t},b_{t})$ with $a_{t}(s)\in P$ and $b_{t}(s)\in T_{a_{t}(s)}P$ for all $s\in[0,1]$ . Notice that

\displaystyle\frac{D}{\mathrm{d}s}\bigg{|}_{s=0}b_{t}(s)=\nabla_{\dot{a}_{t}(0)}b_{t}(0)=\nabla_{z^{\prime}(t)}v(t).

Then, by the definition of the Sasaki metric $g_{1}$ in (4.5), we obtain

\displaystyle\|\gamma^{\prime}(t)\|_{g_{1}}=\|\dot{a}_{t}(0)\|_{g}+\left\|\frac{D}{\mathrm{d}s}\bigg{|}_{s=0}b_{t}(s)\right\|_{g}=\|z^{\prime}(t)\|_{g}+\left\|\nabla_{z^{\prime}(t)}v(t)\right\|_{g}.

Let $y_{1}=(z_{1},X_{f}(z_{1}))$ and $y_{2}=(z_{2},X_{f}(z_{2}))$ for $z_{1},z_{2}\in P$ . Choose the $\{z(t),t\in[0,1]\}$ be the geodesic curve in $P$ . Let $\{v(t),t\in[0,1]\}$ be a curve which is parallel along the curve $\{z^{\prime}(t):t\in[0,1]\}$ , i.e., $\nabla_{z^{\prime}(t)}v(t)=0$ for all $t\in[0,1]$ . Then we obtain that $\|\gamma^{\prime}(t)\|_{g_{1}}=\|z^{\prime}(t)\|_{g}$ . Therefore,

\displaystyle X_{h}[f](z_{1})-X_{h}[f](z_{2})=Dh(y_{1})-Dh(y_{2})\leq\|D^{2}h\|_{\infty}\int_{0}^{1}\|z^{\prime}(t)\|_{g}\mathrm{d}t\leq\|h\|_{C_{b}^{2}}d(z_{1},z_{2}),

where $d$ is the geodesic distance on the Riemannian manifold $(P,g)$ . The result follows. ∎

Proposition 4.10 (From discrete data to continuous-time flows).

Let $\widehat{H}=\widehat{h}_{\lambda,N}$ be the structure-preserving kernel estimator of $H$ using a kernel $K\in C_{b}^{5}(P\times P)$ . Let $F:[0,T]\times P\longrightarrow P$ and $\widehat{F}:[0,T]\times P\longrightarrow P$ be the flows over the time interval $[0,T]$ of the Hamilton equations associated to the Hamiltonian functions $H$ and $\widehat{H}$ , respectively. Suppose that the geodesic distance $d(\cdot,y)$ of the Poisson manifold $P$ is in $C_{b}^{1}(P)$ for each $y\in P$ and that the $C_{b}^{1}$ norm of the family $\{d(\cdot,y),y\in p\}$ is uniformly bounded, that is, $C_{1}=\sup_{y\in P}\|d(\cdot,y)\|_{C_{b}^{1}}<\infty$ . Then, for any initial condition $z\in P$ , we have that

\displaystyle d^{\infty}\left(F({z}),\widehat{F}({z})\right):=\max_{t\in[0,T]}d\left(F_{t}({z}),\widehat{F}_{t}({z})\right)\leq\widetilde{C}\|H-\widehat{H}\|_{\mathcal{H}_{K}},

with the constant $\widetilde{C}=\kappa CC_{1}T\exp\{\|\widehat{H}\|_{C_{b}^{2}}T\}$ , where $C$ is the constant introduced in Assumption 3.7.

Proof.

Let $f\in C_{b}^{1}(P)$ . For each $t\in[0,T]$ and ${z}\in P$ , we have

\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}(f\circ F_{t})(z)=\mathbf{d}f(F_{t}(z))\cdot X_{H}(F_{t}(z))=X_{H}[f](F_{t}(z)).

By the equation (A.4), we have that for any $z\in P$ ,

\displaystyle X_{H}[f](z)-X_{\widehat{H}}[f](z)

\displaystyle=Df(z)\cdot(X_{H}(z)-X_{\widehat{H}}(z))\leq\|f\|_{C_{b}^{1}}\|X_{H-\widehat{H}}(z)\|_{g}\leq\kappa C\|f\|_{C_{b}^{1}}\|H-\widehat{H}\|_{\mathcal{H}_{K}}.

Since $K\in C_{b}^{5}(P\times P)$ , Theorem 3.3 implies that $\widehat{H}\in C_{b}^{2}(P)$ . Then by Lemma 4.9, for any $z_{1},z_{2}\in P$ ,

\displaystyle X_{\widehat{H}}[f](z_{1})-X_{\widehat{H}}[f](z_{2})\leq\|\widehat{H}\|_{C_{b}^{2}}d(z_{1},z_{2}).

Furthermore, for each $t\in[0,T]$ , we can choose $y_{t}\in P$ such that the point $\widehat{F}_{t}({z})$ is in the geodesic curve connecting $F_{t}({z})$ and $y_{t}$ . Now, by hypothesis, $f:=d(\cdot,y_{t})\in C_{b}^{1}(P)$ . Therefore, we obtain that

	$\displaystyle d\left(F_{t}({z}),\widehat{F}_{t}({z})\right)$	$\displaystyle=f(F_{t}(z))-f(\widehat{F}_{t}(z))=\int_{0}^{t}X_{H}[f](F_{s}({z}))-X_{\widehat{H}}[f](\widehat{F}_{s}({z}))\mathrm{d}s$
		$\displaystyle\leq\int_{0}^{t}\left\|X_{H}[f](F_{s}({z}))-X_{\widehat{H}}[f](F_{s}({z}))\right\|+\left\|X_{\widehat{H}}[f](F_{s}({z}))-X_{\widehat{H}}[f](\widehat{F}_{s}({z}))\right\|\mathrm{d}s$
		$\displaystyle\leq\kappa C\\|f\\|_{C_{b}^{1}}\\|H-\widehat{H}\\|_{\mathcal{H}_{K}}t+\\|\widehat{H}\\|_{C_{b}^{2}}\int_{0}^{t}d\left(F_{s}(z),\widehat{F}_{s}(z)\right)\mathrm{d}s.$

Then by the integral form of Grönwall’s inequality, for $d\left(F_{t}({z}),\widehat{F}_{t}({z})\right)$ , we obtain

\displaystyle d\left(F_{t}({z}),\widehat{F}_{t}({z})\right)\leq\kappa C\|d(\cdot,y_{t})\|_{C_{b}^{1}}\|H-\widehat{H}\|_{\mathcal{H}_{K}}t\exp\left\{\|\widehat{H}\|_{C_{b}^{2}}t\right\}.

Therefore, we obtain

\displaystyle d^{\infty}\left(F({z}),\widehat{F}({z})\right):=\max_{t\in[0,T]}d\left(F_{t}({z}),\widehat{F}_{t}({z})\right)\leq\widetilde{C}\|H-\widehat{H}\|_{\mathcal{H}_{K}},

where $\widetilde{C}=\kappa CC_{1}T\exp\{\|\widehat{H}\|_{C_{b}^{2}}T\}$ . The result follows. ∎

5 Test examples and numerical experiments

This section spells out the global structure-preserving kernel estimator and illustrates its performance on two important instances of Poisson systems, namely, Hamiltonian systems on non-Euclidean symplectic manifolds and Lie-Poisson systems on the dual of a Lie algebra. More precisely, we shall derive explicit expressions for the global estimator $\widehat{h}_{\lambda,N}$ in (3.16) in these two scenarios.

5.1 Learning of Lie-Poisson systems

The Lie-Poisson systems introduced in Example 2.3 are very important examples of Poisson systems. In this subsection, we shall make explicit the estimator $\widehat{h}_{\lambda,N}$ in this case. Let $\mathfrak{g}$ be a Lie algebra and $\mathfrak{g}^{*}$ be its dual equipped with the Lie-Poisson bracket $\{\cdot,\cdot\}_{+}$ defined in (2.2). A straightforward computation shows that the Lie-Poisson dynamical system associated with a Hamiltonian function $H:\mathfrak{g}^{*}\to\mathbb{R}$ is

\displaystyle\dot{\mu}=X_{H}(\mu)=\operatorname{ad}^{*}_{\frac{\delta H}{\delta\mu}}\mu,

(5.1)

where $\delta H/\delta\mu\in\mathfrak{g}$ is the functional derivative introduced in (2.3). Since $\mathfrak{g}$ and $\mathfrak{g}^{\ast}$ are vector spaces, natural Riemannian metrics can be obtained out of a non-degenerate inner product $\left\langle\cdot,\cdot\right\rangle_{\mathfrak{g}}:\mathfrak{g}\times\mathfrak{g}\rightarrow\mathbb{R}$ on $\mathfrak{g}$ . The non-degeneracy of $\left\langle\cdot,\cdot\right\rangle_{\mathfrak{g}}$ allows us to define musical isomorphisms $\flat:\mathfrak{g}\rightarrow\mathfrak{g}^{\ast}$ and $\sharp:\mathfrak{g}^{\ast}\rightarrow\mathfrak{g}$ by $\xi^{\flat}:=\left\langle\xi,\cdot\right\rangle_{\mathfrak{g}}$ , $\xi\in\mathfrak{g}$ , and $\sharp:=\flat^{-1}$ , as well as a natural non-degenerate inner product $\left\langle\cdot,\cdot\right\rangle_{\mathfrak{g}^{\ast}}:\mathfrak{g}^{\ast}\times\mathfrak{g}^{\ast}\rightarrow\mathbb{R}$ on $\mathfrak{g}^{\ast}$ given by $\left\langle\xi^{\flat},\eta^{\flat}\right\rangle_{\mathfrak{g}^{\ast}}=\left\langle\xi,\eta\right\rangle_{\mathfrak{g}}$ , $\xi,\eta\in\mathfrak{g}$ . Notice that

\left\langle\xi^{\flat},\eta\right\rangle=\left\langle\left\langle\xi,\cdot\right\rangle_{\mathfrak{g}},\eta\right\rangle=\left\langle\xi,\eta\right\rangle_{\mathfrak{g}}=\left\langle\xi^{\flat},\eta^{\flat}\right\rangle_{\mathfrak{g}^{\ast}}\quad\mbox{for all $\xi,\eta\in\mathfrak{g}$.}

(5.2)

Using these relations and (2.3), we note that

\left\langle\nabla H(\mu),\nu\right\rangle_{\mathfrak{g}^{\ast}}=DH(\mu)\cdot\nu=\left\langle\nu,\frac{\delta H}{\delta\mu}\right\rangle=\left\langle\nu,\frac{\delta H}{\delta\mu}^{\flat}\right\rangle_{\mathfrak{g}^{\ast}},\mbox{ for all $\nu,\mu\in\mathfrak{g}^{\ast}$.}

(5.3)

Consequently,

\nabla H(\mu)=\frac{\delta H}{\delta\mu}^{\flat}=\left\langle\frac{\delta H}{\delta\mu},\cdot\right\rangle_{\mathfrak{g}}\quad\mbox{and hence $\frac{\delta H}{\delta\mu}=\nabla H(\mu)^{\sharp}$.}

(5.4)

We now compute the compatible structure $J$ defined in (2.12) associated with the Lie-Poisson systems (5.1). Note first that, for any $\mu\in\mathfrak{g}^{*}$ and $\eta\in\mathfrak{g}$ , by (5.1) and (5.4), we have that

\displaystyle J(\mu)\nabla H(\mu)=X_{H}(\mu)=\operatorname{ad}^{*}_{\frac{\delta H}{\delta\mu}}\mu=\operatorname{ad}^{*}_{\nabla H(\mu)^{\sharp}}\mu.

As this can be done for any Hamiltonian function and their gradients span $\mathfrak{g}^{\ast}$ , this yields

\displaystyle J(\mu)\nu=\operatorname{ad}^{*}_{\nu^{\sharp}}\mu,\quad\text{ for all }\nu\in{\mathfrak{g}}^{\ast}.

(5.5)

To compute the estimator to $\widehat{h}_{\lambda,N}$ one first evaluates $J$ explicitly using (5.5), we then set

\mathbb{J}_{N}=\operatorname{diag}\{J(\mathbf{Z}^{(1)}),\ldots,J(\mathbf{Z}^{(N)})\},

and substitute it into (3.16). More explicitly, note that the generalized differential Gram matrix $G_{N}$ defined in (3.15) reduces in this case to

\displaystyle G_{N}\mathbf{c}=X_{g_{N}(\mathbf{c},X_{K_{\cdot}}(\mathbf{Z}_{N}))}(\mathbf{Z}_{N})=X_{\mathbf{c}^{\top}\mathbb{J}_{N}\nabla_{1}K(\mathbf{Z}_{N},\cdot)}(\mathbf{Z}_{N})=\mathbb{J}_{N}\nabla_{1,2}K(\mathbf{Z}_{N},\mathbf{Z}_{N})\mathbb{J}_{N}^{\top}\mathbf{c}.

for all $\mathbf{c}\in T_{\mathbf{Z}_{N}}\mathfrak{g}^{*}$ . Therefore, the estimator in Theorem 3.13 turns to

	$\displaystyle\widehat{h}_{\lambda,N}$	$\displaystyle=g_{N}((G_{N}+\lambda NI)^{-1}\mathbf{X}_{\sigma^{2},N},X_{K_{\cdot}}(\mathbf{Z}_{N},\cdot))$
		$\displaystyle=\mathbf{X}_{\sigma^{2},N}^{\top}\left(\mathbb{J}_{N}\nabla_{1,2}K(\mathbf{Z}_{N},\mathbf{Z}_{N})\mathbb{J}_{N}^{\top}+\lambda NI\right)^{-1}\mathbb{J}_{N}\nabla_{1}K(\mathbf{Z}_{N},\cdot).$		(5.6)

In what follows, we demonstrate the effectiveness of this learning scheme using two examples. The first one is the three-dimensional classical rigid body, which models the motion of solid bodies without deformation. The second example is the dynamics of the underwater vehicle, which has a nine-dimensional underlying phase space. We emphasize that the state spaces of Lie-Poisson systems are duals of Lie algebras and, hence, are isomorphic to Euclidean spaces; in such a framework, the Gaussian kernel is a natural choice of universal kernels. As such, even when $H$ does not belong to $\mathcal{H}_{K}$ , we can hope for learning a proxy function in $\mathcal{H}_{K}$ that is close to a ground-truth Hamiltonian $H$ .

We will see in these two examples that $\widehat{h}_{\lambda,N}$ can only recover the ground-truth Hamiltonian function up to a Casimir function. We emphasize that the estimator $\widehat{h}_{\lambda,N}$ is always unique, and the possible difference by a Casimir function is present because $\widehat{h}_{\lambda,N}\in\mathcal{H}^{\bot}_{\mathrm{null}}$ , while the ground-truth $H$ may not belong to $\mathcal{H}^{\bot}_{\mathrm{null}}$ , and hence the difference between $H$ and $\widehat{h}_{\lambda,N}$ could converge not to zero, but to a Casimir function. An analytic reason for this is that the Hamiltonian functions that we shall be using for data generation are quadratic and violate the source condition (4.2) (see Example 2.9 in [Hu 24] for a characterization of the RKHS of the Gaussian kernel). To better illustrate this point, we provide a third example in which the Hamiltonian $H\in\mathcal{H}^{\bot}_{\mathrm{null}}$ satisfies the source condition (4.2) for sufficiently large $S>0$ ( $S$ is the constant introduced in (4.2); see also Proposition 4.5). In that case, according to Theorem 4.8, $\widehat{h}_{\lambda,N}$ converges to $H$ in the $\mathcal{H}_{K}$ norm, and the difference converges to $0$ as $N\rightarrow\infty$ . We will see that, numerically, $\widehat{h}_{\lambda,N}$ accurately recovers the ground-truth Hamiltonian.

5.1.1 The rigid body

Rigid body motion has been introduced in Example 2.5. Using (5.5), the compatible structure $J$ associated with the Euclidean metric in ${\mathbb{R}}^{3}$ can be written in this case as

\displaystyle J(\Pi)=\hat{\Pi}.

Numerical results We pick $\mathbb{I}=\operatorname{diag}\{1,10,0.1\}$ and consider the standard Gaussian kernel with parameter $\eta$ on $\mathfrak{so}(3)^{*}\cong\mathbb{R}^{3}$ . In order to generate the data, we set $N=500$ and sample $N$ states from the uniform distribution on the cube $[-1,1]^{3}$ . We then obtain the Hamiltonian vector fields at these states that, this time, we do not perturb with noise ( $\sigma=0$ ). In view of (4.3), we set the regularization parameter $\lambda=c\cdot N^{-\alpha}$ with $\alpha=0.4$ . We then perform a grid search for the parameters $\eta$ and $c$ using $5$ -fold cross-validation, which results in the optimal parameters $\eta=2.5$ and $c=2.5\cdot 10^{-5}$ . Given that $\mathfrak{so}(3)^{*}\cong\mathbb{R}^{3}$ is $3$ -dimensional, we fix the third dimension, that is, $\Pi_{3}=0$ , and visualize both the ground-truth Hamiltonian functions $H$ and the reconstructed Hamiltonian functions $\widehat{h}_{\lambda,N}$ in the dimensions $(\Pi_{1},\Pi_{2})$ , see Figure 5.1 (a) and (b). The ground-truth Hamiltonian $H$ is quadratic and violates the source condition (4.2), so $\widehat{h}_{\lambda,N}$ and $H$ could differ by a Casimir function.

From Example 2.5, all possible Casimir functions are of the form $\Phi(\|\Pi\|^{2})$ for some function $\Phi:\mathbb{R}\rightarrow\mathbb{R}$ . We try to compensate for the difference via a simple Casimir of the form $a\cdot\|\Pi\|^{2}$ , for some constant $a\in\mathbb{R}$ . After trial and error, we select $a=1.85$ and visualize the “Casimir-corrected” estimator, that is $\widehat{h}_{\lambda,N}+1.85\cdot\|\Pi\|^{2}$ , see Figure 5.1 (c). We also visualize the mean-square error of the predicted Hamiltonian vector fields in a heatmap, see Figure 5.1 (d).

Refer to caption — Figure 5.1: Rigid body dynamics: (a) Ground-truth Hamiltonian (b) Learned Hamiltonian with $N=500$ (c) Learned Hamiltonian adjusted by a Casimir function (d) Squared error of the predicted Hamiltonian vector field

Analysis In this setting, the Hamiltonian function $H$ is a quadratic function, which, according to Example 2.9 in [Hu 24], does not belong to $\mathcal{H}_{K}$ . In particular, the source condition (4.2) cannot be satisfied. Therefore, even though the Gaussian kernel is universal, the proxy function in $\mathcal{H}_{K}$ , which approximates the ground-truth Hamiltonian $H$ , can only be learned at best up to a Casimir function. That is why we choose to manually correct the estimator by adding a Casimir function. From the numerical results, we see that without the “Casimir correction”, $H$ and $\widehat{h}_{\lambda,N}$ can be very different, despite the fact that the Hamiltonian vector fields are very well replicated.

5.1.2 Underwater vehicle dynamics

The Lie-Poisson underwater vehicle dynamics has been introduced in Example 2.6. Using also a Euclidean metric, (5.5), the compatible structure $J$ is

\displaystyle J(\Pi,\mathrm{Q},\Gamma)=\left(\begin{array}[]{ccc}\hat{\Pi}&\hat{\mathrm{Q}}&\hat{\Gamma}\\ \hat{\mathrm{Q}}&0&0\\ \hat{\Gamma}&0&0\end{array}\right).

Numerical results We pick $m=1$ , $g=9.8$ , $I=\operatorname{diag}\{1,2,3\}$ , $M=\operatorname{diag}\{3,2,1\}$ and $r_{G}=(1,1,1)^{\top}$ . We consider the standard Gaussian kernel with parameter $\eta$ on $\mathfrak{so}(3)^{*}\times{\mathbb{R}^{3}}^{*}\times{\mathbb{R}^{3}}^{*}\cong\mathbb{R}^{9}$ . We set $N=400$ and sample $N$ states from the uniform distribution on the cube $[-1,1]^{9}$ , and then obtain the Hamiltonian vector fields at these states. We set the regularization parameter $\lambda=c\cdot N^{-\alpha}$ with $\alpha=0.4$ . We then perform a grid search of the parameters $\eta$ and $c$ using $5$ -fold cross-validation, which results in the optimal parameter $\eta=5$ and $c=7.5\cdot 10^{-6}$ . Given that $\mathfrak{so}(3)^{*}\times{\mathbb{R}^{3}}^{*}\times{\mathbb{R}^{3}}^{*}$ is $9$ -dimensional, we fix the 3rd to 9th dimensions, that is, $\Pi_{3}=0$ , $Q=\Gamma={\bf 0}$ , and visualize both the ground-truth Hamiltonian functions $H$ and the reconstructed Hamiltonian functions $\widehat{h}_{\lambda,N}$ in the dimensions $(\Pi_{1},\Pi_{2})$ , see Figure 5.2 (a) and (b). The ground-truth Hamiltonian $H$ does not satisfy the source condition (4.2), so $\widehat{h}_{\lambda,N}$ and $H$ generically differ by a Casimir function. We also visualize the mean-square error of the predicted Hamiltonian vector fields with a heatmap, see Figure 5.2 (c).

Analysis In this setting, the Hamiltonian function $H$ is a quadratic function, which, according to Example 2.9 in [Hu 24], does not belong to $\mathcal{H}_{K}$ . In particular, the source condition (4.2) cannot be satisfied. Therefore, even though the Gaussian kernel is universal, the proxy function in $\mathcal{H}_{K}$ , which approximates the ground-truth Hamiltonian $H$ , can only be learned at best up to a Casimir function. From Example 2.6, all possible Casimir functions are of the form $\Phi\circ(C_{1},\ldots,C_{6})$ for some function $\Phi:\mathbb{R}^{6}\rightarrow\mathbb{R}$ . To determine such a function $\Phi$ , one might consider running a polynomial regression to approximate $\Phi$ . From the $H$ -axis, we see that, although the shapes of the Hamiltonian look alike, the exact Hamiltonian value still differs significantly by a constant, which is a trivial Casimir function. On the other hand, we notice that the Hamiltonian vector fields are very well replicated.

5.1.3 Exact recovery of Hamiltonians in the RKHS

In the rigid body and the underwater vehicle examples, we saw that even though the Hamiltonian vector field can be learned very well, the data-generating Hamiltonian function can only be estimated up to a Casimir function due to the violation of the source condition. We now aim to manually pick a ground-truth Hamiltonian that satisfies the source condition and reconstruct the Hamiltonian function accurately from the data. To find such a Hamiltonian, we use [Minh 10], which provides an orthonormal basis of the RKHS of the Gaussian kernel. It turns out that exactly one of the basis functions is a Casimir and any other basis function belongs to $\mathcal{H}^{\bot}_{\mathrm{null}}$ according to Example 2.5. Moreover, by Proposition 4.5, the basis functions, except the Casimir basis element, satisfies the source condition for sufficiently large $S>0$ . In view of the above, we consider the standard Gaussian kernel on $\mathbb{R}^{3}$ with fixed parameter $\eta=2$ and choose the Hamiltonian function as

H({\bf x})=x_{1}^{2}x_{2}^{3}e^{-\frac{\|{\bf x}\|^{2}}{\eta^{2}}}.

(5.7)

The Lie-Poisson structure is the same as the one we used for the rigid body, with the only difference in the experiment being the Hamiltonian function.

Numerical results To guarantee the matching of the RKHS, we consider the standard Gaussian kernel with the fixed parameter $\eta=2$ (the same as the ground-truth) on $\mathbb{R}^{3}$ . We set $N=500$ and sample $N$ states from the uniform distribution on the cube $[-1,1]^{3}$ , and then obtain the Hamiltonian vector fields at these states. In view of (4.3), we set the regularization parameter $\lambda=c\cdot N^{-\alpha}$ with $\alpha=0.4$ . We then perform a grid search of the parameter $c$ using $5$ -fold cross-validation, which results in the optimal parameter $c=7.5\cdot 10^{-6}$ . Given that $\mathfrak{so}(3)^{*}\cong\mathbb{R}^{3}$ is $3$ -dimensional, we fix the third dimension, that is, $x_{3}=0$ , and visualize both the ground-truth Hamiltonian functions $H$ and the reconstructed Hamiltonian functions $\widehat{h}_{\lambda,N}$ in the dimensions $(x_{1},x_{2})$ , see Figure 5.3 (a) and (b).

Analysis In this setting, the Hamiltonian function $H$ is chosen to satisfy the source condition (4.2). In this case, both the ground-truth Hamiltonian $H$ and the estimator $\widehat{h}_{\lambda,N}$ belong to $\mathcal{H}^{\bot}_{\mathrm{null}}$ . The possibility of the presence of Casimir functions is eliminated, and hence by Theorem 4.8, $\widehat{h}_{\lambda,N}$ converges to $H$ in the RKHS norm. As shown in Figure 5.3 (c), the Hamiltonian function can replicated even without “Casimir correction”, in contrast with the rigid body example.

5.2 Learning on non-Euclidean symplectic manifolds

As mentioned in Example 2.2, any symplectic manifold is a Poisson manifold. Let $M$ be a symplectic manifold with a symplectic form $\omega$ . Then the Hamiltonian vector field associated with Hamiltonian function $H$ defined in (2.5) is

\displaystyle X_{H}=\omega^{\sharp}(\bm{\mathrm{d}}H)=J\nabla H,

where $\omega^{\sharp}:T^{*}M\to TM$ is the bundle isomorphism induced by the symplectic form $\omega$ and $J:TM\to TM$ given by

\displaystyle J(x)v:=\omega_{x}^{\sharp}\left(g_{x}^{\flat}(v)\right),\quad\forall x\in M,v\in T_{x}M,

(5.8)

with the bundle isomorphism $g^{\flat}:TM\to T^{*}M$ induced by the Riemannian metric $g$ . Unlike in the more general Poisson setting, $J$ is a bundle isomorphism in the symplectic setup. We aim to learn the Hamiltonian function $H$ out of the equation

\displaystyle\dot{z}=J\nabla H(z).

(5.9)

Following Theorem 3.13, the kernel estimator in (3.16) is

\widehat{h}_{\lambda,N}=g_{N}(\widehat{\bf c},X_{K_{\cdot}}(\mathbf{Z}_{N},\cdot))=g_{N}(\widehat{\bf c},\mathbb{J}_{N}\nabla K(\mathbf{Z}_{N},\cdot)),

where $\widehat{\mathbf{c}}\in T_{\mathbf{Z}_{N}}M$ is given by

\displaystyle\widehat{\mathbf{c}}=(G_{N}+\lambda NI)^{-1}\mathbf{X}_{\sigma^{2},N}.

In applications, it is important to compute the kernel estimator $\widehat{h}_{\lambda,N}$ in pre-determined local patches that cover the manifold $M$ . By the local representations (2.17) and (3.21), the compatible structure $J$ and the generalized Gram matrix $G_{N}$ of the Hamiltonian system (5.9) are

\displaystyle J^{ij}=-\omega^{ik}g_{kj},\quad G_{N}^{i,j}=\omega(\mathbf{Z}^{(i)})\partial_{1,2}K(\mathbf{Z}^{(i)},\mathbf{Z}^{(j)})\omega^{\top}(\mathbf{Z}^{(j)})g^{(j)}.

In what follows, we demonstrate the effectiveness of this learning scheme on two examples that share the same symplectic manifold $M=S^{2}\times S^{2}$ with the same symplectic form $\omega$ and the same Riemannian metric $g$ , while governed by two different Hamiltonian functions $H$ . The first Hamiltonian function is the restricted Euclidean 3-norm to the manifold $M$ , and the second one describes the two-vortex dynamics. We will see that the first Hamiltonian is “well-behaved” and the second Hamiltonian exhibits singularities.

We emphasize that, technically, there are three challenges in relation to learning these systems. First, it is difficult to verify that a given Hamiltonian function belongs to the RKHS $\mathcal{H}_{K}$ of a pre-determined kernel $K$ on the manifold due to the lack of results in the literature on the characterization of functions in this type of RKHS on manifolds. Second, in case the Hamiltonian function does not lie in the RKHS of the selected kernel $K$ , it is important to investigate the universality properties of kernels on manifolds to hope for at least a good approximation. Third, the Hamiltonian function of the two-vortex system exhibits singularities, which pose even more difficulties to recover $H$ from $X_{H}$ , since the vector field diverges at singularities. We will see that even in this example, $\widehat{h}_{\lambda,N}$ still recovers the qualitative behaviors of the Hamiltonian function.

5.2.1 Metric, symplectic, and compatible structures in $S^{2}\times S^{2}$

Let $\overline{\mathcal{P}}=S^{2}\times S^{2}$ , the product of $2$ copies of unit $2$ -spheres in $\mathbb{R}^{3}$ . The phase space for the $2$ -vortex problem is $\mathcal{P}=\overline{\mathcal{P}}\backslash\Delta$ , where $\Delta$ is the diagonal, that is,

\displaystyle\Delta=\left\{\mathbf{x}=(x_{1},x_{2})\mid x_{1}=x_{2}\right\}.

The symplectic structure on $\mathcal{P}$ is given by

\displaystyle\omega_{\mathcal{P}}=\lambda_{1}\pi_{1}^{*}\omega_{S^{2}}+\lambda_{2}\pi_{2}^{*}\omega_{S^{2}},

(5.10)

where $\pi_{j}$ is the Cartesian projection on to the $j$ -th factor, $\omega_{S^{2}}$ the natural symplectic form on $S^{2}$ , and $\lambda_{j}$ the vorticity of the $j$ th vortex. The Poisson structure is given by

\{f,g\}=\lambda_{1}^{-1}\left[\mathrm{~{}d}_{1}f,\mathrm{~{}d}_{1}g,x_{1}\right]+\lambda_{2}^{-1}\left[\mathrm{~{}d}_{2}f,\mathrm{~{}d}_{2}g,x_{2}\right],

where $\left[\mathrm{d}_{j}f,\mathrm{~{}d}_{j}g,x_{j}\right]$ is the triple product of $[a,b,c]=a\cdot(b\times c)$ for $a,b,c\in\mathbb{R}^{3}$ .

Local coordinates. We first consider the local structure of the unit sphere $S^{2}$ as a Riemannian manifold. We parametrize the unit sphere $S^{2}$ in terms of two angles, that is, $\theta$ (the colatitude) and $\varphi$ (the longitude), and use two patches $(V_{1},\phi_{1})$ and $(V_{2},\phi_{2})$ to cover the whole sphere $S^{2}$ . The first patch is given by

\begin{cases}V_{1}=\{(\theta,\varphi)\mid\theta\in(0,\pi),\varphi\in(0,2\pi)\},\\ \phi_{1}:(\theta,\varphi)\mapsto(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta),\end{cases}

which implies that the semicircle from the north pole to the south pole lying in the $xz$ -plane is not covered by this chart. In order to cover the whole sphere, we need another chart obtained by deleting the semicircle in the $xy$ -plane from $(0,1,0)$ to $(0,1,0)$ and with $x\leq 0$ . In view that the patch $(V_{1},\phi_{1})$ is already topologically dense in $S^{2}$ , we omit the explicit parameterization of $(V_{2},\phi_{2})$ for convenience and focus on the first patch.

The unit sphere $S^{2}$ admits a natural symplectic form given by its area form that is locally expressed by $\omega_{S^{2}}=\sin\theta\mathrm{d}\theta\wedge\mathrm{d}\varphi$ and a Riemannian metric $g_{S^{2}}=\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\varphi^{2}$ inherited from the ambient Euclidean metric in ${\mathbb{R}}^{3}$ . The corresponding matrix forms in the coordinate patch $(V_{1},\phi_{1})$ can be expressed by

\omega_{S^{2}}(\theta,\varphi)=\begin{bmatrix}0&\sin\theta\\ -\sin\theta&0\end{bmatrix},\quad\text{ and }\quad g_{S^{2}}(\theta,\varphi)=\begin{bmatrix}1&0\\ 0&\sin^{2}\theta\end{bmatrix}.

These structures induce local symplectic and Riemannian metrics on the patch $(V_{1}\times V_{1},\phi_{1}\times\phi_{1})$ (that is, two copies of $(V_{1},\phi_{1})$ ) of the product manifold $S^{2}\times S^{2}$ given by

\displaystyle\omega_{\mathcal{P}}(\theta_{1},\varphi_{1},\theta_{2},\varphi_{2})=\operatorname{diag}\{\lambda_{1}\omega_{S^{2}}(\theta_{1},\varphi_{1}),\lambda_{N}\omega_{S^{2}}(\theta_{2},\varphi_{2})\},

where $(\theta_{1},\varphi_{1},\theta_{2},\varphi_{2})\in V_{1}\times V_{1}$ . The matrix form of the product metric $g_{\mathcal{P}}$ is

\displaystyle g_{\mathcal{P}}(\theta_{1},\varphi_{1},\theta_{2},\varphi_{2})=\operatorname{diag}\{g_{S^{2}}(\theta_{1},\varphi_{1}),g_{S^{2}}(\theta_{2},\varphi_{2})\},

hence, the compatible structure $J_{\mathcal{P}}$ of $2$ -vortex system is

\displaystyle J_{\mathcal{P}}(\theta_{1},\varphi_{1},\theta_{2},\varphi_{2})=\operatorname{diag}\{\lambda_{1}\omega_{S^{2}}(\theta_{1},\varphi_{1})g_{S^{2}}(\theta_{1},\varphi_{1}),\lambda_{2}\omega_{S^{2}}(\theta_{2},\varphi_{2})g_{S^{2}}(\theta_{2},\varphi_{2})\}.

5.2.2 A universal kernel and data sampling on $S^{2}\times S^{2}$

Universal kernel. Consider the RKHS ${\mathcal{H}}_{K}$ defined on a Hausdorff topological space ${\cal X}$ . Let ${\cal Z}\subset{\cal X}$ be an arbitrary but fixed compact subset of ${\cal X}$ and denote by ${\mathcal{H}}_{K}({\cal Z})\subset{\mathcal{H}}_{K}$ the completion in the RKHS norm of the span of kernel sections determined by the elements of ${\cal Z}$ . We write this as:

{\mathcal{H}}_{K}({\cal Z})=\overline{{\rm span}\left\{K_{z}\mid z\in{\cal Z}\right\}}.

(5.11)

Denote now by $\overline{\mathcal{H}_{K}({\cal Z})}$ the uniform closure of ${\mathcal{H}}_{K}({\cal Z})$ . A kernel $K$ is called universal if for any compact subset $\mathcal{Z}\subset{\cal X}$ , we have that $\overline{\mathcal{H}_{K}({\cal Z})}=C({\cal Z})$ , with $C({\cal Z})$ the set of real-valued continuous functions on ${\cal Z}$ . Equivalently, this implies that for any $\varepsilon>0$ , and any function $f\in C(\mathcal{Z})$ , there exists a function $g\in\mathcal{H}_{K}({\cal Z})$ , such that $\|f-g\|_{\infty}<\varepsilon$ . Many kernels that are used in practice are indeed universal [Micc 06, Stei 01], e.g., the Gaussian kernel on Euclidean space. For non-Euclidean spaces there is still a lack of literature in this field. The universality of kernels is a crucial feature in recovering unknown functions using the corresponding RKHS. In our learning approach, we shall consider the kernel on $S^{2}\times S^{2}$ defined by

\displaystyle K_{\eta}(x,y):=\exp\left\{-\frac{1}{\eta^{2}}|\iota(x)-\iota(y)|^{2}\right\},\quad\text{ for all}\quad x,y\in S^{2}\times S^{2},

(5.12)

where $\eta>0$ is a constant and $\iota:S^{2}\times S^{2}\to\mathbb{R}^{3}\times\mathbb{R}^{3}$ is the injection map. We emphasize that since $\iota$ is continuous and injective, by [Chri 10, Theorem 2.2], the kernel $K_{\eta}$ in (5.12) is universal for all $\eta>0$ on $S^{2}\times S^{2}$ .

Sampling method. Note that the total surface of the unit sphere is

\displaystyle\int_{0}^{2\pi}\int_{0}^{\pi}\sin\theta\mathrm{d}\theta\mathrm{d}\varphi=4\pi.

The cumulative distribution function can be computed as $\Theta,\Phi$ such that

\displaystyle\mathbb{P}(\Theta\leq\theta,\Phi\leq\varphi)=\frac{1}{4\pi}\int_{0}^{\varphi}\int_{0}^{\theta}\sin t~{}dtds=\frac{1}{4\pi}\varphi(1-\cos\theta),

which factorizes into the product of the margins

\displaystyle\mathbb{P}(\Theta\leq\theta)=\frac{1}{2}(1-\cos\theta),\quad\mathbb{P}(\Phi\leq\varphi)=\frac{1}{2\pi}\varphi.

One can hence obtain a uniform distribution on $S^{2}$ by sampling $U_{1},U_{2}$ uniformly on the interval $(0,1)$ , and compute

\displaystyle\Theta=\arccos(2U_{1}-1),\quad\Phi=2\pi U_{2}.

This extends to a uniform distribution on the product space $\mathcal{P}$ , which we adopt to generate training data.

5.2.3 Exact recovery of Hamiltonians using a universal kernel

In our first example, we consider as Hamiltonian function on $\overline{\mathcal{P}}$ the Euclidean 3-norm on $\mathbb{R}^{3}\times\mathbb{R}^{3}$ .

Formulation We choose the globally defined Hamiltonian function

H^{\prime}\left(x_{1},x_{2}\right)=(\|x_{1}\|^{3}_{3}+\|x_{2}\|^{3}_{3})^{\frac{1}{3}},

where $\|\cdot\|_{3}$ is the usual 3-norm on $\mathbb{R}^{3}$ defined as $\|(a,b,c)\|_{3}=(a^{3}+b^{3}+c^{3})^{\frac{1}{3}}$ . In local coordinates, we focus on the patch $(V_{1}\times V_{1},\phi_{1}\times\phi_{1})$ and so that the Hamiltonian become

\displaystyle H(\theta_{1},\varphi_{1},\theta_{2},\varphi_{2})=H^{\prime}(\phi_{1}(\theta_{1},\varphi_{1}),\phi_{1}(\theta_{2},\varphi_{2})).

We shall call any Hamiltonian function constructed in this way a spherical 3-norm because the domain is on the products of spheres. The restricted Gaussian kernel $K_{\eta}$ as in (5.12) will be adopted for learning, also focused on the patch $(V_{1}\times V_{1},\phi_{1}\times\phi_{1})$ on $\overline{\mathcal{P}}$ .

Numerical results We set $N=1200$ and sample $N$ states from the uniform distribution on $\overline{\mathcal{P}}$ , and then obtain the Hamiltonian vector fields at these states. In view of (4.3), we set the regularization parameter $\lambda=c\cdot N^{-\alpha}$ with $\alpha=0.4$ . We perform hyperparameter tunning and selected parameter $\eta=0.9$ and $c=1e^{-4}$ . Given that $\overline{\mathcal{P}}$ is $4$ -dimensional, we fix two of the dimensions, that is, $(\theta_{2},\varphi_{2})=(\frac{\pi}{2},0)$ , $(\theta_{1},\theta_{2})=(\frac{\pi}{3},\frac{\pi}{3})$ , and visualize both the ground-truth Hamiltonian functions $H$ and the reconstructed Hamiltonian functions $\widehat{h}_{\lambda,N}$ in the remaining $2$ dimensions, see Figure 5.4 (a)(b) and (c)(d). We also visualize the absolute error of the predicted Hamiltonian function in a heatmap, after a vertical shift to offset any potential constant, see Figure 5.4 (e)(f). We also picture the Hamiltonian function globally on the first sphere and the second sphere respectively as heatmaps, while fixing $(\theta_{2},\varphi_{2})=(\frac{\pi}{2},0)$ and $(\theta_{1},\varphi_{1})=(\frac{\pi}{2},0)$ respectively on the other sphere, see Figure 5.5. Lastly, we picture the absolute error of the learned Hamiltonian in the setting of only one sphere and $N$ increased to 2000, see Figure 5.6, which demonstrates the enhanced quantitative performance of the algorithm as the number of data becomes large.

5.2.4 Recovery of two-vortex dynamics

The Hamiltonian function of the two-vortex dynamics is

H^{\prime}\left(x_{1},x_{2}\right)=-\lambda_{1}\lambda_{2}\log\left(1-x_{1}\cdot x_{2}\right),

(5.13)

for $x_{1},x_{2}\in\mathbb{R}^{3}$ , and the Hamiltonian vector field $X_{H^{\prime}}$ is given by

\dot{x}_{j}=X_{H^{\prime}}(\mathbf{x})_{j}=\sum_{i\neq j}\lambda_{i}\frac{x_{i}\times x_{j}}{1-x_{i}\cdot x_{j}}.

Moreover, in local coordinate, the Hamiltonian can be written as

\displaystyle H(\theta_{1},\varphi_{1},\theta_{2},\varphi_{2})=-\lambda_{1}\lambda_{2}\log(1-\sin(\theta_{1})\sin(\theta_{2})\cos(\varphi_{1}-\varphi_{2})-\cos(\theta_{1})\cos(\theta_{2})),

and, as we have seen in (2.13), the Hamiltonian vector field can be locally written as

\displaystyle X_{H}(\theta_{1},\varphi_{1},\theta_{2},\varphi_{2})=J_{\mathcal{P}}(\theta_{1},\varphi_{1},\theta_{2},\varphi_{2})\nabla H(\theta_{1},\varphi_{1},\theta_{2},\varphi_{2}).

Numerical results For the numerical experiment, since the patch $(V_{1}\times V_{1},\phi_{1}\times\phi_{1})$ is topologically dense in the space $\mathcal{P}$ , we shall assume that all the training data will fall into this patch for computational convenience. We consider the standard Gaussian kernel with parameter $\eta$ on $\mathbb{R}^{6}$ but restricted to $\mathcal{P}$ . Since the Gaussian kernel is positive-definite on the entire $\mathbb{R}^{6}$ , it is also positive-definite when restricted to the manifold $\mathcal{P}$ , and hence constitutes a valid choice of a kernel on $\mathcal{P}$ . We set $N=1200$ and sample $N$ states from the uniform distribution on $\mathcal{P}$ as elaborated above, and then obtain the Hamiltonian vector fields at these states. In view of (4.3), we set the regularization parameter $\lambda=c\cdot N^{-\alpha}$ with $\alpha=0.4$ . We performed hyperparameter tunning and set $\eta=0.7$ and $c=0.01$ . Given that $\mathcal{P}$ is $4$ -dimensional, we shall respectively fix two of the dimensions, that is, $(\theta_{2},\varphi_{2})=(\frac{\pi}{2},0)$ , $(\theta_{1},\varphi_{1})=(\frac{\pi}{2},\frac{\pi}{2})$ and $(\theta_{1},\theta_{2})=(\frac{\pi}{3},\frac{\pi}{3})$ , and respectively visualize the mean-square error of the predicted Hamiltonian vector fields in a heatmap, see Figure 5.7 (a)(b)(c). In the end, we picture the Hamiltonian function globally on the first sphere as a heatmap, while fixing $(\theta_{2},\varphi_{2})=(\frac{\pi}{2},0)$ on the second sphere, see Figure 5.8.

Analysis It can be seen that the estimator $\widehat{h}_{\lambda,N}$ fully captures the qualitative behavior of the ground-truth Hamiltonian $H$ . Since $H$ exhibits singularities, by Theorem 3.2 (iii), $H$ cannot be in $\mathcal{H}_{K}$ . In this scenario, even though the restricted Gaussian kernel to the manifold $\mathcal{P}$ is universal, it only guarantees good approximation properties within the category of continuous functions, and not for functions with singularities.

6 Conclusion

This paper proposes a novel kernel-based machine learning method that offers a closed-form solution to the inverse problem of recovering a globally defined and potentially high-dimensional and nonlinear Hamiltonian function on Poisson manifolds. The method uses noisy observations of Hamiltonian vector fields on the phase space as data input. The approach is formulated as a kernel ridge regression problem, where the optimal candidate is sought within a reproducing kernel Hilbert space (RKHS) on the manifold, minimizing the discrepancy between the observed and candidate Hamiltonian vector fields (measured with respect to a chosen Riemannian metric) and an RKHS-based regularization term. Despite the complexities associated with optimization on a manifold, the problem remains convex, leading to an explicit solution derived by leveraging a “differential” version of the reproducing property on the manifold and considering the variations of the loss function.

The kernel method exhibits various key advantages. First, as we demonstrated in Section 3.1, the differential reproducing property allows the kernel ridge regression framework to naturally incorporate the differential equation constraints necessary for structure preservation, even on manifolds. Second, the method provides a closed-form solution, eliminating the computational burdens typically associated with iterative, gradient-based optimization. Additionally, the strict convexity of the Tikhonov-regularized kernel regression ensures a unique solution, addressing the ill-posedness that is inherent in recovering the Hamiltonian function in the presence of Poisson degeneracy due to the potential presence of Casimir functions. Finally, we introduced an operator-theoretic framework and a kernel estimator representation that enables a rigorous analysis of estimation and approximation errors. Additionally, when the target function exhibits a symmetry that is a priori known, as it is common in mechanical systems, an appropriate kernel can be selected to satisfy the symmetry constraints, further enhancing the method’s versatility, applicability, and structure-preservation features.

This paper tackles just one aspect of the broader class of structure-preserving inverse problems, focusing on extending prior approaches from Euclidean spaces to non-Euclidean manifolds. A significant challenge still lies in addressing systems defined in infinite-dimensional spaces, such as the space of measures. A notable example is the Schrödinger equation, which can be interpreted as a transport equation for measures under the Hamiltonian flow. Interestingly, this transport equation can also be viewed as a Hamiltonian flow within the Wasserstein space, referred to as the Wasserstein Hamiltonian flow. Extending the kernel-based framework to identify the potential function in the Schrödinger equation, thereby accommodating infinite-dimensional underlying spaces, is the focus of the authors’ research agenda.

The methodology in this paper assumes that Hamiltonian vector field data is directly available, which is a reasonable assumption when the observed data has a measurable physical interpretation, such as acceleration or force. However, in more realistic scenarios, the collected data might consist of discrete-time trajectories rather than continuous vector field observations. In this context, using variational integrators and other manifold-preserving techniques becomes crucial, yet these methods have not been fully explored in the existing literature in connection with this problem. Addressing this gap could open new avenues for applying structure-preserving kernel methods to a broader range of real-world problems.

As a final remark, a great wealth of knowledge has been accumulated on the qualitative behavior of Hamiltonian dynamical systems based on their geometry [Abra 78], symmetries [Mars 74, Orte 04], stability properties [Orte 05], or persistence [Mont 97b, Mont 97a, Orte 97] and bifurcation phenomena [Chos 02, Chos 03]. All these concepts surely have an interplay in relation to learnability that needs to be explored.

Table 1: Comparison among symplectic vector space, symplectic manifold, and Poisson manifold setups

Symplectic vector space Symplectic manifold Poisson manifold Space $(V,\omega_{\mathrm{can}})$ $(M,\omega)$ $(P,\{\cdot,\cdot\})$ Hamiltonian vector field $X_{h}=J_{\mathrm{can}}\nabla h$ $i_{X_{h}}=\mathbf{d}h$ $X_{h}=B^{\sharp}\mathbf{d}h$ Compatible structure $J_{\mathrm{can}}:TV\to TV$ $J=\omega^{\sharp}g^{\flat}:TM\to TM$ $J=B^{\sharp}g^{\flat}:TP\to TP$ Differential Gram matrix $G_{N}=\mathbb{J}_{\mathrm{can}}\nabla_{1,2}K(\mathbb{Z}_{N},\mathbb{Z}_{N})\mathbb{J}_{\mathrm{can}}^{\top}$ $G_{N}\mathbf{c}:=X_{g_{N}(\mathbf{c},X_{K_{\cdot}}(\mathbf{Z}_{N}))}(\mathbf{Z}_{N})$ $G_{N}\mathbf{c}:=X_{g_{N}(\mathbf{c},X_{K_{\cdot}}(\mathbf{Z}_{N}))}(\mathbf{Z}_{N})$ Estimator $\mathbf{X}_{\sigma^{2},N}^{\top}\left(G_{N}+\lambda NI\right)^{-1}\mathbb{J}_{\mathrm{can}}\nabla_{1}K(\mathbf{Z}_{N},\cdot)$ $g_{N}\left((G_{N}+\lambda NI)^{-1}\mathbf{X}_{\sigma^{2},N},X_{K_{\cdot}}(\mathbf{Z}_{N})\right)$ $g_{N}\left((G_{N}+\lambda NI)^{-1}\mathbf{X}_{\sigma^{2},N},X_{K_{\cdot}}(\mathbf{Z}_{N})\right)$

Appendix A Appendices

A.1 The proof of Theorem 3.3

We prove (i) and (ii) simultaneously by induction on $s$ . The case $s=0$ is trivial since that means for any $x\in M$ , $D^{(0,0)}K(x,\cdot)=K_{x}$ satisfies the standard reproducing property in $\mathcal{H}_{K}$ . Let now $0\leq l\leq s-1$ . Suppose that $D^{(l,0)}K(y,\cdot)\cdot v\in\mathcal{H}_{K}$ and (3.2) holds for any $y\in T^{l-1}M$ and $v\in T_{y}T^{l-1}M$ . Then (3.2) implies that for any $1\leq k,k^{\prime}\leq l-1$ , $x\in T^{k}M$ , $u\in T_{x}T^{k}M$ , $x^{\prime}\in T^{k^{\prime}}M$ , and $u^{\prime}\in T_{x^{\prime}}T^{k^{\prime}}M$ ,

\langle D^{(k,0)}K(x,\cdot)\cdot u,D^{(k^{\prime},0)}K(x^{\prime},\cdot)\cdot u^{\prime}\rangle_{\mathcal{H}_{K}}=D^{(k,k^{\prime})}K(x,x^{\prime})\cdot(u,u^{\prime}).

(A.1)

Denote $K^{k}(y,\cdot)=D^{(k,0)}K(x,\cdot)\cdot u$ and $K^{(k,k^{\prime})}(y,y^{\prime}):=D^{(k,k^{\prime})}K(x,x^{\prime})\cdot(u,u^{\prime})$ with with $y=(x,u)$ and $y^{\prime}=(x^{\prime},u^{\prime})$ .

Now, we turn to the case $l+1\leq s$ . In the following proof, let $\phi:[-1,1]\to T^{l}M$ be a differentiable curve with $\phi(0)=y\in T^{l}M$ and $\phi^{\prime}(0)=v\in T_{y}T^{l}M$ . Since $T^{l}M$ with Riemannian metric $g_{l}$ is complete (see Remark 2.14), by Hopf-Rinow theorem, we can always choose $\phi$ to be geodesics with constant speed 1. We prove that (i) and (ii) hold for $l+1$ in three steps.

Step 1: Proving $D^{(l+1,0)}K(y,\cdot)\cdot v\in\mathcal{H}_{K}$ for all $y\in T^{l}M$ and $v\in T_{y}T^{l}M$ . Equations (2.10) and (A.1) imply that the set $\{\frac{1}{t}(K^{l}(\phi(t))-K^{l}(y)):|t|\leq 1\}$ of functions in $\mathcal{H}_{K}$ satisfies

	$\displaystyle\Big{\\|}\frac{1}{t}($	$\displaystyle K^{l}(\phi(t),\cdot)-K^{l}(y,\cdot))\Big{\\|}^{2}_{\mathcal{H}_{K}}=\frac{1}{t^{2}}\Big{(}K^{(l,l)}(\phi(t),\phi(t))-K^{(l,l)}(y,\phi(t))-K^{(l,l)}(\phi(t),y)+K^{(l,l)}(y,y)\Big{)}$
		$\displaystyle=\frac{1}{t^{2}}\int_{0}^{t}\Big{(}D^{(1,0)}K^{(l,l)}(\phi(s),\phi(t))\cdot\phi^{\prime}(s)-D^{(1,0)}K^{(l,l)}(\phi(s),y)\cdot\phi^{\prime}(s)\Big{)}\mathrm{d}s$
		$\displaystyle=\frac{1}{t^{2}}\int_{0}^{t}\int_{0}^{t}D^{(1,1)}K^{(l,l)}(\phi(s),\phi(\tau))\cdot(\phi^{\prime}(s),\phi^{\prime}(\tau))\mathrm{d}\tau\mathrm{d}s$
		$\displaystyle\leq\frac{1}{t^{2}}\\|K\\|_{C_{b}^{2(l+1)}}\int_{0}^{t}\int_{0}^{t}\\|\phi^{\prime}(\tau)\\|_{l}\\|\phi^{\prime}(s)\\|_{l}\mathrm{d}\tau\mathrm{d}s=\frac{1}{2}\\|K\\|_{C_{b}^{2(l+1)}},\quad\forall~{}\|t\|\leq 1,$

where $\|\cdot\|_{l}$ is induced by the Riemannian metric on $T^{l}M$ and we have used the assumption that $K\in C_{b}^{2s+1}(M\times M)$ . This means that $\{\frac{1}{t}(K^{l}(\phi(t))-K^{l}(y)):|t|\leq 1\}$ lies in a closed ball of the Hilbert space $\mathcal{H}_{K}$ with a finite radius. Since this ball is sequentially weakly compact (see [Conw 07, Theorem 4.2]), there is a sequence $\{t_{i}\}_{i=1}^{\infty}$ with $|t_{i}|\leq 1$ and $\lim_{i\rightarrow\infty}t_{i}=0$ such that $\{\frac{1}{t_{i}}(K^{l}(\phi(t_{i}),\cdot)-K^{l}(y,\cdot)):|t_{i}|\leq 1\}$ converges weakly to an element $h_{y}$ of $\mathcal{H}_{K}$ as $i\rightarrow\infty$ . The weak convergence tells us that

\lim\limits_{i\rightarrow\infty}\Big{\langle}\frac{1}{t_{i}}(K^{l}(\phi(t_{i}),\cdot)-K^{l}(y,\cdot)),f\Big{\rangle}_{\mathcal{H}_{K}}=\langle h_{y},f\rangle_{\mathcal{H}_{K}},\quad\forall f\in\mathcal{H}_{K}.

(A.2)

In particular, by taking $f=K(x,\cdot)$ with $x\in M$ , it holds that

	$\displaystyle h_{y}(x)$	$\displaystyle=\lim\limits_{i\rightarrow\infty}\Big{\langle}\frac{1}{t_{i}}(K^{l}(\phi(t_{i}),\cdot)-K^{l}(y,\cdot)),K(x,\cdot)\Big{\rangle}_{\mathcal{H}_{K}}$
		$\displaystyle=\lim\limits_{i\rightarrow\infty}\frac{1}{t_{i}}(K^{(l,0)}(\phi(t_{i}),x)-K^{(l,0)}K(y,x))=D^{(1,0)}K^{(l,0)}(y,x)\cdot v=D^{(l+1,0)}K(y,x)\cdot v.$

This is true for an arbitrary point $x\in M$ . Hence $D^{(l+1,0)}K(y,\cdot)\cdot v=h_{y}$ as functions on $M$ . Since $h_{y}\in\mathcal{H}_{K}$ , we know $D^{(l+1,0)}K(y,\cdot)\cdot v\in\mathcal{H}_{K}$ .

Step 2: Proving the convergence

\displaystyle\frac{1}{t}\Big{(}K^{l}(\phi(t),\cdot)-K^{l}(y,\cdot)\Big{)}\longrightarrow D^{(l+1,0)}K(y,\cdot)\cdot v,\quad\text{in}\quad\mathcal{H}_{K}.

(A.3)

Using the notation introduced above, $K^{l+1}((y,v),\cdot)=D^{(l+1,0)}K(y,\cdot)\cdot v\in\mathcal{H}_{K}$ . Then equation (A.2) yields

	$\displaystyle\left\langle K^{l+1}((y,v),\cdot),K^{l+1}((y,v),\cdot)\right\rangle_{\mathcal{H}_{K}}$	$\displaystyle=\lim\limits_{i\rightarrow\infty}\frac{1}{t_{i}}\left\langle K^{l+1}((y,v),\cdot),K^{l}(\phi(t_{i}),\cdot)-K^{l}(y,\cdot)\right\rangle_{\mathcal{H}_{K}}$
	$\displaystyle=\$	$\displaystyle\lim\limits_{i\rightarrow\infty}\frac{1}{t_{i}}\left(D^{(1,0)}K^{(l,l)}(y,\phi(t_{i}))\cdot v-D^{(1,0)}K^{(l,l)}(y,y)\cdot v\right)$
	$\displaystyle=\$	$\displaystyle D^{(1,1)}K^{(l,l)}(y,y)\cdot(v,v)=D^{(l+1,l+1)}K(y,y)\cdot(v,v).$

By equation (2.10), we obtain that

		$\displaystyle\Big{\\|}\frac{1}{t}\Big{(}K^{l}(\phi(t),\cdot)-K^{l}(y,\cdot)\Big{)}-K^{l+1}((y,v),\cdot)\Big{\\|}^{2}_{\mathcal{H}_{K}}$
	$\displaystyle=~{}$	$\displaystyle\frac{1}{t^{2}}\Big{(}K^{(l,l)}(\phi(t),\phi(t))-2K^{(l,l)}(\phi(t),y)+K^{(l,l)}(y,y)\Big{)}+K^{(l+1,l+1)}((y,v),(y,v))$
		$\displaystyle-\frac{2}{t}\Big{(}K^{(l+1,l)}((y,v),\phi(t))-K^{(l+1,l)}((y,v),y)\Big{)}$
	$\displaystyle=~{}$	$\displaystyle\frac{1}{t^{2}}\int_{0}^{t}\int_{0}^{t}K^{(l+1,l+1)}((\phi(s),\phi^{\prime}(s)),(\phi(\tau),\phi^{\prime}(\tau)))\mathrm{d}\tau\mathrm{d}s+K^{(l+1,l+1)}((y,v),(y,v))$
		$\displaystyle-\frac{2}{t}\int_{0}^{t}K^{(l+1,l+1)}((\phi(\tau),\phi^{\prime}(\tau)),(u,v))~{}\mathrm{d}\tau$
	$\displaystyle=~{}$	$\displaystyle\frac{1}{t^{2}}\int_{0}^{t}\int_{0}^{t}\left[K^{(l+1,l+1)}((\phi(s),\phi^{\prime}(s)),(\phi(\tau),\phi^{\prime}(\tau))-K^{(l+1,l+1)}((\phi(s),\phi^{\prime}(s)),(y,v))\right]\mathrm{d}\tau\mathrm{d}s$
		$\displaystyle-\frac{1}{t}\int_{0}^{t}\left[K^{(l+1,l+1)}((\phi(\tau),\phi^{\prime}(\tau)),(y,v))-K^{(l+1,l+1)}((y,v),(y,v))\right]\mathrm{d}\tau$
	$\displaystyle=~{}$	$\displaystyle\frac{1}{t^{2}}\int_{0}^{t}\int_{0}^{t}\int_{0}^{\tau}D^{(0,1)}K^{(l+1,l+1)}((\phi(s),\phi^{\prime}(s)),(\phi(z),\phi^{\prime}(z)))\cdot\psi^{\prime}(z)\mathrm{d}z\mathrm{d}\tau\mathrm{d}s$
		$\displaystyle-\frac{1}{t}\int_{0}^{t}\int_{0}^{\tau}D^{(1,0)}K^{(l+1,l+1)}((\phi(z),\phi^{\prime}(z)),(y,v))\cdot\psi^{\prime}(z)\mathrm{d}z\mathrm{d}\tau$
	$\displaystyle=~{}$	$\displaystyle\frac{1}{t^{2}}\int_{0}^{t}\int_{0}^{t}\int_{0}^{\tau}D^{(l+1,l+2)}K(x,x)\cdot((\phi(s),\phi^{\prime}(s)),(\phi(z),\phi^{\prime}(z)))\cdot\psi^{\prime}(z)\mathrm{d}z\mathrm{d}\tau\mathrm{d}s$
		$\displaystyle-\frac{1}{t}\int_{0}^{t}\int_{0}^{\tau}D^{(l+2,l+1)}K(x,x)\cdot((\phi(z),\phi^{\prime}(z)),(y,v))\cdot\psi^{\prime}(z)\mathrm{d}z\mathrm{d}\tau$
	$\displaystyle\leq~{}$	$\displaystyle\frac{1}{t^{2}}\int_{0}^{t}\int_{0}^{t}\int_{0}^{\tau}\\|D^{(l+1,l+2)}K\\|_{\infty}\\|\psi^{\prime}(z)\\|_{l+1}\mathrm{d}z\mathrm{d}\tau\mathrm{d}s+\frac{1}{t}\int_{0}^{t}\int_{0}^{s}\\|D^{(l+2,l+1)}K\\|_{\infty}\\|\psi^{\prime}(z)\\|_{l+1}\mathrm{d}z\mathrm{d}s$
	$\displaystyle=~{}$	$\displaystyle\frac{1}{6}\\|D^{(l+1,l+2)}K\\|_{\infty}~{}t+\frac{1}{2}\\|D^{(l+2,l+1)}K\\|_{\infty}~{}t\leq\frac{1}{2}\\|K\\|_{C_{b}^{2l+3}}~{}t\leq\frac{1}{2}\\|K\\|_{C_{b}^{2s+1}}~{}t,$

where $\psi:[0,t]\to T^{l+1}M$ is chosen as a geodesic with constant speed 1, together with $\psi(0)=(y,v)$ and $\psi(z)=(\phi(z),\phi^{\prime}(z))$ , the second and fourth qualities are due to the formula (2.10), and the last two inequalities thank to $K\in C_{b}^{2s+1}(M\times M)$ and $2l+3\leq 2s+1$ . As $t\rightarrow 0$ , (A.3) follows.

Step 3: Proving (3.2). Let $f\in\mathcal{H}_{K}$ . By (A.3) we have

	$\displaystyle D^{l+1}f(y)\cdot v$	$\displaystyle=\lim\limits_{t\rightarrow 0}\frac{1}{t}\left(D^{l}f(\phi(t))-D^{l}f(y)\right)=\lim\limits_{t\rightarrow 0}\Big{\langle}\frac{1}{t}\left(K^{l}(\phi(t),\cdot)-K^{l}(y,\cdot)\right),f\Big{\rangle}_{\mathcal{H}_{K}}$
		$\displaystyle=\langle K^{(l+1,0)}((y,v),\cdot),f\rangle_{\mathcal{H}_{K}}=\langle D^{(l+1,0)}K(y,\cdot)\cdot v,f\rangle_{\mathcal{H}_{K}}.$

That is, $D^{l+1}f(y)\cdot v$ exists and equals $\langle(D^{(l+1,0)}K(y,\cdot)\cdot v,f\rangle_{\mathcal{H}_{K}}$ . This verifies (3.2) for $l+1$ . Hence, the parts (i) and (ii) follow.
We conclude by proving part (iii) using (3.2) and (A.1). For $f\in\mathcal{H}_{K}$ , and $x\in M$ , the Cauchy-Schwarz inequality implies that,

\displaystyle|f(x)|=|\langle K_{x},f\rangle_{\mathcal{H}_{K}}|\leq\sqrt{K(x,x)}~{}\|f\|_{\mathcal{H}_{K}}\leq\sqrt{\|K\|_{\infty}}~{}\|f\|_{\mathcal{H}_{K}},

and for all $1\leq k\leq s$ , $y\in T^{k-1}M$ and $v\in T_{y}T^{k-1}M$ ,

\displaystyle|D^{k}f(y)\cdot v|

\displaystyle=|\langle(D^{(k,0)}K(y,\cdot)\cdot v,f\rangle_{\mathcal{H}_{K}}|\leq\sqrt{D^{(k,k)}K(y,y)\cdot(v,v)}~{}\|f\|_{\mathcal{H}_{K}}\leq\sqrt{\|D^{(k,k)}K\|_{\infty}}~{}\|v\|_{k-1}\|f\|_{\mathcal{H}_{K}}.

Therefore, $f\in C^{k}_{b}(M)$ and $\|D^{k}f\|_{\infty}\leq\sqrt{\|D^{(k,k)}K\|_{\infty}}~{}\|f\|_{\mathcal{H}_{K}}$ . Denote $\kappa=\sqrt{(s+1)\|K\|_{C_{b}^{2}}}$ . It follows that

	$\displaystyle\\|f\\|_{C^{s}_{b}}$	$\displaystyle=\\|f\\|_{\infty}+\sum_{k=1}^{s}\\|D^{k}f\\|_{\infty}\leq\left(\sqrt{\\|K\\|_{\infty}}+\sum_{k=1}^{s}\sqrt{\\|D^{(k,k)}K\\|_{\infty}}\right)~{}\\|f\\|_{\mathcal{H}_{K}}$
		$\displaystyle\leq\sqrt{(s+1)\left(\\|K\\|_{\infty}+\sum_{k=1}^{s}\\|D^{(k,k)}K\\|_{\infty}\right)}~{}\\|f\\|_{\mathcal{H}_{K}}\leq\sqrt{(s+1)\\|K\\|_{C_{b}^{2}}}~{}\\|f\\|_{\mathcal{H}_{K}}=\kappa\\|f\\|_{\mathcal{H}_{K}}.$

A.2 The proof of Proposition 3.8

Since $K\in C_{b}^{3}(P\times P)$ , Theorem 3.3 yields that for any $h\in\mathcal{H}_{K}$ and $z\in P$ ,

\displaystyle g(z)(\nabla h(z),\nabla h(z))\leq\|h\|_{C_{b}^{1}}\left(g(z)(\nabla h(z),\nabla h(z))\right)^{\frac{1}{2}}\leq\kappa\|h\|_{\mathcal{H}_{K}}\left(g(z)(\nabla h(z),\nabla h(z))\right)^{\frac{1}{2}},

which implies that $g(z)(\nabla h(z),\nabla h(z))\leq\kappa^{2}\|h\|_{\mathcal{H}_{K}}^{2}$ . Then by Assumption 3.7 and equation (2.13), we have that

\displaystyle g(z)(X_{h}(z),X_{h}(z))=g(z)(J(z)\nabla h(z),J(z)\nabla h(z))\leq\gamma(z)g(z)(\nabla h(z),\nabla h(z))\leq\kappa^{2}\gamma(z)\|h\|_{\mathcal{H}_{K}}^{2}.

(A.4)

Thus, the operator $A$ is bounded linear from $\mathcal{H}_{K}$ to $L^{2}(P,\mu_{\mathbf{Z}})$ since

\displaystyle\|Ah\|^{2}_{L^{2}(\mu_{\mathbf{Z}})}=\int_{P}g(z)(X_{h}(z),X_{h}(z))\mathrm{d}\mu_{\mathbf{Z}}(z)\leq\int_{P}\kappa^{2}\gamma(z)\|h\|_{\mathcal{H}_{K}}^{2}\mu_{\mathbf{Z}}(z)\leq\kappa^{2}C^{2}\|h\|_{\mathcal{H}_{K}}^{2},

Furthermore, for any vector field $Y\in L^{2}(P,\mu_{\mathbf{Z}})$ , Corollary 3.5 and Proposition 2.18 imply that

	$\displaystyle\langle Ah,Y\rangle_{L^{2}(\mu_{\mathbf{Z}})}$	$\displaystyle=\int_{P}g(z)(X_{h}(z),Y(z))\mathrm{d}\mu_{\mathbf{Z}}(z)=\int_{P}g(z)\left(\nabla h(z),-J(z)Y(z)\right)\mathrm{d}\mu_{\mathbf{Z}}(z)$
		$\displaystyle=\int_{P}(-JY)[h](z)\mathrm{d}\mu_{\mathbf{Z}}(z)=\int_{P}(-JY)[\langle h,K_{\cdot}\rangle_{\mathcal{H}_{K}}](z)\mathrm{d}\mu_{\mathbf{Z}}(z)$
		$\displaystyle=\int_{P}\frac{d}{dt}\mid_{t=0}[\langle h,K_{F_{t}(z)}\rangle_{\mathcal{H}_{K}}]\mathrm{d}\mu_{\mathbf{Z}}(z)=\left\langle h,\int_{P}(-JY)[K_{\cdot}](z)\mathrm{d}\mu_{\mathbf{Z}}(z)\right\rangle_{\mathcal{H}_{K}}.$

where $F_{t}$ is the flow of the vector field $-JY$ and Remark 3.4(i) shows that $(-JY)[K_{\cdot}](z)\in\mathcal{H}_{K}$ . This leads to

\displaystyle(A^{*}Y)(y)=\int_{P}(-JY)[K_{y}](z)\mathrm{d}\mu_{\mathbf{Z}}(z)=\int_{P}g(z)(X_{K_{y}}(z),Y(z))\mathrm{d}\mu_{\mathbf{Z}}(z).

Note that hence, for each $h\in\mathcal{H}_{K}$ ,

\displaystyle(Qh)(y):

\displaystyle=(A^{*}Ah)(y)=(A^{*}X_{h})(y)=\int_{P}g(z)(X_{K_{y}}(z),X_{h}(z))\mathrm{d}\mu_{\mathbf{Z}}(z).

We now show that $Q$ is a trace class operator; that is, we show that $\operatorname{Tr}(|Q|)<\infty$ , where $|Q|=\sqrt{Q^{*}Q}$ . The form of the operator $Q=A^{*}A$ automatically guarantees that it is positive semidefinite, and hence we have that $|Q|=Q$ . Therefore, it is equivalent to show that $\operatorname{Tr}(Q)<\infty$ . To do that, we choose a countable spanning orthonormal set $\left\{e_{n}\right\}_{n\in\mathbb{N}}$ for ${\mathcal{H}}_{K}$ whose existence is guaranteed by the continuity of the canonical feature map associated to $K$ that we established in [Hu 24, Lemma A.3] and [Owha 17, Theorem 2.4]. Then,

	$\displaystyle\operatorname{Tr}(Q)$	$\displaystyle=\operatorname{Tr}\left(A^{}A\right)=\sum_{n}\left\langle A^{}Ae_{n},e_{n}\right\rangle_{\mathcal{H}_{K}}=\sum_{n}\left\langle Ae_{n},Ae_{n}\right\rangle_{L^{2}\left(\mu_{\mathbf{Z}}\right)}$
		$\displaystyle=\sum_{n}\int_{P}g(z)(X_{e_{n}}(z),X_{e_{n}}(z))\mathrm{d}\mu_{\mathbf{Z}}(\mathbf{x})$
		$\displaystyle\leq C^{2}\sum_{n}\int_{P}g(z)(\nabla{e_{n}}(z),\nabla{e_{n}}(z))\mathrm{d}\mu_{\mathbf{Z}}(\mathbf{x})$
		$\displaystyle\leq dC^{2}\\|K\\|_{C_{b}^{2}}=\frac{1}{2}dC^{2}\kappa^{2}<\infty,$

where the sixth inequality is derived by choosing an arbitrary orthonormal basis $\{v_{i}\}_{i=1}^{d}$ of the vector space $T_{z}P$ as follows.

	$\displaystyle\sum_{n}g(z)(\nabla{e_{n}}(z),\nabla{e_{n}}(z))$	$\displaystyle=\sum_{n}De_{n}(z)\cdot\nabla{e_{n}}(z)=\sum_{n}\langle D^{(1,0)}K(z,\cdot)\cdot\nabla{e_{n}}(z),e_{n}\rangle_{\mathcal{H}_{K}}$
		$\displaystyle=\sum_{n}\left\langle D^{(1,0)}K(z,\cdot)\cdot\left(\sum_{i=1}^{d}g(z)(\nabla{e_{n}}(z),v_{i})v_{i}\right),e_{n}\right\rangle_{\mathcal{H}_{K}}$
		$\displaystyle=\sum_{n}\sum_{i=1}^{d}\left\langle D^{(1,0)}K(z,\cdot)\cdot v_{i},g(z)(\nabla{e_{n}}(z),v_{i})e_{n}\right\rangle_{\mathcal{H}_{K}}$
		$\displaystyle=\sum_{i=1}^{d}\left\langle D^{(1,0)}K(z,\cdot)\cdot v_{i},\sum_{n}(D{e_{n}}(z)\cdot v_{i})e_{n}\right\rangle_{\mathcal{H}_{K}}$
		$\displaystyle=\sum_{i=1}^{d}\left\langle D^{(1,0)}K(z,\cdot)\cdot v_{i},\sum_{n}\langle D^{(1,0)}K(z,\cdot)\cdot v_{i},e_{n}\rangle_{\mathcal{H}_{K}}e_{n}\right\rangle_{\mathcal{H}_{K}}$
		$\displaystyle=\sum_{i=1}^{d}\left\langle D^{(1,0)}K(z,\cdot)\cdot v_{i},D^{(1,0)}K(z,\cdot)\cdot v_{i}\right\rangle_{\mathcal{H}_{K}}$
		$\displaystyle=\sum_{i=1}^{d}D^{(1,1)}K(z,z)\cdot(v_{i},v_{i})\leq d\\|K\\|_{C_{b}^{2}}.$

A.3 The proof of Proposition 3.12

Firstly, we prove the symmetry of $G_{N}$ . It is sufficient to show that $g_{N}(\mathbf{c},G_{N}\mathbf{c}^{\prime})=g_{N}(\mathbf{c}^{\prime},G_{N}\mathbf{c})$ for $\mathbf{c},\mathbf{c}^{\prime}\in T_{\mathbf{Z}_{N}}P$ arbitrary. Denote $\mathbb{J}_{N}=\operatorname{diag}\{J(\mathbf{Z}^{(1)}),\ldots,J(\mathbf{Z}^{(N)})\}$ . Then by Remark 3.4(i), we obtain that for any $\mathbf{c}\in T_{\mathbf{Z}_{N}}P$ ,

\displaystyle g_{N}(\mathbf{c},X_{K_{\cdot}}(\mathbf{Z}_{N}))=D^{(1,0)}K(\mathbf{Z}_{N},\cdot)\cdot\left(-\mathbb{J}_{N}\mathbf{c}\right).

Therefore, again by Remark 3.4(i), we have that

\displaystyle g_{N}(\mathbf{c},G_{N}\mathbf{c}^{\prime})=g_{N}(\mathbf{c},X_{D^{(1,0)}K(\mathbf{Z}_{N},\cdot)\cdot\left(-\mathbb{J}_{N}\mathbf{c}^{\prime}\right)}(\mathbf{Z}_{N}))=D^{(1,1)}K(\mathbf{Z}_{N},\mathbf{Z}_{N})(-\mathbb{J}_{N}\mathbf{c}^{\prime},-\mathbb{J}_{N}\mathbf{c})

Since $K$ is a symmetric function, $D^{(1,1)}K(\mathbf{Z}_{N},\mathbf{Z}_{N})$ is also symmetric by definition. This proves that $G_{N}$ is symmetric.

We now show the positive semi-definiteness of $G_{N}$ . Let $L=dN$ . Since $G_{N}$ is real symmetric, there exists an orthonormal matrix $O\in T_{\mathbf{Z}_{N}}P\times T_{\mathbf{Z}_{N}}P$ that diagonalizes $G_{N}$ , that is

\displaystyle G_{N}=ODO^{\top}=\begin{bmatrix}|&|&\dots&|\\ f_{1}&f_{2}&\dots&f_{L}\\ |&|&\dots&|\end{bmatrix}\begin{bmatrix}d_{1}&0&\dots&0\\ 0&d_{2}&\dots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&d_{L}\end{bmatrix}\begin{bmatrix}|&|&\dots&|\\ f_{1}&f_{2}&\dots&f_{L}\\ |&|&\dots&|\end{bmatrix}^{\top},

where $\left\{d_{i}\right\}_{i=1}^{L}$ and $\left\{f_{i}\right\}_{i=1}^{L}$ are the real eigenvalues and the corresponding eigenvectors of $G_{N}$ , respectively. We now define $\widetilde{e}_{i}=g_{N}\left(f_{i},X_{K_{\cdot}}(\mathbf{Z}_{N})\right)$ . By Remark 3.4(i), we obtain that $\widetilde{e}_{i}\in\mathcal{H}_{K}$ and that

	$\displaystyle\\|\widetilde{e}_{i}\\|_{\mathcal{H}_{K}}^{2}$	$\displaystyle=\left\langle g_{N}\left(f_{i},X_{K_{\cdot}}(\mathbf{Z}_{N})\right),g_{N}\left(f_{i},X_{K_{\cdot}}(\mathbf{Z}_{N})\right)\right\rangle_{\mathcal{H}_{K}}$
		$\displaystyle=\left\langle D^{(1,0)}K(\mathbf{Z}_{N},\cdot)\cdot(-\mathbb{J}_{N}f_{i}),D^{(1,0)}K(\mathbf{Z}_{N},\cdot)\cdot(-\mathbb{J}_{N}f_{i})\right\rangle_{\mathcal{H}_{K}}$
		$\displaystyle=D^{(1,1)}K(\mathbf{Z}_{N},\mathbf{Z}_{N})(-\mathbb{J}_{N}f_{i},-\mathbb{J}_{N}f_{i})=g_{N}(f_{i},G_{N}f_{i})=d_{i}.$

Hence, $d_{i}\geq 0$ for all $i=1,\ldots,L$ and we can conclude that $G_{N}$ is positive semi-definite.

A.4 The proof of Lemma 4.1

Proof.

Since $K\in C_{b}^{3}(P\times P)$ , then for any $h\in\mathcal{H}_{K}$ , we have

	$\displaystyle\\|Q_{N}h\\|_{\mathcal{H}_{K}}^{2}$	$\displaystyle=\frac{1}{N}\\|g_{N}(X_{h}(\mathbf{Z}_{N}),X_{K_{\cdot}}(\mathbf{Z}_{N}))\\|^{2}_{\mathcal{H}_{K}}=\frac{1}{N}\\|D^{(1,0)}K(\mathbf{Z}_{N},\cdot)\cdot(-JX_{h})(\mathbf{Z}_{N})\\|^{2}_{\mathcal{H}_{K}}$
		$\displaystyle=\frac{1}{N}D^{(1,1)}K(\mathbf{Z}_{N},\mathbf{Z}_{N})\cdot((-JX_{h})(\mathbf{Z}_{N}),(-JX_{h})(\mathbf{Z}_{N}))$
		$\displaystyle\leq\frac{1}{N}\\|D^{(1,1)}K\\|_{\infty}g_{N}(-JX_{h}(\mathbf{Z}_{N}),-JX_{h}(\mathbf{Z}_{N}))\leq\frac{C^{2}}{2N}\kappa^{2}g_{N}(X_{h}(\mathbf{Z}_{N}),X_{h}(\mathbf{Z}_{N}))$
		$\displaystyle\leq\frac{1}{2}C^{4}\kappa^{4}\\|h\\|^{2}_{\mathcal{H}_{K}}<\infty,$

which shows that $Q_{N}h$ are bounded random variables in $\mathcal{H}_{K}$ for all $N\geq 1$ . Moreover, we have $\mathbb{E}\left[\|Q_{N}h\|^{2}_{\mathcal{H}_{K}}\right]\leq\frac{1}{2}C^{4}\kappa^{4}\|h\|^{2}_{\mathcal{H}_{K}}$ . Define now the $\mathcal{H}_{K}$ -valued random variables

\xi^{(n)}=g(\mathbf{Z}^{(n)})(X_{h}(\mathbf{Z}^{(n)}),X_{K_{\cdot}}(\mathbf{Z}^{(n)})),\quad\mbox{$n=1,\ldots,N$.}

Note that the random variables $\{\xi^{(n)}\}_{n=1}^{N}$ are IID and that $Q_{N}h-Qh=\frac{1}{N}\sum_{n=1}^{N}(\xi^{(n)}-\mathbb{E}(\xi^{(n)})).$ The result follows by applying [De V 05, Lemma 8] or [Yuri 95, Lemma A.2] to $\{\xi^{(n)}\}_{n=1}^{N}$ . ∎

A.5 The proof of Lemma 4.2

Proof.

We first notice that, by the symmetry of $G_{N}$ , for any $\mathbf{c},\mathbf{c}^{\prime}\in T_{\mathbf{Z}_{N}}P$ , it holds that

g_{N}((G_{N}+\lambda NI)^{-1}\mathbf{c},\mathbf{c}^{\prime})=g_{N}(\mathbf{c},(G_{N}+\lambda NI)^{-1}\mathbf{c}^{\prime}).

(A.5)

Following an approach similar to the one in Theorem 3.13, we obtain that

	$\displaystyle\big{\\|}\widehat{h}_{\lambda,N}-\widetilde{h}_{\lambda,N}\big{\\|}_{\mathcal{H}_{K}}^{2}$	$\displaystyle=\\|g_{N}(X_{K_{\cdot}}(\mathbf{Z}_{N}),(G_{N}+\lambda NI)^{-1}\mathbf{E}_{N})\\|_{\mathcal{H}_{K}}^{2}$
		$\displaystyle=D^{(1,1)}K(\mathbf{Z}_{N},\mathbf{Z}_{N})(-\mathbb{J}_{N}(G_{N}+\lambda NI)^{-1}\mathbf{E}_{N},-\mathbb{J}_{N}(G_{N}+\lambda NI)^{-1}\mathbf{E}_{N})$
		$\displaystyle:=g_{N}(\mathbf{E}_{N},\Sigma_{N}\mathbf{E}_{N}),$

where $G_{N}$ is defined in (3.15) and $\Sigma_{N}$ is a matrix in $T_{\mathbf{Z}_{N}}P\times T_{\mathbf{Z}_{N}}P$ defined by

\displaystyle\Sigma_{N}\mathbf{c}:=(G_{N}+\lambda NI)^{-1}X_{g_{N}((G_{N}+\lambda NI)^{-1}\mathbf{c},X_{K_{\cdot}}(\mathbf{Z}_{N}))}(\mathbf{Z}_{N}),\quad\mathbf{c}\in T_{\mathbf{Z}_{N}}P.

Indeed, in light of (A.5), the above definition can be verified since

	$\displaystyle g_{N}(\mathbf{E}_{N},\Sigma_{N}\mathbf{E}_{N})$	$\displaystyle=g_{N}(\mathbf{E}_{N},(G_{N}+\lambda NI)^{-1}X_{g_{N}((G_{N}+\lambda NI)^{-1}\mathbf{E}_{N},X_{K_{\cdot}}(\mathbf{Z}_{N}))}(\mathbf{Z}_{N}))$
		$\displaystyle=g_{N}((G_{N}+\lambda NI)^{-1}\mathbf{E}_{N},X_{g_{N}((G_{N}+\lambda NI)^{-1}\mathbf{E}_{N},X_{K_{\cdot}}(\mathbf{Z}_{N}))}(\mathbf{Z}_{N}))$
		$\displaystyle=D^{(1,1)}K(\mathbf{Z}_{N},\mathbf{Z}_{N})(-\mathbb{J}_{N}(G_{N}+\lambda NI)^{-1}\mathbf{E}_{N},-\mathbb{J}_{N}(G_{N}+\lambda NI)^{-1}\mathbf{E}_{N}).$

Choose a basis of the vector space $T_{\mathbf{Z}_{N}}P$ consisting of orthonormal eigenvectors $\{f_{i}\}_{i=1}^{dN}$ of $G_{N}$ . Note that $\{f_{i}\}_{i=1}^{dN}$ are also eigenvectors of $\Sigma_{N}$ . Then,

	$\displaystyle\mathrm{Tr}(\Sigma_{N})$	$\displaystyle=\sum_{i=1}^{dN}g_{N}(f_{i},\Sigma_{N}f_{i})\leq\frac{1}{\lambda^{2}N^{2}}\sum_{i=1}^{dN}D^{(1,1)}K(\mathbf{Z}_{N},\mathbf{Z}_{N})(-\mathbb{J}_{N}f_{i},-\mathbb{J}_{N}f_{i})$
		$\displaystyle\leq\frac{1}{\lambda^{2}N^{2}}\sum_{i=1}^{dN}\\|D^{(1,1)}K\\|_{\infty}\cdot g_{N}(-\mathbb{J}_{N}f_{i},-\mathbb{J}_{N}f_{i})$
		$\displaystyle\leq\frac{1}{\lambda^{2}N^{2}}\\|K\\|_{C_{b}^{2}}\sum_{i=1}^{dN}g_{N}(-\mathbb{J}_{N}f_{i},-\mathbb{J}_{N}f_{i})$
		$\displaystyle\leq\frac{1}{\lambda^{2}N^{2}}\\|K\\|_{C_{b}^{2}}\sum_{i=1}^{dN}\left(\sum_{j=1}^{N}g(-Jf^{(j)}_{i},-Jf^{(j)}_{i})\right)$
		$\displaystyle\leq\frac{1}{\lambda^{2}N^{2}}\\|K\\|_{C_{b}^{2}}\sum_{i=1}^{dN}\left(\sum_{j=1}^{N}C^{2}\cdot g(f^{(j)}_{i},f^{(j)}_{i})\right)$
		$\displaystyle=\frac{C^{2}}{\lambda^{2}N^{2}}\\|K\\|_{C_{b}^{2}}\sum_{i=1}^{dN}g_{N}(f_{i},f_{i})=\frac{dC^{2}}{\lambda^{2}N}\\|K\\|_{C_{b}^{2}}=\frac{\kappa^{2}C^{2}d}{2\lambda^{2}N},$

with $\kappa^{2}=2\|K\|_{C_{b}^{2}}$ as introduced in the statement of Proposition 3.3. The rest of the proof follows from the same argument as in [Hu 24, Appendix C]. ∎

Acknowledgments

The authors thank Lyudmila Grigoryeva for helpful discussions and remarks and acknowledge partial financial support from the School of Physical and Mathematical Sciences of the Nanyang Technological University. DY is funded by a Nanyang President’s Graduate Scholarship of Nanyang Technological University.

References

[Abra 78] R. Abraham and J. E. Marsden. Foundations of Mechanics. Addison-Wesley, Reading, MA, 2nd Ed., 1978.
[Aron 50] N. Aronszajn. “Theory of reproducing kernels”. Transactions of the American mathematical society, Vol. 68, No. 3, pp. 337–404, 1950.
[Baja 23] J. Bajārs. “Locally-symplectic neural networks for learning volume-preserving dynamics”. Journal of Computational Physics, Vol. 476, p. 111911, 2023.
[Beck 22] T. Beckers, J. Seidman, P. Perdikaris, and G. J. Pappas. “Gaussian Process Port-Hamiltonian Systems: Bayesian Learning with Physics Prior”. In: 2022 IEEE 61st Conference on Decision and Control (CDC), pp. 1447–1453, 2022.
[Boro 20] V. Borovitskiy, A. Terenin, P. Mostowsky, and M. P. Deisenroth. “Matérn Gaussian processes on riemannian manifolds”. In: Proceedings of the 34th International Conference on Neural Information Processing Systems, Curran Associates Inc., Red Hook, NY, USA, 2020.
[Bran 23] B. Brantner and M. Kraus. “Symplectic Autoencoders for Model Reduction of Hamiltonian Systems”. 2023.
[Brun 16] S. L. Brunton, J. L. Proctor, and J. N. Kutz. “Discovering governing equations from data by sparse identification of nonlinear dynamical systems”. Proceedings of the National Academy of Sciences, Vol. 113, No. 15, pp. 3932–3937, 2016.
[Chen 19] Z. Chen, J. Zhang, M. Arjovsky, and L. Bottou. “Symplectic recurrent neural networks”. arXiv preprint arXiv:1909.13334, 2019.
[Chen 21] R. Chen and M. Tao. “Data-driven Prediction of General Hamiltonian Dynamics via Learning Exactly-Symplectic Maps”. In: M. Meila and T. Zhang, Eds., Proceedings of the 38th International Conference on Machine Learning, pp. 1717–1727, PMLR, 18–24 Jul 2021.
[Chen 22] X. Chen, J. Duan, J. Hu, and D. Li. “Data-driven method to learn the most probable transition pathway and stochastic differential equation”. Physica D: Nonlinear Phenomena, 2022.
[Chos 02] P. Chossat, J.-P. Ortega, and T. S. Ratiu. “Hamiltonian Hopf bifurcation with symmetry”. Archive for Rational Mechanics and Analysis, Vol. 163, pp. 1–33, 2002.
[Chos 03] P. Chossat, D. Lewis, J.-P. Ortega, and T. S. Ratiu. “Bifurcation of relative equilibria in mechanical systems with symmetry”. Advances in Applied Mathematics, Vol. 31, pp. 10–45, 2003.
[Chri 10] A. Christmann and I. Steinwart. “Universal Kernels on Non-Standard Input Spaces”. In: J. Lafferty, C. Williams, J. Shawe-Taylor, R. Zemel, and A. Culotta, Eds., Advances in Neural Information Processing Systems, p. , Curran Associates, Inc., 2010.
[Conw 07] J. B. Conway. A Course in Functional Analysis. Springer Science & Business Media, second Ed., 2007.
[Cran 20] M. Cranmer, S. Greydanus, S. Hoyer, P. Battaglia, D. Spergel, and S. Ho. “Lagrangian neural networks”. arXiv preprint arXiv:2003.04630, 2020.
[Da C 23a] N. Da Costa, C. Mostajeran, and J.-P. Ortega. “The Gaussian kernel on the circle and spaces that admit isometric embeddings of the circle”. In: International Conference on Geometric Science of Information, pp. 426–435, Springer, 2023.
[Da C 23b] N. Da Costa, C. Mostajeran, J.-P. Ortega, and S. Said. “Geometric Learning with Positively Decomposable Kernels”. arXiv preprint arXiv:2310.13821, 2023.
[Davi 23] M. David and F. Méhats. “Symplectic learning for Hamiltonian neural networks”. Journal of Computational Physics, Vol. 494, p. 112495, 2023.
[De V 05] E. De Vito, L. Rosasco, A. Caponnetto, U. De Giovannini, F. Odone, and P. Bartlett. “Learning from Examples as an Inverse Problem.”. Journal of Machine Learning Research, Vol. 6, No. 5, 2005.
[Doum 23] N. Doumèche, G. Biau, and C. Boyer. “Convergence and error analysis of PINNs”. 2023.
[Eldr 24a] C. Eldred, F. Gay-Balmaz, and V. Putkaradze. “CLPNets: Coupled Lie-Poisson Neural Networks for Multi-Part Hamiltonian Systems with Symmetries”. arXiv preprint arXiv:2408.16160, 2024.
[Eldr 24b] C. Eldred, F. Gay-Balmaz, S. Huraka, and V. Putkaradze. “Lie–Poisson Neural Networks (LPNets): Data-based computing of Hamiltonian systems with symmetries”. Neural Networks, Vol. 173, p. 106162, 2024.
[Emer 12] M. Émery. Stochastic calculus in manifolds. Springer Science & Business Media, 2012.
[Feng 23] J. Feng, C. Kulick, Y. Ren, and S. Tang. “Learning particle swarming models from data with Gaussian processes”. Mathematics of Computation, 2023.
[Fera 15] A. Feragen, F. Lauze, and S. Hauberg. “Geodesic exponential kernels: When curvature and linearity conflict”. In: Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 3032–3042, 2015.
[Ferr 12] J. C. Ferreira and V. A. Menegatto. “Reproducing properties of differentiable Mercer-like kernels”. Mathematische Nachrichten, Vol. 285, No. 8-9, pp. 959–973, 2012.
[Gins 12] D. Ginsbourger, X. Bay, O. Roustant, and L. Carraro. “Argumentwise invariant kernels for the approximation of invariant functions”. In: Annales de la Faculté des sciences de Toulouse: Mathématiques, pp. 501–527, 2012.
[Grey 19] S. Greydanus, M. Dzamba, and J. Yosinski. “Hamiltonian neural networks”. Advances in neural information processing systems, Vol. 32, 2019.
[Holm 09] D. D. Holm, T. Schmah, and C. Stoica. Geometric mechanics and symmetry: from finite to infinite dimensions. Vol. 12, Oxford University Press, 2009.
[Holm 11a] D. D. Holm. Geometric Mechanics. Part I: Dynamics and Symmetry. World Scientific Publishing Company, second Ed., 2011.
[Holm 11b] D. D. Holm. Geometric Mechanics. Part II: Rotating, Translating and Rolling. World Scientific, second Ed., 2011.
[Hsu 02] E. P. Hsu. Stochastic analysis on manifolds. American Mathematical Soc., 2002.
[Hu 24] J. Hu, J.-P. Ortega, and D. Yin. “A Structure-Preserving Kernel Method for Learning Hamiltonian Systems”. arXiv preprint arXiv:2403.10070, 2024.
[Jin 20] P. Jin, Z. Zhang, A. Zhu, Y. Tang, and G. E. Karniadakis. “SympNets: Intrinsic structure-preserving symplectic networks for identifying Hamiltonian systems”. Neural Networks, Vol. 132, pp. 166–179, 2020.
[Jin 22] P. Jin, Z. Zhang, I. G. Kevrekidis, and G. E. Karniadakis. “Learning Poisson systems and trajectories of autonomous systems via Poisson neural networks”. IEEE Transactions on Neural Networks and Learning Systems, 2022.
[Kost 09] B. Kostant. “Orbits, symplectic structures and representation theory”. Collected Papers, pp. 482–482, 2009.
[Leok 12] M. Leok and T. Shingel. “General techniques for constructing variational integrators”. Frontiers of Mathematics in China, Vol. 7, pp. 273–303, 2012.
[Leon 97] N. E. Leonard and J. E. Marsden. “Stability and drift of underwater vehicle dynamics: mechanical systems with rigid motion symmetry”. Physica D: Nonlinear Phenomena, Vol. 105, No. 1-3, pp. 130–162, 1997.
[Li 23] S. Li. “Gaussian kernels on nonsimply connected closed Riemannian manifolds are never positive definite”. Bulletin of the London Mathematical Society, 2023.
[Libe 12] P. Libermann and C.-M. Marle. Symplectic Geometry and Analytical Mechanics. Vol. 35, Springer Science & Business Media, 2012.
[Lie 93] S. Lie and F. Engel. Theorie der Transformationsgruppen. Vol. 3, Teubner, 1893.
[Lu 19] F. Lu, M. Zhong, S. Tang, and M. Maggioni. “Nonparametric inference of interaction laws in systems of agents from trajectory data”. Proceedings of the National Academy of Sciences, Vol. 116, No. 29, pp. 14424–14433, 2019.
[Magg 21] M. Maggioni, J. Miller, H. Qiu, and M. Zhong. “Learning Interaction Kernels for Agent Systems on Riemannian Manifolds”. 2021.
[Mars 01] J. E. Marsden and M. West. “Discrete mechanics and variational integrators”. Acta numerica, Vol. 10, pp. 357–514, 2001.
[Mars 74] J. E. Marsden and A. Weinstein. “Reduction of symplectic manifolds with symmetry”. Reports on Mathematical Physics, Vol. 5, No. 1, pp. 121–130, 1974.
[Mars 99] J. E. Marsden and T. S. Ratiu. Introduction to Mechanics and Symmetry. Springer-Verlag, New York, second Ed., 1999.
[Micc 06] C. A. Micchelli, Y. Xu, and H. Zhang. “Universal Kernels.”. Journal of Machine Learning Research, Vol. 7, No. 12, 2006.
[Minh 10] H. Q. Minh. “Some properties of Gaussian reproducing kernel Hilbert spaces and their implications for function approximation and learning theory”. Constructive Approximation, Vol. 32, pp. 307–338, 2010.
[Mont 97a] J. A. Montaldi. “Persistance d’orbites périodiques relatives dans les systèmes hamiltoniens symétriques”. C. R. Acad. Sci. Paris Sér. I Math., Vol. 324, pp. 553–558, 1997.
[Mont 97b] J. A. Montaldi. “Persistence and stability of relative equilibria”. Nonlinearity, Vol. 10, pp. 449–466, 1997.
[Mrou 15] Y. Mroueh, S. Voinea, and T. A. Poggio. “Learning with Group Invariant Features: A Kernel Perspective.”. Advances in neural information processing systems, Vol. 28, 2015.
[Noet 18] E. Noether. “Invariante Variationsprobleme”. Nachr. d. König. Gesellsch. d. Wiss. zu Göttingen Math. Physik, Vol. 2, pp. 235–257, 1918.
[Nomi 61] K. Nomizu and H. Ozeki. “The existence of complete Riemannian metrics”. Proceedings of the American Mathematical Society, Vol. 12, No. 6, pp. 889–891, 1961.
[Nova 18] E. Novak, M. Ullrich, H. Woźniakowski, and S. Zhang. “Reproducing kernels of Sobolev spaces on $\mathbb{R}^{d}$ and applications to embedding constants and tractability”. Analysis and Applications, Vol. 16, No. 05, pp. 693–715, 2018.
[Offe 24] C. Offen. “Machine learning of continuous and discrete variational ODEs with convergence guarantee and uncertainty quantification”. 2024.
[Orte 04] J.-P. Ortega and T. S. Ratiu. Momentum Maps and Hamiltonian Reduction. Birkhauser Verlag, 2004.
[Orte 05] J.-P. Ortega, V. Planas-Bielsa, and T. S. Ratiu. “Asymptotic and Lyapunov stability of constrained and Poisson equilibria”. Journal of Differential Equations, Vol. 214, No. 1, pp. 92–127, jul 2005.
[Orte 24] J.-P. Ortega and D. Yin. “Learnability of Linear Port-Hamiltonian Systems”. Journal of Machine Learning Research, Vol. 25, No. 68, pp. 1–56, 2024.
[Orte 97] J.-P. Ortega and T. S. Ratiu. “Persistence and smoothness of critical relative elements in Hamiltonian systems with symmetry”. Comptes Rendus de l’Académie des Sciences - Series I - Mathematics, Vol. 325, No. 10, pp. 1107–1111, nov 1997.
[Owha 17] H. Owhadi and C. Scovel. “Separability of reproducing kernel spaces”. Proceedings of the American Mathematical Society, Vol. 145, No. 5, pp. 2131–2138, 2017.
[Penn 06] X. Pennec. “Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements”. Journal of Mathematical Imaging and Vision, Vol. 25, pp. 127–154, 2006.
[Penn 19] X. Pennec, S. Sommer, and T. Fletcher. Riemannian Geometric Statistics in Medical Image Analysis. Academic Press, 2019.
[Pfö 24] M. Pförtner, I. Steinwart, P. Hennig, and J. Wenger. “Physics-Informed Gaussian Process Regression Generalizes Linear PDE Solvers”. 2024.
[Rais 19] M. Raissi, P. Perdikaris, and G. Karniadakis. “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations”. Journal of Computational Physics, Vol. 378, pp. 686–707, 2019.
[Rath 21] K. Rath, C. G. Albert, B. Bischl, and U. von Toussaint. “Symplectic Gaussian process regression of maps in Hamiltonian systems”. Chaos: An Interdisciplinary Journal of Nonlinear Science, Vol. 31, No. 5, p. 053121, 05 2021.
[Rude 13] M. Rudelson, R. Vershynin, et al. “Hanson-wright inequality and sub-Gaussian concentration”. Electronic Communications in Probability, Vol. 18, 2013.
[Shin 20] J. Shin, YeonjongDarbon and G. Em Karniadakis. “On the Convergence of Physics Informed Neural Networks for Linear Second-Order Elliptic and Parabolic Type PDEs”. Communications in Computational Physics, Vol. 28, No. 5, pp. 2042–2074, 2020.
[Sour 65] J.-M. Souriau. Géométrie de l’espace de phases: calcul des variations et mécanique quantique. Faculté des sciences de Marseille, 1965.
[Sour 66] J.-M. Souriau. “Quantification géométrique”. Comm. Math. Phys., Vol. 1, pp. 374–398, 1966.
[Sour 69] J.-M. Souriau. Structure des Systèmes Dynamiques. Dunod, Paris, 1969.
[Stei 01] I. Steinwart. “On the influence of the kernel on the consistency of support vector machines”. Journal of machine learning research, Vol. 2, No. Nov, pp. 67–93, 2001.
[Tu 08] L. W. Tu. An Introduction to Manifolds. Springer, 2008.
[Vais 12] I. Vaisman. Lectures on the Geometry of Poisson Manifolds. Vol. 118, Birkhäuser, 2012.
[Van 12] C. Van Den Berg, J. P. R. Christensen, and P. Ressel. Harmonic Analysis on Semigroups: Theory of Positive Definite and Related Functions. Vol. 100, Springer Science & Business Media, 2012.
[Wein 83] A. Weinstein. “The local structure of Poisson manifolds”. Journal of differential geometry, Vol. 18, No. 3, pp. 523–557, 1983.
[Yuri 95] V. Yurinsky. Sums and Gaussian Vectors. Lecture Notes in Mathematics. Springer Berlin, Heidelberg, 1995.
[Zhou 08] D.-X. Zhou. “Derivative reproducing properties for kernel methods in learning theory”. Journal of Computational and Applied Mathematics, Vol. 220, No. 1-2, pp. 456–463, 2008.
[Zhu 20] A. Zhu, P. Jin, and Y. Tang. “Deep Hamiltonian networks based on symplectic integrators”. 2020.

	$\displaystyle\widehat{R}_{\lambda,N}(\widehat{h}_{\lambda,N}+h)$	$\displaystyle:=\frac{1}{N}\sum_{n=1}^{N}\\|X_{\widehat{h}_{\lambda,N}+h}(\mathbf{Z}^{(n)})-\mathbf{X}^{(n)}_{\sigma^{2}}\\|_{g}^{2}+\lambda\\|\widehat{h}_{\lambda,N}+h\\|_{\mathcal{H}_{K}}^{2}$
		$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\\|X_{\widehat{h}_{\lambda,N}}(\mathbf{Z}^{(n)})-\mathbf{X}^{(n)}_{\sigma^{2}}\\|_{g}^{2}+\lambda\left(\\|\widehat{h}_{\lambda,N}\\|_{\mathcal{H}_{K}}^{2}+\\|h\\|^{2}_{\mathcal{H}_{K}}+2\langle\widehat{h}_{\lambda,N},h\rangle_{\mathcal{H}_{K}}\right)$
		$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\\|X_{\widehat{h}_{\lambda,N}}(\mathbf{Z}^{(n)})-\mathbf{X}^{(n)}_{\sigma^{2}}\\|^{2}+\lambda\left(\\|\widehat{h}_{\lambda,N}\\|_{\mathcal{H}_{K}}^{2}+\\|h\\|^{2}_{\mathcal{H}_{K}}\right).$

	$\displaystyle\Big{\\|}\frac{1}{t}($	$\displaystyle K^{l}(\phi(t),\cdot)-K^{l}(y,\cdot))\Big{\\|}^{2}_{\mathcal{H}_{K}}=\frac{1}{t^{2}}\Big{(}K^{(l,l)}(\phi(t),\phi(t))-K^{(l,l)}(y,\phi(t))-K^{(l,l)}(\phi(t),y)+K^{(l,l)}(y,y)\Big{)}$
		$\displaystyle=\frac{1}{t^{2}}\int_{0}^{t}\Big{(}D^{(1,0)}K^{(l,l)}(\phi(s),\phi(t))\cdot\phi^{\prime}(s)-D^{(1,0)}K^{(l,l)}(\phi(s),y)\cdot\phi^{\prime}(s)\Big{)}\mathrm{d}s$
		$\displaystyle=\frac{1}{t^{2}}\int_{0}^{t}\int_{0}^{t}D^{(1,1)}K^{(l,l)}(\phi(s),\phi(\tau))\cdot(\phi^{\prime}(s),\phi^{\prime}(\tau))\mathrm{d}\tau\mathrm{d}s$
		$\displaystyle\leq\frac{1}{t^{2}}\\|K\\|_{C_{b}^{2(l+1)}}\int_{0}^{t}\int_{0}^{t}\\|\phi^{\prime}(\tau)\\|_{l}\\|\phi^{\prime}(s)\\|_{l}\mathrm{d}\tau\mathrm{d}s=\frac{1}{2}\\|K\\|_{C_{b}^{2(l+1)}},\quad\forall~{}\|t\|\leq 1,$

	$\displaystyle\\|Q_{N}h\\|_{\mathcal{H}_{K}}^{2}$	$\displaystyle=\frac{1}{N}\\|g_{N}(X_{h}(\mathbf{Z}_{N}),X_{K_{\cdot}}(\mathbf{Z}_{N}))\\|^{2}_{\mathcal{H}_{K}}=\frac{1}{N}\\|D^{(1,0)}K(\mathbf{Z}_{N},\cdot)\cdot(-JX_{h})(\mathbf{Z}_{N})\\|^{2}_{\mathcal{H}_{K}}$
		$\displaystyle=\frac{1}{N}D^{(1,1)}K(\mathbf{Z}_{N},\mathbf{Z}_{N})\cdot((-JX_{h})(\mathbf{Z}_{N}),(-JX_{h})(\mathbf{Z}_{N}))$
		$\displaystyle\leq\frac{1}{N}\\|D^{(1,1)}K\\|_{\infty}g_{N}(-JX_{h}(\mathbf{Z}_{N}),-JX_{h}(\mathbf{Z}_{N}))\leq\frac{C^{2}}{2N}\kappa^{2}g_{N}(X_{h}(\mathbf{Z}_{N}),X_{h}(\mathbf{Z}_{N}))$
		$\displaystyle\leq\frac{1}{2}C^{4}\kappa^{4}\\|h\\|^{2}_{\mathcal{H}_{K}}<\infty,$

	$\displaystyle\mathrm{Tr}(\Sigma_{N})$	$\displaystyle=\sum_{i=1}^{dN}g_{N}(f_{i},\Sigma_{N}f_{i})\leq\frac{1}{\lambda^{2}N^{2}}\sum_{i=1}^{dN}D^{(1,1)}K(\mathbf{Z}_{N},\mathbf{Z}_{N})(-\mathbb{J}_{N}f_{i},-\mathbb{J}_{N}f_{i})$
		$\displaystyle\leq\frac{1}{\lambda^{2}N^{2}}\sum_{i=1}^{dN}\\|D^{(1,1)}K\\|_{\infty}\cdot g_{N}(-\mathbb{J}_{N}f_{i},-\mathbb{J}_{N}f_{i})$
		$\displaystyle\leq\frac{1}{\lambda^{2}N^{2}}\\|K\\|_{C_{b}^{2}}\sum_{i=1}^{dN}g_{N}(-\mathbb{J}_{N}f_{i},-\mathbb{J}_{N}f_{i})$
		$\displaystyle\leq\frac{1}{\lambda^{2}N^{2}}\\|K\\|_{C_{b}^{2}}\sum_{i=1}^{dN}\left(\sum_{j=1}^{N}g(-Jf^{(j)}_{i},-Jf^{(j)}_{i})\right)$
		$\displaystyle\leq\frac{1}{\lambda^{2}N^{2}}\\|K\\|_{C_{b}^{2}}\sum_{i=1}^{dN}\left(\sum_{j=1}^{N}C^{2}\cdot g(f^{(j)}_{i},f^{(j)}_{i})\right)$
		$\displaystyle=\frac{C^{2}}{\lambda^{2}N^{2}}\\|K\\|_{C_{b}^{2}}\sum_{i=1}^{dN}g_{N}(f_{i},f_{i})=\frac{dC^{2}}{\lambda^{2}N}\\|K\\|_{C_{b}^{2}}=\frac{\kappa^{2}C^{2}d}{2\lambda^{2}N},$

A global structure-preserving kernel method for the learning of Poisson systems

Abstract

1 Introduction

Main Results.

2 Preliminaries on Poisson mechanics

2.1 Hamiltonian systems on Poisson manifolds

Definition 2.1.

Example 2.2 (Symplectic bracket).

Example 2.3 (Lie-Poisson bracket).

Definition 2.4.

Example 2.5 (Rigid body Casimirs [Mars 99]).

Example 2.6 (Underwater vehicle Casimirs [Leon 97]).

Definition 2.7.

Example 2.8 (The linear momentum).

Example 2.9 (The angular momentum).

Example 2.10 (Lifted actions on cotangent bundles).

2.2 High-order differentials and random variables on Riemannian manifolds

Definition 2.11 (Differentials).

Definition 2.12 (High-order differentials).

Remark 2.13.

Remark 2.14 (Functions with bounded differentials).

Definition 2.15 (Random variables on Riemannian manifold).

Remark 2.16 (Random variables on tangent spaces).

2.3 Compatible structures and local coordinate representations

Remark 2.17.

Proposition 2.18.

Proof.

Theorem 2.19 (Lie-Darboux-Weinstein coordinates).

3 Structure-preserving kernel regressions on Poisson manifolds

3.1 RKHS on Riemannian manifolds

Proposition 3.1.

Definition 3.2 (RKHS).

Theorem 3.3 (Differential reproducing property on Riemannian manifolds).

Remark 3.4.

Corollary 3.5.

Proof.

3.2 The structure-preserving learning problem

Remark 3.6.

3.3 Operator representations of the kernel estimator

Assumption 3.7.

Proposition 3.8.

Remark 3.9.

Proposition 3.10.

Proof.

Proposition 3.11.

3.4 The global solution of the learning problem

Proposition 3.12.

Theorem 3.13 (Differential Representer Theorem on Poisson manifolds).

Proof.

Remark 3.14 (On the uniqueness of the estimated Hamiltonian).

3.5 Symmetries and momentum conservation by kernel estimated Hamiltonians

Definition 3.15.

Proposition 3.16 (Noether’s Theorem for the structure-preserving kernel estimator).

Proof.

4 Estimation and approximation error analysis

Lemma 4.1.

Lemma 4.2 (Error bounds of noisy sampling error).

Assumption 4.3 (Source condition).

Proposition 4.4.

Proof.

Proposition 4.5.

Proof.

Theorem 4.6 (PAC bounds of the total reconstruction error).

Definition 4.7 (Coercivity condition).

Theorem 4.8 (Convergence upper rates for the total reconstruction error).

Lemma 4.9.

Proof.

Proposition 4.10 (From discrete data to continuous-time flows).

Proof.

5 Test examples and numerical experiments

5.1 Learning of Lie-Poisson systems

5.1.1 The rigid body

5.1.2 Underwater vehicle dynamics

5.1.3 Exact recovery of Hamiltonians in the RKHS

5.2 Learning on non-Euclidean symplectic manifolds

5.2.1 Metric, symplectic, and compatible structures in S2×S2S^{2}\times S^{2}

5.2.2 A universal kernel and data sampling on S2×S2S^{2}\times S^{2}

5.2.3 Exact recovery of Hamiltonians using a universal kernel

5.2.4 Recovery of two-vortex dynamics

6 Conclusion

5.2.1 Metric, symplectic, and compatible structures in $S^{2}\times S^{2}$

5.2.2 A universal kernel and data sampling on $S^{2}\times S^{2}$