# Learning Hamiltonian Dynamics with Reproducing Kernel Hilbert Spaces and
Random Features111This work was funded by the Norwegian Research Council under
Project Number 237896, SFI Offshore Mechatronics.

Torbjørn Smith<EMAIL_ADDRESS>Olav Egeland<EMAIL_ADDRESS>
###### Abstract

A method for learning Hamiltonian dynamics from a limited and noisy dataset is
proposed. The method learns a Hamiltonian vector field on a reproducing kernel
Hilbert space (RKHS) of inherently Hamiltonian vector fields, and in
particular, odd Hamiltonian vector fields. This is done with a symplectic
kernel, and it is shown how the kernel can be modified to an odd symplectic
kernel to impose the odd symmetry. A random feature approximation is developed
for the proposed kernel to reduce the problem size. This includes random
feature approximations for odd kernels. The performance of the method is
validated in simulations for three Hamiltonian systems. It is demonstrated
that the use of an odd symplectic kernel improves prediction accuracy, and
that the learned vector fields are Hamiltonian and exhibit the imposed odd
symmetry characteristics.

###### keywords:

Machine learning , System identification , Reproducing kernel Hilbert space

††journal: European Journal of Control

[aff1] organization = Department of Mechanical and Industrial Engineering,
Norwegian University of Science and Technology (NTNU), addressline = Richard
Birkelands vei 2B, city = Trondheim, postcode = 7491, country = Norway

##  1 Introduction

Learning dynamical systems is an important area of research in robotics and
control engineering, and data-driven methods have emerged as a robust approach
for system identification, where classical analytical methods may be
impractical. The aim is to utilize machine learning to derive a model of the
underlying dynamical system from a set of measurement data [1]. The efficiency
of data-driven methods is dependent on the quality of the training dataset.
Notably, these methods may encounter challenges such as limited generalization
beyond the provided dataset [2] and may be susceptible to overfitting in cases
where the data set is limited or noisy [3]. Assembling a viable dataset may
also be labor intensive or even impractical in real-world or online
applications. Furthermore, as the data set grows, the computational cost of
learning the model increases, as does the inference time of the final learned
model for some learning methods [2] [4]. Ensuring that the learned models are
stable and robust is vital, particularly for safety-critical and control
applications [5]. In response to these challenges, researchers have developed
different strategies to guide or restrict the learning of dynamical systems
using prior information, yielding satisfactory results even with limited
datasets.

###  1.1 Related work

In [6], financial price prediction was explored using a data-driven approach
with functions in a reproducing kernel Hilbert space (RKHS). By designing an
odd reproducing kernel which imposed an odd symmetry constraint on the price
action, the learned model demonstrated improved prediction accuracy and
reduced overfitting compared to the unconstrained method.

Inverse dynamics was learned in [7] by capturing the system’s Lagrangian as a
function in a RKHS. Polynomial basis functions were used in [8] to investigate
control-oriented learning of Lagrangian and Hamiltonian systems. It was
demonstrated that accurate and generalized learning from a limited number of
trajectories could be achieved by learning these functions. The learning-from-
demonstrations problem was addressed in [2] where the focus was on copying
human demonstrations using a data-driven approach. A dynamical system was
learned in an RKHS with random Fourier features for dimensionality reduction,
allowing the imitation of human-drawn shapes. The learned dynamical system
included desired equilibrium points, and point-wise contraction constraints
were enforced along the trajectory to create a contraction region around the
desired path, conditioning the learned vector field.

The work in [9] focused on learning the Hamiltonian dynamics of energy-
conserving systems using neural networks, by learning the Hamiltonian as a
parametric function. This approach significantly enhanced the predictive
accuracy of the learned system. Building upon this, [10] further refined the
method by eliminating the need for higher-order derivatives of the generalized
coordinate and incorporating the option for energy-based control. In [11], the
work in [9] was extended by using the symplectic Leapfrog integrator to
integrate the partial derivatives of the learned Hamiltonian. The loss was
then back-propagated through the integrator over multiple time steps,
resulting in improved learning performance for more complex and noisy
Hamiltonian systems.

Learning nonlinear dynamics with a stabilizability constraint was investigated
in [4]. The dynamics were learned using a contraction constraint, and the
model was evaluated using a planar drone. The method enhanced trajectory
generation, tracking, and data efficiency. The model was learned in an RKHS,
and utilized random Fourier features for dimensionality reduction. In [12],
they performed nonlinear system identification by incorporating constraints
enforcing prior knowledge of the region of attraction. The stability region
was enforced using a Lyapunov function, and the hypothesis space for the
learned model was an RKHS. [13] explores learning dynamical systems in an
RKHS, incorporating a bias term in the regularized least squares cost to embed
prior knowledge, improving data efficiency and out-of-sample generalization.
In [14] the identification of nonlinear input-output operators in an RKHS is
studied. Nonexpansive operators are introduced to identify operators that
satisfy a wide range of dissipativity and integral quadratic constraints.

A method for optimization-based learning of vector fields with side
constraints from a limited number of trajectories was presented in [3]. The
method utilized polynomial basis functions and the learning problem was solved
using semidefinite programming. They showed how the side constraints improved
learning accuracy and model generalization. The considered side constraints
included interpolation at a finite set of points, sign symmetry, gradient and
Hamiltonian dynamics, coordinate non-negativity, directional monotonicity, and
invariance of sets.

###  1.2 Contribution

In this paper, we show how Hamiltonian dynamical systems with odd vector
fields can be learned in an RKHS by selecting a kernel that ensures that the
learned vector fields are Hamiltonian and odd symmetric. The proposed kernel
is approximated using random Fourier features for dimensionality reduction. We
also include a novel approximation of the odd and even kernels using random
Fourier features. Encoding the constraints in the kernel improves learning
time as the straightforward closed-form solution of the learning problem is
retained. Three illustrative simulation examples demonstrate that the
generalization properties of the learned model when using out-of-sample data
points are improved through the additional constraints, and that energy
preservation and odd symmetry are encoded in the final model.

A preliminary version of the proposed learning algorithm was presented in
[15]. In the present paper, we extend the method in [15] by incorporating
random Fourier features to approximate the proposed kernel. Furthermore, the
simulation experiments are expanded to encompass more sophisticated
Hamiltonian systems that are common in the system identification literature.

The paper is organized as follows: Section 2 outlines the problem addressed in
this work. Section 3 provides a review of the relevant theory related to
reproducing kernel Hilbert spaces, learning dynamical systems, random Fourier
features, and Hamiltonian mechanics. The proposed method is detailed in
Section 4, where an odd symplectic kernel is developed and approximated using
random features. Section 6 presents the numerical experiments used to validate
the proposed method. Finally, Section 7 presents the conclusion and future
work.

##  2 Problem formulation

This paper explores learning the dynamics of an unknown system from limited
data. The system dynamics are given by the vector field

$\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})$ (1)

where ${\boldsymbol{x}\in{\mathbb{R}}^{n}}$ is the state vector,
${\dot{\boldsymbol{x}}\in{\mathbb{R}}^{n}}$ is the time derivative of the
state vector, and
${\boldsymbol{f}:{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}^{n}}$ are the system
dynamics. It is assumed that $\boldsymbol{y}=\dot{\boldsymbol{x}}$ is
available as a measurement or from numerical differentiation. Given a set of
${N}$ data points
${\\{(\boldsymbol{x}_{i},\boldsymbol{y}_{i})\in{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\\}_{i=1}^{N}}$
from simulations or measurements, the aim is to learn a function
${\boldsymbol{f}^{*}\in{\mathcal{F}}}$, where ${\mathcal{F}}$ is a class of
functions. The class of functions ${\mathcal{F}}$ will be the functions of a
reproducing kernel Hilbert space (RKHS) determined by a reproducing kernel
[16]. The function $\boldsymbol{f}^{*}$ is found by the regularized
minimization problem [17]

$\boldsymbol{f}^{*}=\operatorname*{\arg\min}_{\boldsymbol{f}\in{\mathcal{F}}}\frac{1}{N}\sum_{i=1}^{N}\|\boldsymbol{f}(\boldsymbol{x}_{i})-\boldsymbol{y}_{i}\|^{2}+\lambda\|\boldsymbol{f}\|^{2}_{{\mathcal{F}}}$
(2)

where ${\lambda>0}$ is the regularization parameter. It is well known that
this approach may lead to inaccurate generalization beyond the data set used
to learn the dynamical model. Furthermore, if the trajectories in the data set
are limited and noisy, the learned dynamical model may fail to capture the
dynamics of the underlying system due to overfitting.

It is assumed that there is some information about the physical properties of
the dynamical system. This type of side information about the system was
treated in [3] where the function class ${\mathcal{F}}$ was polynomial
functions. The additional information about the dynamics was included as side
constraints in [3] by defining a subset
${\mathcal{S}}_{i}\subset{\mathcal{F}}$ for each side constraint $i$, so that
the function $\boldsymbol{f}^{*}$ satisfies the side constraint whenever
$\boldsymbol{f}^{*}\in{\mathcal{S}}_{i}$. The learning problem including the
side constraints can be formulated as

$\boldsymbol{f}^{*}=\operatorname*{\arg\min}_{\boldsymbol{f}\in{\mathcal{F}}\cap{\mathcal{S}}_{1}\cap\dots\cap{\mathcal{S}}_{k}}\frac{1}{N}\sum_{i=1}^{N}\|\boldsymbol{f}(\boldsymbol{x}_{i})-\boldsymbol{y}_{i}\|^{2}$
(3)

In this paper the side constraints are instead handled by defining a
reproducing kernel which ensures that the RKHS function class ${\mathcal{F}}$
inherently satisfies the relevant side constraints. It is well-known that this
can be done to have a RKHS where the vector field $\boldsymbol{f}^{*}$ is
curl-free or divergence-free [18], symplectic [19], odd or even [6]. It is
also possible to impose additional side constraints like contraction [2] or
stabilizability [4] along the trajectories of the dataset, but this will not
be addressed in this paper.

In this paper the function class ${\mathcal{F}}$ will be a reproducing kernel
Hilbert space (RKHS). The side constraints are that the state dynamics are
symplectic, and, in addition, odd in the sense that
$\boldsymbol{f}(-\boldsymbol{x})=-\boldsymbol{f}(\boldsymbol{x})$, and this is
ensured by selecting an appropriate reproducing kernel.

##  3 Preliminaries

###  3.1 Reproducing kernel Hilbert space

This paper deals with the learning of functions in a reproducing kernel
Hilbert space (RKHS) [16]. Learning real-valued scalar functions
$f:{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}$ is performed by assuming
$f\in\mathcal{H}_{k}$, where $\mathcal{H}_{k}$ is an RKHS defined in terms of
a reproducing kernel $k$ which is a real-valued scalar function
${k:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}}$. Given a
sample of $N$ data-points $\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{N}$ the
function value of a function $f\in\mathcal{H}_{k}$ is given by

$f(\boldsymbol{x})=\sum_{i=1}^{N}k(\boldsymbol{x},\boldsymbol{x}_{i})a_{i}\in{\mathbb{R}}$
(4)

where $a_{1},\ldots,a_{N}\in{\mathbb{R}}$. Let the function
$k_{x}:{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}$ be defined by
$k_{x}(\boldsymbol{z})=k(\boldsymbol{x},\boldsymbol{z})$ for all
$\boldsymbol{x},\boldsymbol{z}\in{\mathbb{R}}^{n}$. Then the $\mathcal{H}_{k}$
is defined as the closure of

$\operatorname{span}\\{k_{x}:k_{x}(\boldsymbol{z})=k(\boldsymbol{x},\boldsymbol{z})\\}\subseteq\mathcal{H}_{k}$
(5)

and the reproducing property is $f(\boldsymbol{x})={\left\langle
f,k_{x}\right\rangle}_{\mathcal{H}_{k}}$ where
${\left\langle\cdot,\cdot\right\rangle}_{\mathcal{H}_{k}}$ is the inner
product in $\mathcal{H}_{k}$. It is noted that the reproducing property
implies that

${\left\langle
k_{x},k_{z}\right\rangle}_{\mathcal{H}_{k}}=k_{x}(\boldsymbol{z})=k(\boldsymbol{x},\boldsymbol{z})$
(6)

According to the Moore-Aronszajn theorem [16] the reproducing kernel ${k}$ is
a positive definite kernel, which means that
$k(\boldsymbol{x},\boldsymbol{z})=k(\boldsymbol{z},\boldsymbol{x})$ for all
$\boldsymbol{x},\boldsymbol{z}\in{\mathbb{R}}^{n}$, and

$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}c_{i}c_{j}k(\boldsymbol{x}_{i},\boldsymbol{x}_{j})$
$\displaystyle\geq 0$ (7)

for any set of points ${\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{N}\in X}$ and
any set of real numbers ${c_{1},\dots,c_{N}\in{\mathbb{R}}}$.

The extension of RKHS to vector-valued functions
${\boldsymbol{f}:{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}^{m}}$ is found in
[17]. Then the RKHS $\mathcal{H}_{K}$ with inner product
${\left\langle\cdot,\cdot\right\rangle}_{\mathcal{H}_{K}}$ is defined for the
reproducing kernel
$\boldsymbol{K}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\rightarrow{\mathbb{R}^{m\times
m}}$, and a function $\boldsymbol{f}\in\mathcal{H}_{K}$ is given by

$\boldsymbol{f}(\boldsymbol{x})=\sum_{i=1}^{N}\boldsymbol{K}(\boldsymbol{x},\boldsymbol{x}_{i})\boldsymbol{a}_{i}\in{\mathbb{R}}^{m}$
(8)

where $\boldsymbol{a}_{1},\ldots,\boldsymbol{a}_{N}\in{\mathbb{R}}^{m}$.
Define the function
${\boldsymbol{K}_{x}\boldsymbol{y}:{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}^{m}}$
as the function which for all $\boldsymbol{f}\in\mathcal{H}_{K}$ satisfies the
reproducing property

${\left\langle\boldsymbol{f}(\boldsymbol{x}),\boldsymbol{y}\right\rangle}={\left\langle\boldsymbol{f},\boldsymbol{K}_{x}\boldsymbol{y}\right\rangle}_{\mathcal{H}_{K}}$
(9)

Then the reproducing kernel $\boldsymbol{K}(\boldsymbol{x},\boldsymbol{z})$ is
defined for every $\boldsymbol{y}$ by
$(\boldsymbol{K}_{z}\boldsymbol{y})(\boldsymbol{x})=\boldsymbol{K}(\boldsymbol{x},\boldsymbol{z})\boldsymbol{y}\in{\mathbb{R}}^{m}$
and ${\mathcal{H}_{K}}$ is the closure of
$\operatorname{span}\\{\boldsymbol{K}_{x}\boldsymbol{y}\\}$.

The reproducing kernel is positive definite in the sense that
$\boldsymbol{K}(\boldsymbol{x},\boldsymbol{z})=\boldsymbol{K}(\boldsymbol{z},\boldsymbol{x})^{\text{T}}$
for all $\boldsymbol{x},\boldsymbol{z}\in{\mathbb{R}}^{n}$, and, in addition,

$\sum_{i=1}^{N}\sum_{j=1}^{N}{\left\langle\boldsymbol{y}_{i},\boldsymbol{K}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})\boldsymbol{y}_{j}\right\rangle}\geq
0$ (10)

for every set of vectors
${\\{\boldsymbol{x}_{i}\\}_{i=1}^{N}\in{\mathbb{R}}^{n}}$ and
${\\{\boldsymbol{y}_{i}\\}_{i=1}^{N}\in{\mathbb{R}}^{m}}$. The reproducing
property implies that

${\left\langle\boldsymbol{K}_{z}\boldsymbol{w},\boldsymbol{K}_{x}\boldsymbol{y}\right\rangle}_{\mathcal{H}_{K}}={\left\langle\boldsymbol{K}_{z}\boldsymbol{w}(\boldsymbol{x}),\boldsymbol{y}\right\rangle}={\left\langle\boldsymbol{K}(\boldsymbol{x},\boldsymbol{z})\boldsymbol{w},\boldsymbol{y}\right\rangle}$
(11)

Then functions in the RKHS $\mathcal{H}_{K}$ can be defined as

$\boldsymbol{f}=\sum_{i=1}^{N}\boldsymbol{K}_{x_{i}}\boldsymbol{y}_{i}\in\mathcal{H}_{K},\quad\boldsymbol{g}=\sum_{j=1}^{N}\boldsymbol{K}_{z_{j}}\boldsymbol{w}_{j}\in\mathcal{H}_{K}$
(12)

with inner product

${\left\langle\boldsymbol{f},\boldsymbol{g}\right\rangle}_{\mathcal{H}_{K}}=\sum_{i=1}^{N}\sum_{j=1}^{N}{\left\langle\boldsymbol{y}_{i},\boldsymbol{K}(\boldsymbol{x}_{i},\boldsymbol{z}_{j})\boldsymbol{w}_{j}\right\rangle}_{\mathcal{H}_{K}}$
(13)

and norm
$\|\boldsymbol{f}\|_{\mathcal{H}_{K}}^{2}={\left\langle\boldsymbol{f},\boldsymbol{f}\right\rangle}_{\mathcal{H}_{K}}$.

It is noted that (9) gives
${\left\langle\boldsymbol{f}(\boldsymbol{x}),\boldsymbol{y}\right\rangle}={\left\langle\boldsymbol{K}_{x}^{*}\boldsymbol{f},\boldsymbol{y}\right\rangle}$,
where $\boldsymbol{K}_{x}^{*}$ is the adjoint of $\boldsymbol{K}_{x}$. It
follows that
$\boldsymbol{f}(\boldsymbol{x})=\boldsymbol{K}_{x}^{*}\boldsymbol{f}$, and

$\|\boldsymbol{f}(\boldsymbol{x})\|=\|\boldsymbol{K}_{x}^{*}\boldsymbol{f}\|\leq\|\boldsymbol{K}_{x}^{*}\|\|\boldsymbol{f}\|_{\mathcal{H}_{K}}\leq\sqrt{\|\boldsymbol{K}(\boldsymbol{x},\boldsymbol{x})\|}\|\boldsymbol{f}\|_{\mathcal{H}_{K}}$
(14)

where $\sqrt{\|\boldsymbol{K}(\boldsymbol{x},\boldsymbol{x})\|}$ is bounded
[17]. This implies

$\|\boldsymbol{f}(\boldsymbol{x})-\boldsymbol{g}(\boldsymbol{x})\|\leq\sqrt{\|\boldsymbol{K}(\boldsymbol{x},\boldsymbol{x})\|}\|\boldsymbol{f}-\boldsymbol{g}\|_{\mathcal{H}_{K}}$
(15)

which means that if ${\|\boldsymbol{f}-\boldsymbol{g}\|_{\mathcal{H}_{K}}}$
converges to zero, then
${\|\boldsymbol{f}(\boldsymbol{x})-\boldsymbol{g}(\boldsymbol{x})}\|$
converges to zero for each ${\boldsymbol{x}}$.

###  3.2 Learning dynamical systems with RKHS

In this paper we aim to learn vector fields given by (1) where
$\boldsymbol{y}=\dot{\boldsymbol{x}}$ is available. The estimation of the
vector field is done using the vector-valued regularized least-squares problem
[17]

$\boldsymbol{f}^{*}=\operatorname*{\arg\min}_{\boldsymbol{f}\in\mathcal{H}_{K}}\frac{1}{N}\sum_{i=1}^{N}\|\boldsymbol{f}(\boldsymbol{x}_{i})-\boldsymbol{y}_{i}\|^{2}+\lambda\|\boldsymbol{f}\|_{\mathcal{H}_{K}}^{2}$
(16)

where
${{\mathcal{Z}}=\\{(\boldsymbol{x}_{i},\boldsymbol{y}_{i})\in{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\\}_{i=1}^{N}}$
is the data used to learn the vector field, and $\lambda>0$ is the
regularization parameter. The optimal solution is [17]

$\boldsymbol{f}^{*}=\sum_{i=1}^{N}\boldsymbol{K}(\cdot,\boldsymbol{x}_{i})\boldsymbol{a}_{i}\in\mathcal{H}_{K}$
(17)

where the coefficients ${\boldsymbol{a}_{i}\in\mathbb{R}^{m}}$ are the unique
solutions of

$\sum_{j=1}^{N}\boldsymbol{K}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})\boldsymbol{a}_{j}+N\lambda\boldsymbol{a}_{i}=\boldsymbol{y}_{i}$
(18)

The function value of the optimal vector field is

$\boldsymbol{f}^{*}(\boldsymbol{x})=\sum_{i=1}^{N}\boldsymbol{K}(\boldsymbol{x},\boldsymbol{x}_{i})\boldsymbol{a}_{i}\in{\mathbb{R}}^{n}$
(19)

A matrix formulation of (18) is found in [20].

###  3.3 Random Fourier features

As the number of trajectories and samples in the training set increases, the
size of the learning problem (16) increases, increasing the training time for
the model coefficients. Furthermore, the inference time of the learned model
(19) increases. To limit training and inference time random features are used
to approximate the kernel functions associated with functions in an RKHS [4].

Let the matrix-valued kernel function
${\boldsymbol{K}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\rightarrow{\mathbb{R}^{m\times
m}}}$ be approximated by the matrix-valued feature map
${\boldsymbol{\Psi}:{\mathbb{R}}^{n}\rightarrow{\mathbb{R}^{d\times m}}}$ [2]
as follows

$\boldsymbol{K}(\boldsymbol{x},\boldsymbol{y})\approx\boldsymbol{\Psi}(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\Psi}(\boldsymbol{y})$
(20)

Then the function value of the vector field in (19) can be parameterized as

$\boldsymbol{f}(\boldsymbol{x})=\sum_{i=1}^{N}\boldsymbol{K}(\boldsymbol{x},\boldsymbol{x}_{i})\boldsymbol{a}_{i}\approx\sum_{i=1}^{N}\boldsymbol{\Psi}(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\Psi}(\boldsymbol{x}_{i})\boldsymbol{c}_{i}=\boldsymbol{\Psi}(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\alpha}$
(21)

where ${\boldsymbol{c}_{i}\in\mathbb{R}^{m}}$, and the new model coefficients
are

$\boldsymbol{\alpha}=\sum_{i=1}^{N}\boldsymbol{\Psi}(\boldsymbol{x}_{i})\boldsymbol{c}_{i}\in{\mathbb{R}}^{d}$
(22)

This makes ${\boldsymbol{f}(\boldsymbol{x})}$ a linear combination of ${d}$
random feature functions, instead of ${N}$ kernel functions. The number of
random features ${d}$ is chosen by balancing the quality of the approximation
and the computational time of the function. Using (13) and (21) the norm of
${\boldsymbol{f}\in\mathcal{H}_{K}}$ can be written as

$\|\boldsymbol{f}\|_{\mathcal{H}_{K}}^{2}=\|\boldsymbol{\alpha}\|^{2}$ (23)

The vector-valued regularized least-squares problem in (16) is then re-
formulated as an optimization problem over
${\boldsymbol{\alpha}\in{\mathbb{R}}^{d}}$

$\min_{\boldsymbol{\alpha}\in{\mathbb{R}}^{d}}\frac{1}{N}\sum_{i=1}^{N}\|\boldsymbol{\Psi}(\boldsymbol{x}_{i})^{{\text{T}}}\boldsymbol{\alpha}-\boldsymbol{y}_{i}\|^{2}+\lambda\boldsymbol{\alpha}^{{\text{T}}}\boldsymbol{\alpha}$
(24)

By defining

$\boldsymbol{\Phi}=\begin{bmatrix}\boldsymbol{\Psi}(\boldsymbol{x}_{1})^{{\text{T}}}\\\
\vdots\\\
\boldsymbol{\Psi}(\boldsymbol{x}_{N})^{{\text{T}}}\end{bmatrix}\in{\mathbb{R}}^{mN\times
d}\;,\qquad\boldsymbol{Y}=\begin{bmatrix}\boldsymbol{y}_{1}\\\ \vdots\\\
\boldsymbol{y}_{N}\end{bmatrix}\in{\mathbb{R}}^{mN}$ (25)

the optimization problem in (24) can be written as

$\min_{\boldsymbol{\alpha}\in{\mathbb{R}}^{d}}\frac{1}{N}\left(\boldsymbol{\Phi}\boldsymbol{\alpha}-\boldsymbol{Y}\right)^{{\text{T}}}\left(\boldsymbol{\Phi}\boldsymbol{\alpha}-\boldsymbol{Y}\right)+\lambda\boldsymbol{\alpha}^{{\text{T}}}\boldsymbol{\alpha}$
(26)

where the minimum is found when the gradient with respect to
${\boldsymbol{\alpha}}$ is zero

$\boldsymbol{\alpha}^{*}={\left(\boldsymbol{\Phi}^{{\text{T}}}\boldsymbol{\Phi}+\lambda
N{\boldsymbol{I}}\right)}^{-1}\boldsymbol{\Phi}^{{\text{T}}}\boldsymbol{Y}$
(27)

and the function value of the optimal vector field is

$\boldsymbol{f}^{*}(\boldsymbol{x})=\boldsymbol{\Psi}(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\alpha}^{*}$
(28)

###  3.4 Hamiltonian dynamics

Consider a holonomic system with generalized coordinates
$\boldsymbol{q}=[q_{1},\ldots,q_{n}]^{\text{T}}$ and Lagrangian [21]

$L(\boldsymbol{q},\dot{\boldsymbol{q}},t)=T(\boldsymbol{q},\dot{\boldsymbol{q}},t)-U(\boldsymbol{q})$
(29)

and the Lagrangian equation of motion

$\frac{d}{dt}\left(\frac{\partial
L(\boldsymbol{q},\dot{\boldsymbol{q}},t)}{\partial\dot{\boldsymbol{q}}}\right)-\frac{\partial
L(\boldsymbol{q},\dot{\boldsymbol{q}},t)}{\partial\boldsymbol{q}}=\boldsymbol{\tau}$
(30)

The momentum vector $\boldsymbol{p}=[p_{1},\ldots,p_{n}]^{\text{T}}$ is
defined by

$\boldsymbol{p}=\frac{\partial
L(\boldsymbol{q},\dot{\boldsymbol{q}},t)}{\partial\dot{\boldsymbol{q}}}^{{\text{T}}}$
(31)

A change of coordinate from ${(\boldsymbol{q},\dot{\boldsymbol{q}})}$ to
${(\boldsymbol{q},\boldsymbol{p})}$ is introduced with the Legendre
transformation [21]

$H(\boldsymbol{q},\boldsymbol{p},t)=\boldsymbol{p}^{{\text{T}}}\boldsymbol{\phi}(\boldsymbol{q},\boldsymbol{p},t)-L(\boldsymbol{q},\boldsymbol{\phi}(\boldsymbol{q},\boldsymbol{p},t),t)$
(32)

where
${\dot{\boldsymbol{q}}=\boldsymbol{\phi}(\boldsymbol{q},\boldsymbol{p},t)}$.
Hamilton’s equations of motion then follow as

$\displaystyle\dot{\boldsymbol{q}}$ $\displaystyle=\frac{\partial
H}{\partial\boldsymbol{p}}^{{\text{T}}}$ (33)
$\displaystyle\dot{\boldsymbol{p}}$ $\displaystyle=-\frac{\partial
H}{\partial\boldsymbol{q}}^{{\text{T}}}+\boldsymbol{\tau}$ (34)

The time derivative of the Hamiltonian is seen to be

$\frac{dH}{dt}=\frac{\partial
H}{\partial\boldsymbol{p}}\dot{\boldsymbol{p}}+\frac{\partial
H}{\partial\boldsymbol{q}}\dot{\boldsymbol{q}}+\frac{\partial H}{\partial
t}=\dot{\boldsymbol{q}}^{{\text{T}}}\boldsymbol{\tau}+\frac{\partial
H}{\partial t}$ (35)

If the Hamiltonian does not depend on time ${t}$ and
${\boldsymbol{\tau}=\boldsymbol{0}}$, then

$\frac{dH(\boldsymbol{q},\boldsymbol{p})}{dt}=0$ (36)

The numerical value of the Hamiltonian $H$ will depend on the definition of
the zero level of the potential $U(\boldsymbol{q})$.

###  3.5 Symplectic property of Hamiltonian dynamics

The phase space is defined as a ${2n}$-dimensional space with state vector
${\boldsymbol{x}=\begin{bmatrix}\boldsymbol{q}^{{\text{T}}},\boldsymbol{p}^{{\text{T}}}\end{bmatrix}^{{\text{T}}}\in{\mathbb{R}^{2n}}}$.
The Hamiltonian dynamics (33), (34) with $\boldsymbol{\tau}={\boldsymbol{0}}$
can be given in the phase space as

$\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})=\boldsymbol{J}\nabla
H(\boldsymbol{x})$ (37)

where

$\boldsymbol{J}=\begin{bmatrix}\boldsymbol{0}&{\boldsymbol{I}}\\\
-{\boldsymbol{I}}&\boldsymbol{0}\end{bmatrix}\in{\mathbb{R}}^{2n\times 2n}$
(38)

is the symplectic matrix. The flow in the phase space with initial condition
$\boldsymbol{x}(0)=\boldsymbol{x}_{0}$ is given by

$\boldsymbol{\phi}_{t}(\boldsymbol{x}_{0})=\boldsymbol{x}(t)$ (39)

where ${\boldsymbol{x}(t)}$ is the solution of the Hamiltonian dynamics (37)
at time ${t}$ given the initial value ${\boldsymbol{x}_{0}}$. Define the
Jacobian
${\boldsymbol{\Psi}(t)=\partial\boldsymbol{\phi}_{t}(\boldsymbol{x}_{0})/\partial\boldsymbol{x}_{0}}$.
Then the system (37) is said to be symplectic if

$\boldsymbol{\Psi}(t)^{{\text{T}}}\boldsymbol{J}\boldsymbol{\Psi}(t)=\boldsymbol{J}$
(40)

for all $t\geq 0$. The system (37) is Hamiltonian if and only if it is
symplectic [22].

##  4 Reproducing kernels

###  4.1 Gaussian kernel

The Gaussian reproducing kernel $k_{\sigma}$ is widely used in applications.
In the scalar case it is defined as the shift-invariant function

$k_{\sigma}(\boldsymbol{x},\boldsymbol{z})=g_{\sigma}(\boldsymbol{x}-\boldsymbol{z})=e^{-\frac{\|\boldsymbol{x}-\boldsymbol{z}\|^{2}}{2\sigma^{2}}}$
(41)

where $\sigma>0$ and the function $g_{\sigma}$ is called the signature of the
shift-invariant kernel [18]. It is noted that the Gaussian kernel satisfies
$k_{\sigma}(\boldsymbol{x},\boldsymbol{z})=k_{\sigma}(-\boldsymbol{x},-\boldsymbol{z})$.
In the case of RKHS for vector-valued functions a separable Gaussian kernel
$\boldsymbol{K}_{\sigma}(\boldsymbol{x},\boldsymbol{z})=\boldsymbol{G}_{\sigma}(\boldsymbol{x}-\boldsymbol{z})$
can be used where [2]

$\boldsymbol{G}_{\sigma}(\boldsymbol{x})=g_{\sigma}(\boldsymbol{x}){\boldsymbol{I}}_{n}\quad\in{\mathbb{R}^{n\times
n}}$ (42)

###  4.2 Curl-free kernel

In the learning of a vector field
$\boldsymbol{f}(\boldsymbol{x})\in{\mathbb{R}}^{n}$ a curl-free reproducing
kernel
$\boldsymbol{K}_{c}(\boldsymbol{x},\boldsymbol{z})=\boldsymbol{G}_{c}(\boldsymbol{x}-\boldsymbol{z})\in{\mathbb{R}}^{n\times
n}$ can be derived from a shift-invariant reproducing kernel
${k(\boldsymbol{x},\boldsymbol{z})=g(\boldsymbol{x}-\boldsymbol{z})\in{\mathbb{R}}}$,
where $k$ is typically the Gaussian kernel [23]. The signature of the curl-
free kernel is given by

$\boldsymbol{G}_{c}(\boldsymbol{x})=-\nabla\nabla^{{\text{T}}}g(\boldsymbol{x})=\begin{bmatrix}-\nabla\frac{\partial
g(\boldsymbol{x})}{\partial x_{1}}\;\dots\;-\nabla\frac{\partial
g(\boldsymbol{x})}{\partial x_{n}}\end{bmatrix}$ (43)

where column $i$ is the gradient of the scalar field $\frac{\partial
g(\boldsymbol{x})}{\partial x_{i}}$. The curl of the gradient of a scalar
field is always zero, which implies that each column of ${\boldsymbol{G}_{c}}$
is curl-free. A function of the RKHS of $\boldsymbol{K}_{c}$ is given by

$\boldsymbol{f}(\boldsymbol{x})=\sum_{i=1}^{N}\boldsymbol{G}_{c}(\boldsymbol{x}-\boldsymbol{x}_{i})\boldsymbol{a}_{i}$
(44)

where $\boldsymbol{f}(\boldsymbol{x})$ is a linear combination of the $n$
columns of the $N$ matrices
$\boldsymbol{G}_{c}(\boldsymbol{x}-\boldsymbol{x}_{i})$. Since each of these
$nN$ columns is curl-free it follows that the vector field
$\boldsymbol{f}(\boldsymbol{x})$ is curl-free. If $k$ is selected as the
Gaussian kernel, then the curl-free kernel is given by [23]

$\boldsymbol{G}_{c}(\boldsymbol{x})=-\nabla\nabla^{{\text{T}}}g_{\sigma}(\boldsymbol{x})=\frac{1}{\sigma^{2}}e^{-\frac{\boldsymbol{x}^{{\text{T}}}\boldsymbol{x}}{2\sigma^{2}}}\left({\boldsymbol{I}}-\frac{\boldsymbol{x}\boldsymbol{x}^{{\text{T}}}}{\sigma^{2}}\right)$
(45)

which is the curl-free kernel used in [2].

###  4.3 Symplectic kernel

In [19] a symplectic kernel was presented for adaptive prediction for
Hamiltonian dynamics. The symplectic kernel is based on the curl-free kernel
(43), and is given by

$\boldsymbol{K}_{s}(\boldsymbol{x},\boldsymbol{z})=\boldsymbol{G}_{s}(\boldsymbol{x}-\boldsymbol{z})=\boldsymbol{J}\boldsymbol{G}_{c}(\boldsymbol{x}-\boldsymbol{z})\boldsymbol{J}^{{\text{T}}}$
(46)

To verify that functions in the resulting RKHS is symplectic it is used that
any function in the RKHS is given by

$\boldsymbol{f}(\boldsymbol{x})=\sum_{i=1}^{N}\boldsymbol{G}(\boldsymbol{x}-\boldsymbol{x}_{i})\boldsymbol{a}_{i}$
(47)

which gives

$\boldsymbol{f}(\boldsymbol{x})=-\sum_{i=1}^{N}\boldsymbol{J}\nabla\nabla^{{\text{T}}}g(\boldsymbol{x}-\boldsymbol{x}_{i})\boldsymbol{J}^{{\text{T}}}\boldsymbol{a}_{i}=-\boldsymbol{J}\nabla\sum_{i=1}^{N}\nabla^{{\text{T}}}g(\boldsymbol{x}-\boldsymbol{x}_{i})\boldsymbol{c}_{i}$
(48)

where ${\boldsymbol{c}_{i}=\boldsymbol{J}^{{\text{T}}}\boldsymbol{a}_{i}}$.
This results in the Hamiltonian dynamics

$\boldsymbol{f}(\boldsymbol{x})=\boldsymbol{J}\nabla H(\boldsymbol{x})$ (49)

where the Hamiltonian is

${H}(\boldsymbol{x})=-\sum_{i=1}^{N}\nabla^{{\text{T}}}g(\boldsymbol{x}-\boldsymbol{x}_{i})\boldsymbol{c}_{i}$
(50)

In [19], the Gaussian kernel (41) was used. Using the curl-free kernel derived
from the Gaussian kernel in (45), the symplectic kernel is written as

$\boldsymbol{G}_{s}(\boldsymbol{x})=\frac{1}{\sigma^{2}}e^{-\frac{\boldsymbol{x}^{{\text{T}}}\boldsymbol{x}}{2\sigma^{2}}}\boldsymbol{J}\left({\boldsymbol{I}}-\frac{\boldsymbol{x}\boldsymbol{x}^{{\text{T}}}}{\sigma^{2}}\right)\boldsymbol{J}^{{\text{T}}}$
(51)

###  4.4 Odd kernel

Consider a shift-invariant reproducing kernel
$k(\boldsymbol{x},\boldsymbol{z})=g(\boldsymbol{x}-\boldsymbol{z})$ which
satisfies
$k(\boldsymbol{x},\boldsymbol{z})=k(-\boldsymbol{x},-\boldsymbol{z})$. Then
$k(-\boldsymbol{x},\boldsymbol{z})=k(\boldsymbol{x},-\boldsymbol{z})$, and an
odd kernel can be defind as

$k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})=\frac{1}{2}(k(\boldsymbol{x},\boldsymbol{z})-k(-\boldsymbol{x},\boldsymbol{z}))\quad\in{\mathbb{R}}$
(52)

which is odd in the variable $\boldsymbol{x}$ since

$k_{\text{odd}}(-\boldsymbol{x},\boldsymbol{z})=\frac{1}{2}(k(-\boldsymbol{x},\boldsymbol{z})-k(\boldsymbol{x},\boldsymbol{z}))=-k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})$
(53)

This means that the kernel $k_{\text{odd}}$ is an odd reproducing kernel with
an associated RKHS [6]. Any function
$f(\boldsymbol{x})=\sum_{i=1}^{N}a_{i}k_{\text{odd}}(\boldsymbol{x},\boldsymbol{x}_{i})$
in the RKHS defined by $k_{\text{odd}}$ will then be odd, since
$f(-\boldsymbol{x})=-f(\boldsymbol{x})$.

An odd kernel will not be shift-invariant, and it will not have a signature.
Instead, it is given by a sum of two signatures as

$k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})=\frac{1}{2}(g(\boldsymbol{x}-\boldsymbol{z})-g(\boldsymbol{x}+\boldsymbol{z}))$
(54)

###  4.5 Odd curl-free and odd symplectic kernel

The odd curl-free kernel was proposed in [15]. It is derived from the odd
Gaussian kernel

$\displaystyle k_{\sigma,o}(\boldsymbol{x},\boldsymbol{z})$
$\displaystyle=\frac{1}{2}\left(e^{-\frac{\|\boldsymbol{x}-\boldsymbol{z}\|^{2}}{2\sigma^{2}}}-e^{-\frac{\|\boldsymbol{x}+\boldsymbol{z}\|^{2}}{2\sigma^{2}}}\right)$
(55)

The odd curl-free kernel is then derived from

$\boldsymbol{K}_{c,o}(\boldsymbol{x},\boldsymbol{z})=-\nabla\nabla^{{\text{T}}}k_{\sigma,o}(\boldsymbol{x},\boldsymbol{z})$
(56)

which results in

$\boldsymbol{K}_{c,o}(\boldsymbol{x},\boldsymbol{z})=\frac{1}{2\sigma^{2}}\biggl{(}e^{-\frac{\|\boldsymbol{x}-\boldsymbol{z}\|^{2}}{2\sigma^{2}}}\biggl{(}{\boldsymbol{I}}-\frac{(\boldsymbol{x}-\boldsymbol{z})^{{\text{T}}}(\boldsymbol{x}-\boldsymbol{z})}{\sigma^{2}}\biggr{)}\\\
-e^{-\frac{\|\boldsymbol{x}+\boldsymbol{z}\|^{2}}{2\sigma^{2}}}\biggl{(}{\boldsymbol{I}}-\frac{(\boldsymbol{x}+\boldsymbol{z})^{{\text{T}}}(\boldsymbol{x}+\boldsymbol{z})}{\sigma^{2}}\biggr{)}\biggr{)}$
(57)

which is the same as using (45) to write

$\boldsymbol{K}_{c,o}(\boldsymbol{x},\boldsymbol{z})=\frac{1}{2}\left(\boldsymbol{G}_{c}(\boldsymbol{x}-\boldsymbol{z})-\boldsymbol{G}_{c}(\boldsymbol{x}+\boldsymbol{z})\right)$
(58)

Finally, the odd symplectic kernel is written combining (46) and (58) as

$\boldsymbol{K}_{s,o}(\boldsymbol{x},\boldsymbol{z})=\frac{1}{2}\left(\boldsymbol{G}_{s}(\boldsymbol{x}-\boldsymbol{z})-\boldsymbol{G}_{s}(\boldsymbol{x}+\boldsymbol{z})\right)$
(59)

where ${\boldsymbol{G}_{s}}$ is defined as in (51).

##  5 Random feature approximation

###  5.1 Random features approximation of the Gaussian separable kernel

Let ${{\mathcal{F}}[\cdot]}$ and ${{{\mathcal{F}}}^{-1}[\cdot]}$ denote the
Fourier transformation and inverse Fourier transformation, respectively.
Defining $p_{\sigma}(\boldsymbol{w})$ as the inverse Fourier transformation of
the Gaussian kernel function ${g_{\sigma}(\boldsymbol{x})}$ results in

$\displaystyle g_{\sigma}(\boldsymbol{x})$
$\displaystyle={\mathcal{F}}[p_{\sigma}(\boldsymbol{w})]=\int_{-\infty}^{\infty}e^{-i\boldsymbol{w}^{{\text{T}}}\boldsymbol{x}}p_{\sigma}(\boldsymbol{w})d\boldsymbol{w}$
(60) $\displaystyle p_{\sigma}(\boldsymbol{w})$
$\displaystyle={{\mathcal{F}}}^{-1}[g_{\sigma}(\boldsymbol{x})]=\left(\frac{1}{2\pi}\right)^{n}\int_{-\infty}^{\infty}e^{i\boldsymbol{w}^{{\text{T}}}\boldsymbol{x}}g_{\sigma}(\boldsymbol{x})d\boldsymbol{x}$
(61)

which gives

$p_{\sigma}(\boldsymbol{w})=\left(\frac{1}{\sqrt{2\pi\rho^{2}}}\right)^{n}e^{-\frac{\|\boldsymbol{w}\|^{2}}{2\rho^{2}}}\sim{\mathcal{N}}\left({\boldsymbol{0}},\rho^{2}{\boldsymbol{I}}\right)$
(62)

where ${\rho^{2}=2/\sigma^{2}}$ [24]. As ${g_{\sigma}(\boldsymbol{x})}$ is a
real-valued function, it follows that

$g_{\sigma}(\boldsymbol{x})=\int_{-\infty}^{\infty}\cos(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})p_{\sigma}(\boldsymbol{w})d\boldsymbol{w}$
(63)

Since ${\cos(\alpha-\beta)=\cos\alpha\cos\beta+\sin\alpha\sin\beta}$, this can
be written

$g_{\sigma}(\boldsymbol{x}-\boldsymbol{z})=\int_{-\infty}^{\infty}\left(\cos(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\cos(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})+\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})\right)p_{\sigma}(\boldsymbol{w})d\boldsymbol{w}$

It follows that $g_{\sigma}(\boldsymbol{x}-\boldsymbol{z})$ is the expected
value

$g_{\sigma}(\boldsymbol{x}-\boldsymbol{z})={\mathbb{E}}_{w}\begin{bmatrix}\cos(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\cos(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})+\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})\end{bmatrix}$
(64)

where $\boldsymbol{w}$ has the distribution $p_{\sigma}(\boldsymbol{w})$. The
Gaussian kernel can then be approximated with the finite sum

$g_{\sigma}(\boldsymbol{x}-\boldsymbol{z})\approx\frac{1}{d}\sum_{i=1}^{d}\left(\cos(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{x})\cos(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{z})+\sin(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{z})\right)$
(65)

where the weights ${\\{\boldsymbol{w}_{i}\\}_{i=1}^{d}\in{\mathbb{R}}^{n}}$
are drawn i.i.d. with distribution ${p_{\sigma}(\boldsymbol{w})}$. Then the
scalar Gaussian kernel can be approximated as

$k_{\sigma}(\boldsymbol{x},\boldsymbol{z})\approx\boldsymbol{\psi}_{\sigma}(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\psi}_{\sigma}(\boldsymbol{z})$
(66)

where the feature map for the scalar Gaussian kernel is

$\boldsymbol{\psi}_{\sigma}(\boldsymbol{x})=\frac{1}{\sqrt{d}}\begin{bmatrix}\cos(\boldsymbol{w}_{1}^{{\text{T}}}\boldsymbol{x})\\\
\vdots\\\ \cos(\boldsymbol{w}_{d}^{{\text{T}}}\boldsymbol{x})\\\
\sin(\boldsymbol{w}_{1}^{{\text{T}}}\boldsymbol{x})\\\ \vdots\\\
\sin(\boldsymbol{w}_{d}^{{\text{T}}}\boldsymbol{x})\end{bmatrix}\in{\mathbb{R}}^{2d}$
(67)

The matrix-valued Gaussian separable kernel function can then be approximated
by the matrix-valued feature map as follows [2]

$\boldsymbol{K}_{\sigma}(\boldsymbol{x},\boldsymbol{z})\approx\boldsymbol{\psi}_{\sigma}(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\psi}_{\sigma}(\boldsymbol{z}){\boldsymbol{I}}_{n}=\boldsymbol{\Psi}_{\sigma}(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\Psi}_{\sigma}(\boldsymbol{z})$
(68)

where

$\boldsymbol{\Psi}_{\sigma}(\boldsymbol{x})=\boldsymbol{\psi}_{\sigma}(\boldsymbol{x})\otimes{\boldsymbol{I}}_{n}\in{\mathbb{R}}^{2dn\times
n}$ (69)

and ${\otimes}$ is the Kronecker product.

###  5.2 Random features approximation of the curl-free kernel

The signature
$\boldsymbol{G}_{c}(\boldsymbol{x})=-\nabla\nabla^{{\text{T}}}g_{\sigma}(\boldsymbol{x})$
of the curl-free kernel is found by differentiation of (63) to be given by

$\boldsymbol{G}_{c}(\boldsymbol{x})=\int_{-\infty}^{\infty}\cos(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\boldsymbol{w}\boldsymbol{w}^{{\text{T}}}p_{\sigma}(\boldsymbol{w})d\boldsymbol{w}$
(70)

It follows that

$\boldsymbol{G}_{c}(\boldsymbol{x}-\boldsymbol{z})={\mathbb{E}}_{w}\begin{bmatrix}\left(\cos(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\cos(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})+\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})\right)\boldsymbol{w}\boldsymbol{w}^{{\text{T}}}\end{bmatrix}$
(71)

The curl-free kernel can then be approximated using the finite sum

$\boldsymbol{G}_{c}(\boldsymbol{x}-\boldsymbol{z})\approx\frac{1}{d}\sum_{i=1}^{d}\left(\cos(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{x})\cos(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{z})+\sin(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{z})\right)\boldsymbol{w}_{i}\boldsymbol{w}_{i}^{{\text{T}}}$
(72)

where the weights ${\\{\boldsymbol{w}_{i}\\}_{i=1}^{d}\in{\mathbb{R}}^{n}}$
are drawn i.i.d. with distribution ${p_{\sigma}(\boldsymbol{w})}$. The matrix-
valued curl-free kernel function can then be approximated by the matrix-valued
feature map as follows

$\boldsymbol{K}_{c}(\boldsymbol{x},\boldsymbol{z})\approx\boldsymbol{\Psi}_{c}(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\Psi}_{c}(\boldsymbol{z})$
(73)

where

$\boldsymbol{\Psi}_{c}(\boldsymbol{x})=\frac{1}{\sqrt{d}}\begin{bmatrix}\cos(\boldsymbol{w}_{1}^{{\text{T}}}\boldsymbol{x})\boldsymbol{w}_{1}^{{\text{T}}}\\\
\vdots\\\
\cos(\boldsymbol{w}_{d}^{{\text{T}}}\boldsymbol{x})\boldsymbol{w}_{d}^{{\text{T}}}\\\
\sin(\boldsymbol{w}_{1}^{{\text{T}}}\boldsymbol{x})\boldsymbol{w}_{1}^{{\text{T}}}\\\
\vdots\\\
\sin(\boldsymbol{w}_{d}^{{\text{T}}}\boldsymbol{x})\boldsymbol{w}_{d}^{{\text{T}}}\end{bmatrix}\in{\mathbb{R}}^{2d\times
n}$ (74)

###  5.3 Random features approximation of the symplectic kernel

The random features approximation for the symplectic kernel
${\boldsymbol{K}_{s}(\boldsymbol{x},\boldsymbol{z})=\boldsymbol{J}\boldsymbol{K}_{c}(\boldsymbol{x},\boldsymbol{z})\boldsymbol{J}^{{\text{T}}}}$
is found from (70) to be

$\displaystyle\boldsymbol{G}_{s}(\boldsymbol{x})$
$\displaystyle=\boldsymbol{J}\boldsymbol{G}_{c}(\boldsymbol{x})\boldsymbol{J}^{{\text{T}}}=\boldsymbol{J}\int_{-\infty}^{\infty}\cos(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\boldsymbol{w}\boldsymbol{w}^{{\text{T}}}p_{\sigma}(\boldsymbol{w})d\boldsymbol{w}\boldsymbol{J}^{{\text{T}}}$
(75)

It follows that

$\boldsymbol{G}_{s}(\boldsymbol{x}-\boldsymbol{z})={\mathbb{E}}_{w}\begin{bmatrix}\boldsymbol{J}\boldsymbol{w}\left(\cos(\boldsymbol{w}^{{\text{T}}}(\boldsymbol{x}-\boldsymbol{z}))\right)(\boldsymbol{J}\boldsymbol{w})^{{\text{T}}}\end{bmatrix}$
(76)

The symplectic kernel can then be approximated using the finite sum

$\boldsymbol{G}_{s}(\boldsymbol{x}-\boldsymbol{z})\approx\frac{1}{d}\sum_{i=1}^{d}\boldsymbol{J}\boldsymbol{w}_{i}\left(\cos(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{x})\cos(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{z})+\sin(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{z})\right)(\boldsymbol{J}\boldsymbol{w}_{i})^{{\text{T}}}$

where the weights ${\\{\boldsymbol{w}_{i}\\}_{i=1}^{d}\in{\mathbb{R}}^{n}}$
are drawn i.i.d. with distribution ${p_{\sigma}(\boldsymbol{w})}$. The matrix-
valued symplectic kernel function can then be approximated by the matrix-
valued feature map as follows

$\boldsymbol{K}_{s}(\boldsymbol{x},\boldsymbol{z})\approx\boldsymbol{\Psi}_{s}(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\Psi}_{s}(\boldsymbol{z})$
(77)

where

$\boldsymbol{\Psi}_{s}(\boldsymbol{x})=\frac{1}{\sqrt{d}}\begin{bmatrix}\cos(\boldsymbol{w}_{1}^{{\text{T}}}\boldsymbol{x})(\boldsymbol{J}\boldsymbol{w}_{1})^{{\text{T}}}\\\
\vdots\\\
\cos(\boldsymbol{w}_{d}^{{\text{T}}}\boldsymbol{x})(\boldsymbol{J}\boldsymbol{w}_{d})^{{\text{T}}}\\\
\sin(\boldsymbol{w}_{1}^{{\text{T}}}\boldsymbol{x})(\boldsymbol{J}\boldsymbol{w}_{1})^{{\text{T}}}\\\
\vdots\\\
\sin(\boldsymbol{w}_{d}^{{\text{T}}}\boldsymbol{x})(\boldsymbol{J}\boldsymbol{w}_{d})^{{\text{T}}}\end{bmatrix}\in{\mathbb{R}}^{2d\times
n}$ (78)

###  5.4 Random features approximation of the odd Gaussian kernel

The odd Gaussian kernel ${k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})}$ is
derived from the Gaussian kernel
$k_{\sigma}(\boldsymbol{x},\boldsymbol{z})={g_{\sigma}(\boldsymbol{x}-\boldsymbol{z})}$
where [6]

$\displaystyle k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})$
$\displaystyle=\frac{1}{2}\left(g_{\sigma}(\boldsymbol{x}-\boldsymbol{z})-g_{\sigma}(\boldsymbol{x}+\boldsymbol{z})\right)$
(79)

The odd Gaussian kernel is not shift-invariant and does not have a signature.
This means that Bochner’s theorem [24] cannot be used directly on the odd
Gaussian kernel. Instead, it is used that the odd Gaussian kernel is the sum
of the two signature functions $g_{\sigma}(\boldsymbol{x}-\boldsymbol{z})$ and
$g_{\sigma}(\boldsymbol{x}+\boldsymbol{z})$, and Bochner’s theorem can be
applied for each of these functions. This gives

$\displaystyle g_{\sigma}(\boldsymbol{x}-\boldsymbol{z})$
$\displaystyle=\int_{-\infty}^{\infty}\cos(\boldsymbol{w}^{{\text{T}}}(\boldsymbol{x}-\boldsymbol{z}))p_{\sigma}(\boldsymbol{w})d\boldsymbol{w}$
(80) $\displaystyle g_{\sigma}(\boldsymbol{x}+\boldsymbol{z})$
$\displaystyle=\int_{-\infty}^{\infty}\cos(\boldsymbol{w}^{{\text{T}}}(\boldsymbol{x}+\boldsymbol{z}))p_{\sigma}(\boldsymbol{w})d\boldsymbol{w}$
(81)

Using ${\cos(\alpha\pm\beta)=\cos\alpha\cos\beta\mp\sin\alpha\sin\beta}$, the
odd kernel can be written as

$\displaystyle k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})$
$\displaystyle=\int_{-\infty}^{\infty}\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})p_{\sigma}(\boldsymbol{w})d\boldsymbol{w}$
(82)

It follows that the odd kernel is the expected value

$k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})={\mathbb{E}}_{w}\begin{bmatrix}\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})\end{bmatrix}$
(83)

where $\boldsymbol{w}$ has distribution $p_{\sigma}(\boldsymbol{w})$. The odd
kernel can then be approximated using the finite sum

$k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})\approx\frac{1}{d}\sum_{i=1}^{d}\sin(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{z})$
(84)

where the weights ${\\{\boldsymbol{w}_{i}\\}_{i=1}^{d}\in{\mathbb{R}}^{n}}$
are drawn i.i.d. with distribution ${p_{\sigma}(\boldsymbol{w})}$. The odd
Gaussian kernel can then be approximated by

$k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})\approx\boldsymbol{\psi}_{\text{odd}}(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\psi}_{\text{odd}}(\boldsymbol{z})$
(85)

where the feature map is

$\boldsymbol{\psi}_{\text{odd}}(\boldsymbol{x})=\frac{1}{\sqrt{d}}\begin{bmatrix}\sin(\boldsymbol{w}_{1}^{{\text{T}}}\boldsymbol{x})\\\
\vdots\\\
\sin(\boldsymbol{w}_{d}^{{\text{T}}}\boldsymbol{x})\end{bmatrix}\in{\mathbb{R}}^{d}$
(86)

It is noted that the approximation for the odd kernel can also be found from

$k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})=\frac{1}{2}\left(k_{\sigma}(\boldsymbol{x},\boldsymbol{z})-k_{\sigma}(-\boldsymbol{x},\boldsymbol{z})\right)$
(87)

where $k_{\sigma}(\boldsymbol{x},\boldsymbol{z})$ is the approximation of the
Gaussian kernel given in (67). The same method can be used to derive the
random features for the even kernel

$k_{\text{even}}(\boldsymbol{x},\boldsymbol{z})=\frac{1}{2}\left(k(\boldsymbol{x},\boldsymbol{z})+k(-\boldsymbol{x},\boldsymbol{z})\right)$
(88)

in [6].

###  5.5 Random features approximation of the odd curl-free kernel

Using the curl-free kernel
${\boldsymbol{K}_{c}(\boldsymbol{x},\boldsymbol{z})=\boldsymbol{G}_{c}(\boldsymbol{x}-\boldsymbol{z})}$
the random features approximation for the odd curl-free kernel is derived in
the following

$\displaystyle\boldsymbol{K}_{c,o}(\boldsymbol{x},\boldsymbol{z})$
$\displaystyle=\frac{1}{2}\left(\boldsymbol{G}_{c}(\boldsymbol{x}-\boldsymbol{z})-\boldsymbol{G}_{c}(\boldsymbol{x}+\boldsymbol{z})\right)$
$\displaystyle=\frac{1}{2}\int_{-\infty}^{\infty}\boldsymbol{w}\boldsymbol{w}^{{\text{T}}}\cos(\boldsymbol{w}^{{\text{T}}}(\boldsymbol{x}-\boldsymbol{z}))p_{\sigma}(\boldsymbol{w})d\boldsymbol{w}$
$\displaystyle\quad-\frac{1}{2}\int_{-\infty}^{\infty}\boldsymbol{w}\boldsymbol{w}^{{\text{T}}}\cos(\boldsymbol{w}^{{\text{T}}}(\boldsymbol{x}+\boldsymbol{z}))p_{\sigma}(\boldsymbol{w})d\boldsymbol{w}$
$\displaystyle=\int_{-\infty}^{\infty}\boldsymbol{w}\boldsymbol{w}^{{\text{T}}}\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})p_{\sigma}(\boldsymbol{w})d\boldsymbol{w}$
(89)

where the trigonometric identity
${\cos(\alpha\pm\beta)=\cos\alpha\cos\beta\mp\sin\alpha\sin\beta}$ is applied.
It follows that

$\boldsymbol{K}_{c,o}(\boldsymbol{x},\boldsymbol{z})={\mathbb{E}}_{w}\begin{bmatrix}\boldsymbol{w}\boldsymbol{w}^{{\text{T}}}\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})\end{bmatrix}$
(90)

The odd curl-free kernel can then be approximated using the finite sum

$\boldsymbol{K}_{c,o}(\boldsymbol{x},\boldsymbol{z})\approx\frac{1}{d}\sum_{i=1}^{d}\boldsymbol{w}_{i}\boldsymbol{w}_{i}^{{\text{T}}}\sin(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{z})$
(91)

where the weights ${\\{\boldsymbol{w}_{i}\\}_{i=1}^{d}\in{\mathbb{R}}^{n}}$
are drawn i.i.d. with distribution ${p_{\sigma}(\boldsymbol{w})}$. The matrix-
valued odd curl-free kernel function can then be approximated by the matrix-
valued feature map as follows

$\boldsymbol{K}_{c,o}(\boldsymbol{x},\boldsymbol{z})\approx\boldsymbol{\Psi}_{c,o}(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\Psi}_{c,o}(\boldsymbol{z})$
(92)

where

$\boldsymbol{\Psi}_{c,o}(\boldsymbol{x})=\frac{1}{\sqrt{d}}\begin{bmatrix}\sin(\boldsymbol{w}_{1}^{{\text{T}}}\boldsymbol{x})\boldsymbol{w}_{1}^{{\text{T}}}\\\
\vdots\\\
\sin(\boldsymbol{w}_{d}^{{\text{T}}}\boldsymbol{x})\boldsymbol{w}_{d}^{{\text{T}}}\end{bmatrix}\in{\mathbb{R}}^{d\times
n}$ (93)

###  5.6 Random features approximation of the odd symplectic kernel

The random features approximation for the odd symplectic kernel is found by
applying Bochner’s theorem to the signature functions in

$\boldsymbol{K}_{s,o}(\boldsymbol{x},\boldsymbol{z})=\frac{1}{2}\left(\boldsymbol{G}_{s}(\boldsymbol{x}-\boldsymbol{z})-\boldsymbol{G}_{s}(\boldsymbol{x}+\boldsymbol{z})\right)$
(94)

which gives

$\displaystyle\boldsymbol{K}_{s,o}(\boldsymbol{x},\boldsymbol{z})$
$\displaystyle=\frac{1}{2}\boldsymbol{J}\left(\boldsymbol{G}_{c}(\boldsymbol{x}-\boldsymbol{z})-\boldsymbol{G}_{c}(\boldsymbol{x}+\boldsymbol{z})\right)\boldsymbol{J}^{{\text{T}}}$
(95)
$\displaystyle=\boldsymbol{J}\left(\int_{-\infty}^{\infty}\boldsymbol{w}\boldsymbol{w}^{{\text{T}}}\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})p_{\sigma}(\boldsymbol{w})d\boldsymbol{w}\right)\boldsymbol{J}^{{\text{T}}}$
(96)

It follows that

$\boldsymbol{K}_{s,o}(\boldsymbol{x},\boldsymbol{z})={\mathbb{E}}_{w}\begin{bmatrix}\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})(\boldsymbol{J}\boldsymbol{w})\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})(\boldsymbol{J}\boldsymbol{w})^{{\text{T}}}\end{bmatrix}$
(97)

The odd symplectic kernel can then be approximated using the finite sum

$\boldsymbol{K}_{s,o}(\boldsymbol{x},\boldsymbol{z})\approx\frac{1}{d}\sum_{i=1}^{d}\sin(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{x})(\boldsymbol{J}\boldsymbol{w}_{i})\sin(\boldsymbol{w}_{i}^{{\text{T}}}\boldsymbol{z})(\boldsymbol{J}\boldsymbol{w}_{i})^{{\text{T}}}$
(98)

where the weights ${\\{\boldsymbol{w}_{i}\\}_{i=1}^{d}\in{\mathbb{R}}^{n}}$
are drawn i.i.d. with distribution ${p_{\sigma}(\boldsymbol{w})}$. The matrix-
valued odd symplectic kernel function can then be approximated by the matrix-
valued feature map as follows

$\boldsymbol{K}_{s,o}(\boldsymbol{x},\boldsymbol{z})\approx\boldsymbol{\Psi}_{s,o}(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\Psi}_{s,o}(\boldsymbol{z})$
(99)

where

$\boldsymbol{\Psi}_{s,o}(\boldsymbol{x})=\frac{1}{\sqrt{d}}\begin{bmatrix}\sin(\boldsymbol{w}_{1}^{{\text{T}}}\boldsymbol{x})(\boldsymbol{J}\boldsymbol{w}_{1})^{{\text{T}}}\\\
\vdots\\\
\sin(\boldsymbol{w}_{d}^{{\text{T}}}\boldsymbol{x})(\boldsymbol{J}\boldsymbol{w}_{d})^{{\text{T}}}\end{bmatrix}\in{\mathbb{R}}^{d\times
n}$ (100)

###  5.7 Learned Hamiltonian using RFF

The Hamiltonian dynamics are given by

$\boldsymbol{f}(\boldsymbol{x})=\boldsymbol{J}\nabla H(\boldsymbol{x})$ (101)

and using the symplectic kernel the learned Hamiltonian dynamics are given by

$\boldsymbol{f}(\boldsymbol{x})=-\boldsymbol{J}\nabla\sum_{i=1}^{N}\nabla^{{\text{T}}}g(\boldsymbol{x}-\boldsymbol{x}_{i})\boldsymbol{c}_{i}$
(102)

where ${\boldsymbol{c}_{i}=\boldsymbol{J}^{{\text{T}}}\boldsymbol{a}_{i}}$.
The learned Hamiltonian is then

$\hat{H}(\boldsymbol{x})=-\sum_{i=1}^{N}\nabla^{{\text{T}}}g(\boldsymbol{x}-\boldsymbol{x}_{i})\boldsymbol{c}_{i}$
(103)

Starting with the real-valued Fourier-transform of the odd Gaussian kernel

$k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})=\frac{1}{2}\int_{-\infty}^{\infty}\left(\cos(\boldsymbol{w}^{{\text{T}}}(\boldsymbol{x}-\boldsymbol{z}))-\cos(\boldsymbol{w}^{{\text{T}}}(\boldsymbol{x}+\boldsymbol{z}))\right)p(\boldsymbol{w})d\boldsymbol{w}$
(104)

taking the gradient with respect to ${\boldsymbol{x}}$ gives

$\nabla^{{\text{T}}}k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})=-\frac{1}{2}\int_{-\infty}^{\infty}\boldsymbol{w}^{{\text{T}}}\left(\sin(\boldsymbol{w}^{{\text{T}}}(\boldsymbol{x}-\boldsymbol{z}))-\sin(\boldsymbol{w}^{{\text{T}}}(\boldsymbol{x}+\boldsymbol{z}))\right)p(\boldsymbol{w})d\boldsymbol{w}$
(105)

The trigonometric identities
${\sin(\alpha\pm\beta)=\sin\alpha\cos\beta\pm\cos\alpha\sin\beta}$ gives

$\nabla^{{\text{T}}}k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})=\int_{-\infty}^{\infty}\boldsymbol{w}^{{\text{T}}}\cos(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})p(\boldsymbol{w})d\boldsymbol{w}$
(106)

which means that

$\nabla^{{\text{T}}}k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})={\mathbb{E}}_{w}\begin{bmatrix}\boldsymbol{w}^{{\text{T}}}\cos(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{z})\end{bmatrix}$
(107)

It is seen that the first derivative of the odd kernel can be approximated
using the finite sum

$\nabla^{{\text{T}}}k_{\text{odd}}(\boldsymbol{x},\boldsymbol{z})\approx\frac{1}{d}\sum_{j=1}^{d}\boldsymbol{w}_{j}^{{\text{T}}}\cos(\boldsymbol{w}_{j}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}_{j}^{{\text{T}}}\boldsymbol{z})$
(108)

where the weights ${\\{\boldsymbol{w}_{j}\\}_{j=1}^{d}\in{\mathbb{R}}^{n}}$
are i.i.d. with distribution ${p(\boldsymbol{w})}$. The expression for the
learned Hamiltonian using RFF is then written as

$\displaystyle\hat{H}(\boldsymbol{x})$
$\displaystyle=-\sum_{i=1}^{N}\frac{1}{d}\sum_{j=1}^{d}\boldsymbol{w}_{j}^{{\text{T}}}\cos(\boldsymbol{w}_{j}^{{\text{T}}}\boldsymbol{x})\sin(\boldsymbol{w}_{j}^{{\text{T}}}\boldsymbol{x}_{i})\boldsymbol{c}_{i}$
(109)
$\displaystyle=-\frac{1}{d}\sum_{j=1}^{d}\cos(\boldsymbol{w}_{j}^{{\text{T}}}\boldsymbol{x})\sum_{i=1}^{N}\boldsymbol{w}_{j}^{{\text{T}}}\sin(\boldsymbol{w}_{j}^{{\text{T}}}\boldsymbol{x}_{i})\boldsymbol{c}_{i}$
(110)
$\displaystyle=-\frac{1}{\sqrt{d}}\sum_{j=1}^{d}\cos(\boldsymbol{w}_{j}^{{\text{T}}}\boldsymbol{x})\alpha_{j}$
(111)

Note that the dataset is included in the model weights as follows

$\alpha_{j}=\sum_{i=1}^{N}\psi(\boldsymbol{x}_{i})\boldsymbol{c}_{i}=\frac{1}{\sqrt{d}}\sum_{i=1}^{N}\boldsymbol{w}_{j}^{{\text{T}}}\sin(\boldsymbol{w}_{j}^{{\text{T}}}\boldsymbol{x}_{i})\boldsymbol{c}_{i}$
(112)

where ${\boldsymbol{c}_{i}=\boldsymbol{J}^{{\text{T}}}\boldsymbol{a}_{i}}$ and

$\psi(\boldsymbol{x})=\frac{1}{\sqrt{d}}\sin(\boldsymbol{w}^{{\text{T}}}\boldsymbol{x})(\boldsymbol{J}\boldsymbol{w})^{{\text{T}}}$
(113)

is the scalar component of the RFF approximation of the odd symplectic kernel
in (99). It is seen that the model weight vector
${\boldsymbol{\alpha}\in{\mathbb{R}}^{d}}$ is made up of ${d}$ scalar
components
${\boldsymbol{\alpha}=\begin{bmatrix}\alpha_{1},\dots,\alpha_{d}\end{bmatrix}^{{\text{T}}}}$,
meaning the expression in (111) is scalar. Finally, the expression for the
learned Hamiltonian is written as

$\hat{H}(\boldsymbol{x})=\Gamma(\boldsymbol{x})^{{\text{T}}}\boldsymbol{\alpha}$
(114)

where

$\Gamma(\boldsymbol{x})=\frac{1}{\sqrt{d}}\begin{bmatrix}\cos(\boldsymbol{w}_{1}^{{\text{T}}}\boldsymbol{x})\\\
\vdots\\\
\cos(\boldsymbol{w}_{d}^{{\text{T}}}\boldsymbol{x})\end{bmatrix}\in{\mathbb{R}}^{d}$
(115)

##  6 Experiments

The proposed kernel was evaluated by learning the Hamiltonian dynamics of
three Hamiltonian systems with odd vector fields. The regularized least-
squares problem in (26) was solved using (27), and the resulting learned
vector field was given by (28). The RFF approximation of the Gaussian
separable kernel in (68) and the RFF approximation of the odd symplectic
kernel in (99) were used and compared for the three systems.

###  6.1 Hyperparameter tuning

The hyperparameters ${\sigma}$ and ${\lambda}$ greatly influence the learned
vector fields. The hyperparameters were tuned using a genetic algorithm [25]
in MATLAB, by minimizing the cross-validation error [26] over the training
set. The training set
${{\mathcal{Z}}=\\{(\boldsymbol{x}_{i},\boldsymbol{y}_{i})\in{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\\}_{i=1}^{N}}$
was split into mutually exclusive subsets
${{\mathcal{Z}}_{1},\dots\mathcal{Z}_{k}}$, and for each iteration
${i\in\\{1,\dots,k\\}}$, the model was trained on the subset
${\hat{{\mathcal{Z}}}_{i}={\mathcal{Z}}\setminus{\mathcal{Z}}_{i}}$ and
evaluated on ${{\mathcal{Z}}_{i}}$. Formally, the hyperparameter optimization
is written as [6]

$\min_{\sigma,\lambda}\frac{1}{k}\sum_{i=1}^{k}\text{MSE}\left(\boldsymbol{f}_{\hat{{\mathcal{Z}}}_{i}},{\mathcal{Z}}_{i}\right)$
(116)

where ${\boldsymbol{f}_{\hat{{\mathcal{Z}}}_{i}}}$ is the learned vector field
trained on
${\hat{{\mathcal{Z}}}_{i}={\mathcal{Z}}\setminus{\mathcal{Z}}_{i}}$, using the
hyperparameters ${\sigma}$ and ${\lambda}$, and MSE is the empirical mean
square error between the learned model and ${{\mathcal{Z}}_{i}}$.

###  6.2 Simple pendulum

A simple pendulum is modeled with a point mass ${m}$ and at the end of a mass-
less rod of length ${l}$. The pendulum angle is ${\theta}$. The equation of
motion is given by

$\ddot{\theta}=-\frac{g}{l}\sin(\theta)$ (117)

where ${g}$ is the acceleration of gravity. The generalized coordinate is
${q=\theta}$, the kinetic energy is ${T=\frac{1}{2}ml^{2}{\dot{q}}^{2}}$ and
the potential energy is ${U=mgl(1-\cos(q))}$. The Lagrangian is defined as
${L=T-U}$ and from (31) the generalized momentum is

$p=\frac{\partial L}{\partial{\dot{q}}}=ml^{2}{\dot{q}}$ (118)

The Hamiltonian is

$H(q,p)=p{\dot{q}}-L=\frac{p^{2}}{2ml^{2}}+mgl(1-\cos(q))$ (119)

The Hamiltonian dynamics are then given by

${\dot{q}}=\dfrac{\partial H}{\partial
p}=\frac{p}{ml^{2}},\qquad{\dot{p}}=-\dfrac{\partial H}{\partial
q}=-mgl\sin(q)$ (120)

Figure 1(a) shows the true system with parameters ${m=1}$, ${l=1}$, and
${g=9.81}$. Three trajectories were generated, and the system was simulated
with three different initial conditions:
${\boldsymbol{x}_{0}=\\{\begin{bmatrix}\frac{2\pi}{5},0\end{bmatrix}^{{\text{T}}},\begin{bmatrix}\frac{4\pi}{5},0\end{bmatrix}^{{\text{T}}},\begin{bmatrix}\frac{19\pi}{20},-4\end{bmatrix}^{{\text{T}}}\\}}$.
The time step was set to ${h=0.1}$ as the system was simulated for
${t\in\begin{bmatrix}0,0.7\end{bmatrix}}$ seconds, giving ${8}$ data points
for each trajectory, and ${N=24}$ total data points. The velocities
${\boldsymbol{y}}$ were sampled at each trajectory point, and zero mean
Gaussian noise with standard deviation ${\sigma_{n}=0.01}$ was added to the
trajectory and velocity data. Figure 1(b) shows the resulting data set.

The ${d=50}$ random features were used for the Gaussian model, and ${d=400}$
random features were used for the odd symplectic model, giving an equal number
of model coefficients ${\boldsymbol{\alpha}}$ for each model. For additional
comparison, the symplectic kernel (77) proposed in [19] was also used to learn
the dynamics of the simple pendulum, and ${d=200}$ random features are used.

(a) True system

(b) Data set

(c) Gaussian model

(d) Symplectic model

(e) Odd sympl. model

Figure 1: Stream and trajectory plots for the simple pendulum and extracted
data set, and the resulting learned models using the separable Gaussian
kernel, symplectic kernel, and the odd symplectic kernel.

(a) Test trajectory

(b) Test trajectory error

Figure 2: Comparison of the two learned models against the simple pendulum
system, using the test trajectory.

Figures 1(c) and 1(e) show phase plots of the learned models using the
separable Gaussian kernel and the odd symplectic kernel, respectively. The
function learned with the Gaussian separable kernel did not give an accurate
representation of the true dynamics from such a limited dataset. The model
learned with the odd symplectic kernel was accurate and gave a good
representation of the vector field of the simple pendulum system. It is seen
from Figure 1(e) that symmetry and energy conservation lead to periodic orbits
like the true system. The symplectic kernel also failed to learn the dynamics
of the simple pendulum accurately (Figure 1(d)). The learned vector field is
similar to the learned Gaussian model, but the symplectic kernel enforced
energy conservation evident from the periodic orbits.

A separate test trajectory was simulated to test the generalized performance
of the learned models. The initial condition was
${\boldsymbol{x}_{0}=\begin{bmatrix}\frac{\pi}{2},0\end{bmatrix}^{{\text{T}}}}$
and the time horizon is ${t\in\begin{bmatrix}0,2\end{bmatrix}}$ seconds. The
error between the true system and the learned model trajectories was defined
as ${\text{Err}=\|\boldsymbol{x}_{gt}-\boldsymbol{x}_{l}\|}$. Figure 2(a)
shows the four resulting trajectories, and Figure 2(b) shows the error for
each time step. The results showed that the odd symplectic model was far more
accurate than both the Gaussian separable model and the symplectic model,
which both failed to generalize beyond the area close to the data set.

###  6.3 Cart-pole

The Cart-Pole [27] is a planar, underactuated mechanical system where the task
is to balance an inverted pendulum starting at an arbitrary initial condition,
using only the linear motion of the cart as the input. For this task of system
identification, the un-actuated system is modeled. The system consists of a
cart moving linearly in the horizontal ${x}$-direction with mass ${m_{c}}$ and
an inverted pendulum with point mass ${m_{p}}$ and massless rod with length
${l}$, connected to the cart through a pivot. The angle between the pendulum
and the vertical axis is denoted by ${\theta}$, which is zero at the upright
position. The kinematics of the system are given by

$\boldsymbol{r}_{c}=\begin{bmatrix}x\\\
0\end{bmatrix},\quad\boldsymbol{r}_{p}=\begin{bmatrix}x+l\sin(\theta)\\\
l\cos(\theta)\end{bmatrix}$ (121)

where ${\boldsymbol{r}_{c}}$ and ${\boldsymbol{r}_{p}}$ are the positions in
the ${xy}$-plane of the cart ${m_{c}}$ and pendulum ${m_{p}}$, respectively.
The velocities are

$\boldsymbol{v}_{c}=\begin{bmatrix}{\dot{x}}\\\
0\end{bmatrix},\quad\boldsymbol{v}_{p}=\begin{bmatrix}{\dot{x}}+l\dot{\theta}\cos(\theta)\\\
-l\dot{\theta}\sin(\theta)\end{bmatrix}$ (122)

The kinetic energy ${T}$ and potential energy ${U}$ of the system are

$\displaystyle T$
$\displaystyle=\frac{1}{2}\left(m_{c}+m_{p}\right){\dot{x}}^{2}+m_{p}l{\dot{x}}\dot{\theta}\cos(\theta)+\frac{1}{2}m_{p}l^{2}\dot{\theta}^{2}$
(123) $\displaystyle U$ $\displaystyle=m_{p}gl\cos(\theta)$ (124)

Defining the generalized coordinate
${\boldsymbol{q}=\begin{bmatrix}x,\theta\end{bmatrix}^{{\text{T}}}}$ and its
time derivative
${\dot{\boldsymbol{q}}=\begin{bmatrix}{\dot{x}},\dot{\theta}\end{bmatrix}^{{\text{T}}}}$,
the Lagrangian ${L=T-U}$ is written as

$L(\boldsymbol{q},\dot{\boldsymbol{q}})=\frac{1}{2}\dot{\boldsymbol{q}}^{{\text{T}}}M(\boldsymbol{q})\dot{\boldsymbol{q}}-U(\boldsymbol{q})$
(125)

where

$M(\boldsymbol{q})=\begin{bmatrix}(m_{c}+m_{p})&m_{p}l\cos(\theta)\\\
m_{p}l\cos(\theta)&m_{p}l^{2}\end{bmatrix},\quad
U(\boldsymbol{q})=m_{p}gl\cos(\theta)$ (126)

The generalized momentum is

$\boldsymbol{p}=\frac{\partial
L}{\partial\dot{\boldsymbol{q}}}=M(\boldsymbol{q})\dot{\boldsymbol{q}}$ (127)

and the Hamiltonian of the system is

$H(\boldsymbol{q},\boldsymbol{p})=\frac{1}{2}\boldsymbol{p}^{{\text{T}}}M(\boldsymbol{q})^{-1}\boldsymbol{p}+U(\boldsymbol{q})$
(128)

Finally, the Hamiltonian dynamics are written as

$\displaystyle\dot{\boldsymbol{q}}$ $\displaystyle=\dfrac{\partial
H}{\partial\boldsymbol{p}}=M(\boldsymbol{q})^{-1}\boldsymbol{p}$ (129)
$\displaystyle\dot{\boldsymbol{p}}$ $\displaystyle=-\dfrac{\partial
H}{\partial\boldsymbol{q}}=-\left(\frac{1}{2}\boldsymbol{p}^{{\text{T}}}\dfrac{\partial
M(\boldsymbol{q})^{-1}}{\partial\boldsymbol{q}}\boldsymbol{p}+\dfrac{\partial
U(\boldsymbol{q})}{\partial\boldsymbol{q}}\right)$ (130)

The parameters of the true system were ${m_{c}=0.8}$, ${m_{p}=0.5}$, ${l=1}$,
and ${g=9.81}$. The training set was generated by uniformly sampling an
increasing number of initial conditions in the interval

$\begin{bmatrix}-2\\\ -\pi\\\ -2\\\ -2\end{bmatrix}\leq\begin{bmatrix}q_{1}\\\
q_{2}\\\ p_{1}\\\ p_{2}\end{bmatrix}\leq\begin{bmatrix}2\\\ \pi\\\ 2\\\
2\end{bmatrix}$ (131)

The number of initial conditions was 15, 31, 63, 127, 255, 511, and 1023. For
each initial condition, the true system was simulated for
${t\in\begin{bmatrix}0,2\end{bmatrix}}$ seconds, with 30 times steps in each
trajectory. The velocities ${\boldsymbol{y}}$ were sampled at each trajectory
point, and zero mean Gaussian noise with standard deviation
${\sigma_{n}=0.01}$ was added to the trajectory and velocity data. A separate
test set was generated uniformly sampling 10 initial conditions in the
interval

$\begin{bmatrix}-2\\\ -\pi\\\ -2\\\ -2\end{bmatrix}<\begin{bmatrix}q_{1}\\\
q_{2}\\\ p_{1}\\\ p_{2}\end{bmatrix}<\begin{bmatrix}2\\\ \pi\\\ 2\\\
2\end{bmatrix}$ (132)

and simulating the true system for ${t\in\begin{bmatrix}0,2\end{bmatrix}}$
seconds, with 30 times steps in each trajectory.

The ${d=50}$ random features were used for the Gaussian model, and ${d=400}$
random features were used for the odd symplectic model, giving an equal number
of model coefficients ${\boldsymbol{\alpha}}$ for each model.

(a) Training MSE

(b) Test MSE

Figure 3: Cart-pole: MSE for the training set and test set, with an increasing
set of training trajectories. Axis are in log-log scale

The final learned models were simulated using the same initial conditions and
time horizon as the true system, and the resulting trajectories were compared
by calculating the MSE for both the training trajectories and the test
trajectories. The experiments showed that the odd symplectic model
outperformed the Gaussian model for both the training set and the test set.
The differences between the two methods was largest when there were few
trajectories in the training set, but the odd symplectic model outperformed
the Gaussian model for all numbers of initial conditions. This can be observed
in Figure 3, where the training MSE and test MSE are shown for each number of
initial conditions in the training set.

###  6.4 Two-link Planar Robot

The two-link planar robot [27], also known as Acrobot or Pendubot when
underactuated, consists of two pendulums linked together. The first link with
uniformly distributed mass ${m_{1}}$ and length ${L_{1}}$ rotates about some
fixed point with angle ${\theta_{1}}$ like a simple pendulum. The second link
with uniformly distributed mass ${m_{2}}$ and length ${L_{2}}$ rotates like a
simple pendulum about the end of link 1 with the angle ${\theta_{2}}$. The
Hamiltonian dynamics are derived for the unactuated system for this task of
system identification. The zero configuration ${\theta_{1}=\theta_{2}=0}$ is
for both links to point directly down. The center of masses of the two links
are at the lengths ${l_{1}}$ and ${l_{2}}$ from their respective pivot points.
In the inertial frame, the positions of the two center of mass points are
given by the kinematics

$\boldsymbol{r}_{1}=\begin{bmatrix}l_{1}\sin\theta_{1}\\\
-l_{1}\cos\theta_{1}\end{bmatrix},\quad\boldsymbol{r}_{2}=\begin{bmatrix}l_{2}\sin\left(\theta_{1}+\theta_{2}\right)+L_{1}\sin\theta_{1}\\\
-l_{2}\cos\left(\theta_{1}+\theta_{2}\right)-L_{1}\cos\theta_{1}\end{bmatrix}$
(133)

where ${\boldsymbol{r}_{1}}$ and ${\boldsymbol{r}_{2}}$ are the positions in
the ${xy}$-plane of the two idealized masses ${m_{1}}$ and ${m_{1}}$,
respectively. The velocities are then

$\displaystyle\boldsymbol{v}_{1}$
$\displaystyle=\begin{bmatrix}l_{1}\omega_{1}\cos\theta_{1}\\\
l_{1}\omega_{1}\sin\theta_{1}\end{bmatrix}$ (134)
$\displaystyle\boldsymbol{v}_{2}$
$\displaystyle=\begin{bmatrix}l_{2}\omega_{1}\cos\left(\theta_{1}+\theta_{2}\right)+l_{2}\omega_{2}\cos\left(\theta_{1}+\theta_{2}\right)+L_{1}\omega_{1}\cos\theta_{1}\\\
l_{2}\omega_{1}\sin\left(\theta_{1}+\theta_{2}\right)+l_{2}\omega_{2}\sin\left(\theta_{1}+\theta_{2}\right)+L_{1}\omega_{1}\sin\theta_{1}\end{bmatrix}$
(135)

Using the mass moment of inertia about the center of mass for a slender rod
given as ${I=\frac{1}{12}mL^{2}}$, the kinetic energy ${T}$ of the system is

$T=T_{1}+T_{2}$ (136)

where

$\displaystyle T_{1}$
$\displaystyle=\frac{1}{2}\left(m_{1}l_{1}^{2}+I_{1}\right)\omega_{1}^{2}$
(137)
$\displaystyle\begin{split}T_{2}&=\frac{1}{2}m_{2}(l_{2}^{2}\omega_{1}^{2}+2l_{2}^{2}\omega_{1}\omega_{2}+2l_{2}L_{1}\omega_{1}^{2}\cos(\theta_{2})+l_{2}^{2}\omega_{2}^{2}\\\
&\quad+2l_{2}L_{1}\omega_{1}\omega_{2}\cos(\theta_{2})+L_{1}^{2}\omega_{1}^{2})+\frac{1}{2}I_{2}(\omega_{1}^{2}+2\omega_{1}\omega_{2}+\omega_{2}^{2})\end{split}$
(138)

The potential energy ${U}$ of the system is then derived using the kinematics
in the vertical direction

$U=g\left(-\left(m_{1}l_{1}+m_{2}L_{1}\right)\cos\left(\theta_{1}\right)-m_{2}l_{2}\cos\left(\theta_{1}+\theta_{2}\right)\right)$
(139)

Defining the generalized coordinate
${\boldsymbol{q}=\begin{bmatrix}\theta_{1},\theta_{2}\end{bmatrix}^{{\text{T}}}}$
and its time derivative
${\dot{\boldsymbol{q}}=\begin{bmatrix}\omega_{1},\omega_{2}\end{bmatrix}^{{\text{T}}}}$,
the Lagrangian ${L=T-U}$ is written as

$L(\boldsymbol{q},\dot{\boldsymbol{q}})=\frac{1}{2}\dot{\boldsymbol{q}}^{{\text{T}}}M(\boldsymbol{q})\dot{\boldsymbol{q}}-U(\boldsymbol{q})$
(140)

where

$M(\boldsymbol{q})=\begin{bmatrix}M_{1}&M_{2}\\\ M_{2}&M_{3}\end{bmatrix}$
(141)

with

$\displaystyle M$
${}_{1}=m_{1}l_{1}^{2}+m_{2}l_{2}^{2}+m_{2}L_{1}^{2}+I_{1}+I_{2}+2m_{2}l_{2}L_{1}\cos(\theta_{2})$
(142) $\displaystyle M$
${}_{2}=m_{2}l_{2}^{2}+I_{2}+m_{2}l_{2}L_{1}\cos(\theta_{2})$ (143)
$\displaystyle M$ ${}_{3}=m_{2}l_{2}^{2}+I_{2}$ (144)

and

$U(\boldsymbol{q})=g\left(-\left(m_{1}l_{1}+m_{2}L_{1}\right)\cos\left(\theta_{1}\right)-m_{2}l_{2}\cos\left(\theta_{1}+\theta_{2}\right)\right)$
(145)

The generalized momentum is

$\boldsymbol{p}=\frac{\partial
L}{\partial\dot{\boldsymbol{q}}}=M(\boldsymbol{q})\dot{\boldsymbol{q}}$ (146)

and the Hamiltonian of the system is

$H(\boldsymbol{q},\boldsymbol{p})=\frac{1}{2}\boldsymbol{p}^{{\text{T}}}M(\boldsymbol{q})^{-1}\boldsymbol{p}+U(\boldsymbol{q})$
(147)

Finally, the Hamiltonian dynamics are written as

$\displaystyle\dot{\boldsymbol{q}}$ $\displaystyle=\dfrac{\partial
H}{\partial\boldsymbol{p}}=M(\boldsymbol{q})^{-1}\boldsymbol{p}$ (148)
$\displaystyle\dot{\boldsymbol{p}}$ $\displaystyle=-\dfrac{\partial
H}{\partial\boldsymbol{q}}=-\left(\frac{1}{2}\boldsymbol{p}^{{\text{T}}}\dfrac{\partial
M(\boldsymbol{q})^{-1}}{\partial\boldsymbol{q}}\boldsymbol{p}+\dfrac{\partial
U(\boldsymbol{q})}{\partial\boldsymbol{q}}\right)$ (149)

The parameters of the true system were ${m_{1}=m_{2}=1}$, ${L_{1}=1}$,
${L_{2}=2}$, ${l_{1}=0.5}$, ${l_{2}=1}$, and ${g=9.81}$. The training set was
generated by uniformly sampling an increasing number of initial conditions in
the interval

$\begin{bmatrix}-\pi\\\ -\pi\\\ -2\\\
-2\end{bmatrix}\leq\begin{bmatrix}q_{1}\\\ q_{2}\\\ p_{1}\\\
p_{2}\end{bmatrix}\leq\begin{bmatrix}\pi\\\ \pi\\\ 2\\\ 2\end{bmatrix}$ (150)

The number of initial conditions was 15, 31, 63, 127, 255, 511, and 1023. For
each initial condition, the true system was simulated for
${t\in\begin{bmatrix}0,1\end{bmatrix}}$ seconds, with 30 times steps in each
trajectory. The velocities ${\boldsymbol{y}}$ were sampled at each trajectory
point, and zero mean Gaussian noise with standard deviation
${\sigma_{n}=0.01}$ was added to the trajectory and velocity data. A separate
test set was generated uniformly sampling 10 initial conditions in the
interval

$\begin{bmatrix}-\pi\\\ -\pi\\\ -2\\\ -2\end{bmatrix}<\begin{bmatrix}q_{1}\\\
q_{2}\\\ p_{1}\\\ p_{2}\end{bmatrix}<\begin{bmatrix}\pi\\\ \pi\\\ 2\\\
2\end{bmatrix}$ (151)

and simulating the true system for ${t\in\begin{bmatrix}0,2\end{bmatrix}}$
seconds, with 30 times steps in each trajectory.

The ${d=100}$ random features were used for the Gaussian model, and ${d=800}$
random features were used for the odd symplectic model, giving an equal number
of model coefficients ${\boldsymbol{\alpha}}$ for each model.

The final learned models were simulated using the same initial conditions and
time horizon as the true system, and the resulting trajectories were compared
by calculating the MSE for both the training trajectories and the test
trajectories. The odd symplectic model outperformed the Gaussian model in the
experiments for both the training set and the test set.

(a) Training MSE

(b) Test MSE

Figure 4: 2-link Robot: MSE for the training set and test set, with an
increasing set of training trajectories. Axis are in log-log scale

As with the cart-pole, there is a trend for the errors to drop as the size of
the training set increases. However, the trend is not as evident, and the
absolute magnitude of the error is larger. This might be due to the chaotic
nature of the 2-link robot combined with the noise added to the training data.
The results from the experiments can be observed in Figure 4, where the
training MSE and test MSE are shown for each number of initial conditions in
the training set.

###  6.5 Varying number of random features

The odd symplectic kernel was compared to its random feature approximation.
The comparison was performed by learning the Hamiltonian dynamics of the
simple pendulum given in (120), using the odd symplectic kernel in (59) and
its random feature approximation in (99).

The training set was generated by randomly sampling ${5000}$ points in the set
${S=\\{\boldsymbol{x}\in{\mathbb{R}}^{2}|q\leq\pi,p\leq 8\\}}$. The velocities
${\boldsymbol{y}}$ were sampled at each point, and zero mean Gaussian noise
with standard deviation ${\sigma_{n}=0.01}$ was added to the trajectory and
velocity data. The learned models were evaluated on the three trajectories
used as the training set in Section 6.2.

The random feature models were learned using an increasing number of random
samples ${\boldsymbol{w}}$, and each model was learned using 50 different sets
of random samples, using the mean trajectory error to evaluate the
performance.

Figure 5: Trajectory error for the odd symplectic kernel and its RFF
approximation. The error for the RFF approximation is the mean error over 20
draws of ${d=\\{10,20,40,80,160,320,640,1280,2560\\}}$ random features.

The results show that the true kernel was more accurate than the random
feature approximation for the trajectories used in the evaluation, and the
error using the random feature decreased exponentially with an increase in the
number of random features $d$. The results are shown in Figure 5. According to
Theorem 12 in [24] the approximation of the two signature functions in (94)
will converge exponentially in $d$, and it follows that the approximation of
the odd symplectic kernel will converge exponentially in $d$. This result
agrees well with the observed exponential convergence of the MSE.

###  6.6 Numerical evaluation

The models were evaluated numerically to investigate the ability of the
learned models to capture the side information of the true systems. The odd
symmetry was evaluated by sampling ${10\,000}$ points in the right half plane
for each of the vector fields, and calculating the odd error given as

$e_{\text{odd}}=\|\boldsymbol{f}(\boldsymbol{x})+\boldsymbol{f}(-\boldsymbol{x})\|$
(152)

where ${\boldsymbol{f}:{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}^{n}}$ is the
dynamical system being evaluated and ${\boldsymbol{x}\in{\mathbb{R}}^{n}}$ is
the sampled point. As the cart-pole and 2-link robot are learned for different
numbers of trajectories, the values corresponding to the maximum mean odd
error were used.

The results in Table 1 document that the learned odd symplectic model enforces
odd symmetry like the true systems, whereas the Gaussian separable model does
not.

Table 1: Odd error ${e_{\text{odd}}}$ for the three dynamical systems | Simple pendulum | Cart-pole | 2-link robot
---|---|---|---
System | Mean | Variance | Mean | Variance | Mean | Variance
True ${e_{\text{odd}}}$ | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
Gaussian sep. ${e_{\text{odd}}}$ | 7.87 | 6.22 | 1.67 | 0.45 | 10.62 | 21.6
Odd symplectic ${e_{\text{odd}}}$ | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00

The learned Hamiltonian in (114) for the learned odd symplectic models were
compared to their corresponding real Hamiltonians in (119), (128), and (147),
over the test trajectories. For the Cart-pole and the 2-link robot, the
presented values were selected by selecting for the maximum variance across
all test trajectories.

The results in Table 2 demonstrate that the value of the learned Hamiltonian
${\hat{H}(\boldsymbol{x})}$ has a constant offset from the true Hamiltonian
${H(\boldsymbol{x})}$. This agrees with the fact that the potential energy’s
zero potential cannot be expected to be the same for the learned and true
systems. It is seen that the value of the learned Hamiltonian is constant in
agreement with (36) since the system is unforced and independent of time. This
is reflected in the variance of both ${H(\boldsymbol{x})}$ and
${\hat{H}(\boldsymbol{x})}$. Noting that these are results from numerical
simulations, the results indicate that the Hamiltonian mechanics are captured
in the learned odd symplectic model.

Table 2: Hamiltonian for the three dynamical systems over the test trajectories | Simple pendulum | Cart-pole | 2-link robot
---|---|---|---
Hamiltonian | Mean | Variance | Mean | Variance | Mean | Variance
Real ${H(\boldsymbol{x})}$ | $9.81$ | $4.99\cdot 10^{-15}$ | $7.64$ | $3.3\cdot 10^{-7}$ | $22.43$ | $1.2\cdot 10^{-10}$
Learned ${\hat{H}(\boldsymbol{x})}$ | $28.93$ | $4.08\cdot 10^{-15}$ | $0.05$ | $1.5\cdot 10^{-10}$ | $11.72$ | $7.8\cdot 10^{-08}$

##  7 Conclusion

A specialized kernel enforcing side information relating to Hamiltonian
dynamics and odd symmetry has been presented. The odd symplectic kernel was
developed, approximated using random features, and utilized in three
comparative experiments. By enforcing the side information through the kernel,
the closed-form solution to the learning problem is retained and the side
information is enforced for the whole domain of the learned function. This
stands in contrast to the approach of enforcing the side information through
the use of constraints in a constrained optimization problem, enforcing the
side information only on selected points. Through comparative experiments we
have demonstrated that the proposed kernel outperforms a more standard kernel
ridge regression and Gaussian kernel, as the error over both the training set
and a separate test set is lower, indicating a more accurate and generalized
learned model.

###  7.1 Future work

A challenge with learning Hamiltonian dynamics is the potential lack of data
for the generalized momenta and their derivatives. As a result, the method
should be extended so that it can be applied using only data for the
generalized coordinates and velocities. Furthermore, control-oriented learning
could be studied using the proposed kernel.

## References

  * [1] S. L. Brunton, J. N. Kutz, Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control, 2nd Edition, Cambridge University Press, 2022.  doi:10.1017/9781009089517.
  * [2] V. Sindhwani, S. Tu, M. Khansari, Learning Contracting Vector Fields For Stable Imitation Learning, arXiv preprint arXiv:1804.04878 [cs.RO] (2018).  doi:10.48550/arXiv.1804.04878.
  * [3] A. A. Ahmadi, B. E. Khadir, Learning Dynamical Systems with Side Information, in: Proceedings of the 2nd Conference on Learning for Dynamics and Control, Vol. 120 of Proceedings of Machine Learning Research, PMLR, 2020, pp. 718–727.
  * [4] S. Singh, S. M. Richards, V. Sindhwani, J.-J. E. Slotine, M. Pavone, Learning stabilizable nonlinear dynamics with contraction-based regularization, The International Journal of Robotics Research 40 (10-11) (2021) 1123–1150.  doi:10.1177/0278364920949931.
  * [5] M. Revay, I. Manchester, Contracting Implicit Recurrent Neural Networks: Stable Models with Improved Trainability, in: Proceedings of the 2nd Conference on Learning for Dynamics and Control, Vol. 120 of Proceedings of Machine Learning Research, PMLR, 2020, pp. 393–403.
  * [6] M. Krejnik, A. Tyutin, Reproducing Kernel Hilbert Spaces With Odd Kernels in Price Prediction, IEEE Transactions on Neural Networks and Learning Systems 23 (10) (2012) 1564–1573.  doi:10.1109/TNNLS.2012.2207739.
  * [7] C.-A. Cheng, H.-P. Huang, Learn the Lagrangian: A Vector-Valued RKHS Approach to Identifying Lagrangian Systems, IEEE Transactions on Cybernetics 46 (12) (2016) 3247–3258.  doi:10.1109/TCYB.2015.2501842.
  * [8] M. Ahmadi, U. Topcu, C. Rowley, Control-Oriented Learning of Lagrangian and Hamiltonian Systems, in: American Control Conference (ACC), 2018, pp. 520–525.  doi:10.23919/ACC.2018.8431726.
  * [9] S. Greydanus, M. Dzamba, J. Yosinski, Hamiltonian Neural Networks, in: Proceedings of the Conference on Neural Information Processing Systems (NeurIPS), 2019.
  * [10] Y. D. Zhong, B. Dey, A. Chakraborty, Symplectic ODE-Net: Learning Hamiltonian Dynamics with Control, in: Proceedings of the International Conference on Learning Representations (ICLR), 2020.
  * [11] Z. Chen, J. Zhang, M. Arjovsky, L. Bottou, Symplectic Recurrent Neural Networks, in: International Conference on Learning Representations, 2020.
  * [12] M. Khosravi, R. S. Smith, Nonlinear System Identification With Prior Knowledge on the Region of Attraction, IEEE Control Systems Letters 5 (3) (2021) 1091–1096.  doi:10.1109/LCSYS.2020.3005163.
  * [13] A. J. Thorpe, C. Neary, F. Djeumou, M. M. K. Oishi, U. Topcu, Physics-Informed Kernel Embeddings: Integrating Prior System Knowledge with Data-Driven Control, arXiv preprint arXiv:2301.03565 [eess.SY] (2023).  doi:10.48550/arXiv.2301.03565.
  * [14] H. J. van Waarde, R. Sepulchre, Kernel-Based Models for System Analysis, IEEE Transactions on Automatic Control 68 (9) (2023) 5317–5332.  doi:10.1109/TAC.2022.3218944.
  * [15] T. Smith, O. Egeland, Learning of Hamiltonian Dynamics with Reproducing Kernel Hilbert Spaces, arXiv preprint arXiv:2312.09734 [cs.RO] (2023).  doi:10.48550/arXiv.2312.09734.
  * [16] N. Aronszajn, Theory of Reproducing Kernels, Transactions of the American Mathematical Society 68 (3) (1950) 337–404.
  * [17] C. A. Micchelli, M. Pontil, On Learning Vector-Valued Functions, Neural Computation 17 (1) (2005) 177–204.  doi:10.1162/0899766052530802.
  * [18] R. Brault, M. Heinonen, F. Buc, Random Fourier Features For Operator-Valued Kernels, in: Proceedings of The 8th Asian Conference on Machine Learning, Vol. 63 of Proceedings of Machine Learning Research, PMLR, 2016, pp. 110–125.
  * [19] N. M. Boffi, S. Tu, J.-J. E. Slotine, Nonparametric adaptive control and prediction: theory and randomized algorithms, Journal of Machine Learning Research 23 (281) (2022) 1–46.
  * [20] H. Q. Minh, V. Sindhwani, Vector-valued Manifold Regularization, in: Proceedings of the 28’th International Conference on Machine Learning, 2011.
  * [21] H. Goldstein, C. Poole, J. Safko, Classical Mechanics, 3rd Edition, Addison Wesley, 2002.
  * [22] E. Hairer, G. Wanner, C. Lubich, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd Edition, Springer Berlin, Heidelberg, 2006.  doi:10.1007/3-540-30666-8.
  * [23] E. J. Fuselier, Refined Error Estimates for Matrix-Valued Radial Basis Functions, PhD Thesis, Texas A&M University (2006).
  * [24] H. Q. Minh, Operator-Valued Bochner Theorem, Fourier Feature Maps for Operator-Valued Kernels, and Vector-Valued Learning, arXiv preprint arXiv:1608.05639 [cs.LG] (2016).  doi:10.48550/arXiv.1608.05639.
  * [25] D. E. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning, 1st Edition, Addison-Wesley Longman Publishing Co., Inc., USA, 1989.
  * [26] R. Kohavi, A Study of Cross-Validation and Bootstrap for Accuracy Estimation and Model Selection, in: Proceedings of the 14th International Joint Conference on Artificial Intelligence, IJCAI-95, 1995, p. 1137–1143.
  * [27] R. Olfati-Saber, Nonlinear Control of Underactuated Mechanical Systems with Application to Robotics and Aerospace Vehicles, PhD Thesis, Massachusetts Institute of Technology (2001).