# Fixed-kinetic Neural Hamiltonian Flows for
enhanced interpretability and reduced complexity

Vincent Souveton  Arnaud Guillin  Jens Jasche  Guilhem Lavaux  Manon Michel

###### Abstract

Normalizing Flows (NF) are Generative models which are particularly robust and
allow for exact sampling of the learned distribution. They however require the
design of an invertible mapping, whose Jacobian determinant has to be
computable. Recently introduced, Neural Hamiltonian Flows (NHF) are based on
Hamiltonian dynamics-based Flows, which are continuous, volume-preserving and
invertible and thus make for natural candidates for robust NF architectures.
In particular, their similarity to classical Mechanics could lead to easier
interpretability of the learned mapping. However, despite being Physics-
inspired architectures, the originally introduced NHF architecture still poses
a challenge to interpretability. For this reason, in this work, we introduce a
fixed kinetic energy version of the NHF model. Inspired by physics, our
approach improves interpretability and requires less parameters than
previously proposed architectures. We then study the robustness of the NHF
architectures to the choice of hyperparameters. We analyze the impact of the
number of leapfrog steps, the integration time and the number of neurons per
hidden layer, as well as the choice of prior distribution, on sampling a
multimodal 2D mixture. The NHF architecture is robust to these choices,
especially the fixed-kinetic energy model. Finally, we adapt NHF to the
context of Bayesian inference and illustrate our method on sampling the
posterior distribution of two cosmological parameters knowing type Ia
supernovae observations.

Generative models, Normalizing flows, Hamiltonian flows, Sampling

##  1 Introduction

Generative models are widely used for sampling high-dimensional probability
distributions with applications from molecular biology (Lopez et al., 2020) to
cosmology (Rodriguez et al., 2018) or medical science (Frazer et al., 2021).
Traditional architectures like Generative Adversarial Networks (Goodfellow et
al., 2014) have shown impressive results in image generation but since their
adversarial loss seeks a saddle point rather than a local minimum, GANs are
notoriously hard to train and may suffer from mode-collapse (Lin et al., 2018;
Arjovsky & Bottou, 2017; Berard et al., 2020). More robust techniques like
Normalizing Flows (NF) were developed (Tabak & Vanden-Eijnden, 2010; Dinh et
al., 2014; Rezende & Mohamed, 2015). NF consist in training a neural network
to map a simple prior distribution onto the desired target through a chain of
invertible transformations. They come with interesting characteristics, such
as stability and correctness, see for example (Papamakarios et al., 2022). The
main limitation comes from the design of an invertible function for the
mapping. In particular, computing the Jacobian determinant in the change of
variable formula may be costly. Furthermore, while explainability is now under
a growing concern within the community (Gilpin et al., 2018), and in
particular regarding applications in natural science, the transformation
learned by NF models is commonly hard to interpret.

Exploiting the Newtonian evolution in classical Mechanics, Neural Hamiltonian
Flows (NHF) (Toth et al., 2020) are NF models that use Hamiltonian
transformations. The Hamiltonian of a system is composed of a kinetic and
potential energy terms, which sets its dynamical evolution, which is
reversible and has a Jacobian determinant equal to one. They come with
performance similar to the ones obtained with Real-NVPs in sampling 2D
distributions (Toth et al., 2020). Being Physics-driven models, they are
expected to enhance interpretability. However, the NHF architecture is made of
four neural networks black-boxes that render difficult the interpretation of
the learned dynamics. Furthermore, even if they have been numerically shown to
transfer multimodality from the target distribution to the potential energy in
some cases (Toth et al., 2020), this property is not guaranteed. However, the
Hamiltonian formalism could offer a powerful formalism to design more easily
interpreted architectures. Thus, some recent works have proposed to exploit
the Hamiltonian properties to design interpretable flows that are invariant
under symmetrical transformations (Jimenez Rezende et al., 2019).

In this work, we focus on the transfer of the negative logarithm of the target
distribution into the learned potential, which should be the case in any
physical systems. This leads us to propose a fixed-kinetic version of NHF
where the kinetic-energy term in the Hamiltonian is set to its true form in
classical mechanics, where momenta follow a Gaussian distribution. We also
discuss the impact of the hyperparameters on such transfer and study the
overall robustness of NHF towards the choice of the numerical integration
scheme and of the prior distributions. In more details, this work provides the
following four main contributions:

  * •

We introduce a fixed-kinetic version of NHF that, thanks to the Hamiltonian
evolution, enhances interpretability of the model and does so at a cheaper
computational price.

  * •

We analyze the effect of multiple parameters of the architecture on sampling a
2-dimensional multimodal distribution and show that NHF is robust to the
choice of hyper-parameters, especially the fixed-kinetic model.

  * •

We show that the choice of prior has an influence on the learned dynamics. The
fixed-kinetic model also allows better robustness to the choice of prior
distribution.

  * •

Finally, aside from Generative modeling, flow-based models have been shown to
be suited as inference algorithms (Rezende & Mohamed, 2015; Winkler et al.,
2019). We test a framework for Bayesian inference using NHF and present
numerical experiments for inferring cosmological parameters from astronomical
observations. The methodology we propose is inspired by Boltzmann generators
(Noé et al., 2019).

The manuscript is organized as follows: in Section 2, we present and review
the related works. In Section 3, we describe the theoretical framework and the
practical implementation of NHF. In Section 4, we discuss the choice of models
for the kinetic energy and introduce a fixed-kinetic version for enhanced
interpretability and reduced complexity. Section 5 discusses the choices one
can make to maximize expressivity given a fixed computational budget, with
tests on a 2D Gaussian mixture and analysis of the impact of the Leapfrog-
hyperparameters and model complexity. We also show how the choice of base
distribution affects the learned energies and thus the interpretability of the
model. Finally, in Section 6, we adapt NHF for Bayesian inference and
illustrate our method on a toy model from cosmology.

##  2 Related works

Generative models. These methods make the assumption that data can be
represented by an underlying probability distribution. The aim of generative
models is to learn such distribution in order to produce original samples
similar to the training dataset. Various architectures have been presented
such as Generative Adversarial Networks (Goodfellow et al., 2014) or diffusion
networks (Sohl-Dickstein et al., 2015). In this paper, we will focus on
Normalizing Flows techniques (Tabak & Vanden-Eijnden, 2010; Dinh et al., 2014;
Rezende & Mohamed, 2015) as a way of smoothly transforming a simple prior
distribution into the target posterior.

Learning Hamiltonians. Learning Hamiltonian potentials, i.e. physical
conserved quantities, is a first step towards a better understanding of the
physical processes that have governed the data generation. Multiple
architectures have been proposed, such as Hamiltonian Neural Networks
(Greydanus et al., 2019) or Hamiltonian Generative Networks (Toth et al.,
2020). These methods parameterize the Hamiltonian of the system with neural
networks and come with useful properties such as exact reversibility and
smoothness. They have inspired applications from domain translation (Menier et
al., 2022) to fault-detection in industry (Shen et al., 2023). It is worth
mentioning that they can be combined with Markov-Chain Monte Carlo (MCMC)
methods, for instance as proposals in the Hamiltonian Monte Carlo (HMC)
algorithm (Duane et al., 1987; Dhulipala et al., 2022). In this work, we learn
artificial Hamiltonians for sampling and our goal is to extract the negative
logarithm of the target distribution into the potential.

Neural Hamiltonian Flows. NHF is a Normalizing Flow architecture derived as a
modification from Hamiltonian Generative Networks. As generative models, they
have been used to sample from 2D distributions (Toth et al., 2020). The
Hamiltonian formalism is also a natural framework for exploiting symmetries
inside data and build interpretable flows (Jimenez Rezende et al., 2019). In
this paper, we discuss its robustness with respect to hyperparameters and
choice of prior distribution. Also, we propose an alternative version of NHF
to enhance interpretability while reducing the complexity of the model.

Inference with NHF. Inferring the posterior distribution of parameters
conditioned on data is an important problem in science. Traditional MCMC
methods (Robert & Casella, 2004) are very popular because they come with
guarantees in terms of convergence and many progress have been made regarding
their tuning (Homan & Gelman, 2014; Carpenter et al., 2017). NF architectures
have also been proposed in this framework (Rezende & Mohamed, 2015; Winkler et
al., 2019). Here, we adapt NHF to sampling Bayesian posterior distributions by
transforming the prior distribution into the posterior with no access to
samples from the target.

Explainable AI. XAI deals with the problem of understanding the decisions made
by an Artificial Intelligence (Samek & Müller, 2019). Indeed, complex
architectures made of multiple (deep) neural network are often easier to train
than to understand. Some solutions involve surrogate techniques (Ribeiro et
al., 2016), local perturbations (Ancona et al., 2022) or meta-explanations
(Lapuschkin et al., 2019). Including physical prior knowledge into neural
networks may be another solution to understand the model (Raissi et al., 2019;
Toth et al., 2020). In this work, we build on that idea and try to make the
model as explainable as possible by fixing the kinetic energy of the model and
thus enforcing a classical Mechanics knowledge into the architecture.

##  3 Normalizing Flows with Hamiltonian transformations

###  3.1 Normalizing Flows

Normalizing flows are generative models that map a complex target distribution
$\pi$ onto a known prior distribution $\pi_{0}$ from which it is easy to
sample (Papamakarios et al., 2022). This mapping is a series of smooth
invertible transformations. Once the model is trained, one can reverse the
learned dynamics to generate samples from the target distribution starting
from the prior. If $X=\mathcal{T}(Z)$, where $\mathcal{T}$ is a
$\mathcal{C}^{1}$-diffeomorphism and $Z\sim\pi_{0}$, then the density followed
by $X$ reads $m(x)=\pi_{0}\left(\mathcal{T}^{-1}(x)\right)\times\left|\det
J_{\mathcal{T}^{-1}}(x)\right|,$ where $\det J_{f^{-1}}(x)$ is the Jacobian
determinant of the inverse transformation of $f$ evaluated in $x$. For more
expressivity, it is possible to chain several invertible transformations
$x=\mathcal{T}_{L}\circ...\circ\mathcal{T}_{1}(z)$. In this case, the change
of variable formula is now:
$m(x)=\pi_{0}\left(\mathcal{T}^{-1}_{1}\circ...\circ\mathcal{T}^{-1}_{L}(x)\right)\times\prod_{k=1}^{L}\left|\det
J_{\mathcal{T}_{k}^{-1}}(x)\right|.$ The model parameters to optimize are
denoted $\Theta$. The goal is to minimize the Kullback-Leibler divergence
between the target distribution $\pi$ and the model distribution $m$ with
respect to $\Theta$, i.e. minimizing:

$\displaystyle\mathcal{L}(\Theta)=\mathbf{E}_{\pi}\left[\log\pi(X)-\log
m(X;\Theta)\right]$
$\displaystyle=-\mathbf{E}_{\pi}\Bigg{[}\log\pi\left(\mathcal{T}_{1}^{-1}\circ...\circ\mathcal{T}_{T}^{-1}(X;\Theta)\right)$
$\displaystyle\qquad\left.+\sum_{k=1}^{T}\left|\det
J_{\mathcal{T}_{k}^{-1}}(X)\right|\right]+C.$

As we seek to find the minimum of this quantity, $C$ can be discarded. Then,
one can use samples $x_{1},...,x_{m}$ from the target distribution for
estimating the above loss via Monte Carlo. The function to minimize becomes:

$\begin{split}\small
L(\Theta)=-\frac{1}{m}\sum_{k=1}^{m}&\Bigg{[}\log\pi\left(\mathcal{T}_{1}^{-1}\circ...\circ\mathcal{T}_{L}^{-1}(x_{k};\Theta)\right)\\\
&\left.+\sum_{k=1}^{T}\left|\det
J_{\mathcal{T}_{k}^{-1}}(x_{k})\right|\right].\end{split}$ (1)

At this point, transformations are to some extent arbitrary. Any smooth
invertible function would be suited but the main computational cost comes from
the Jacobian determinants. The first goal is then to reduce this computational
cost, and the second one to enhance interpretability, by a proper choice of
the transformation, and of its induced inverse.

###  3.2 Neural Hamiltonian Flows

To alleviate these issues, Neural Hamiltonian Flows (NHF, Toth et al., 2020)
is a NF technique that uses a series of Hamiltonian transformations as
normalizing flows. In classical Mechanics, a system is fully described by its
coordinates $(\mathbf{q},\mathbf{p})$ in phase-space. From that description,
it is possible to define a scalar quantity called a Hamiltonian (Landau &
Lifshitz, 1982). It can be seen as the total energy of the system and, in this
paper, we make the assumption that it is written as the sum of a potential
energy $V$, solely depending on the generalized positions $\mathbf{q}$, and a
kinetic energy $K$, solely depending on the momenta $\mathbf{p}$. The system
evolves in phase-space following Hamilton’s equations that read:

$\frac{d\mathbf{q}}{dt}=\frac{\partial H}{\partial\mathbf{p}},\
\frac{d\mathbf{p}}{dt}=-\frac{\partial H}{\partial\mathbf{q}}.$ (2)

Hamiltonian transformations present at least two main advantages that make
them suited for normalizing flows:

  * •

they are invertible by construction and inversion is easy by using a classical
numerical integrator, i.e. just reversing the speed;

  * •

their Jacobian determinant is equal to 1, removing the necessity to compute
such determinant for each transformation.

Numerically, the continuous solution can be approached by a symplectic,
invertible and stable integrator as a Leapfrog:

$\left\\{\begin{array}[]{rll}\mathbf{p}_{n+\frac{1}{2}}&=&\mathbf{p}_{n}-\nabla
V(\mathbf{q}_{n})\times\frac{dt}{2},\\\
\mathbf{q}_{n+1}&=&\mathbf{q}_{n}+\nabla K(\mathbf{p}_{n+\frac{1}{2}})\times
dt,\\\ \mathbf{p}_{n+1}&=&\mathbf{p}_{n+\frac{1}{2}}-\nabla
V(\mathbf{q}_{n+1})\times\frac{dt}{2}.\end{array}\right.$ (3)

NHF is trained on a dataset consisting in realizations from the target
distribution. In order to simulate a Hamiltonian dynamics, one must extend the
position space in which live the samples into the phase space, by adding
artificial momenta: this is the role of the Encoder. The dynamics is
integrated in phase-space with the Leapfrog integrator (Toth et al., 2020).
More precisely, during training, NHF takes batches of $\mathbf{q}_{T}$ from
the training dataset as inputs. For each $\mathbf{q}_{T}$, one
$\mathbf{p}_{T}$ is drawn from a Gaussian distribution whose mean
$\mu(\mathbf{q}_{T})$ and standard deviation $\sigma(\mathbf{q}_{T})$ depend
on the $\mathbf{q}_{T}$. The resulting point in phase-space is then evolving
through a series of $L$ Leapfrog steps with integration timestep $dt$. The
outputs consist in the final position $\mathbf{q}_{0}$ and momenta
$\mathbf{p}_{0}$, as well as the initial mean $\mu(\mathbf{q}_{T})$, standard
deviation $\sigma(\mathbf{q}_{T})$ and $\mathbf{p}_{T}$ that are used in the
loss computation. Once trained, one can easily define a sampling function that
transforms $\mathbf{q}_{0}$, $\mathbf{p}_{0}$ into $\mathbf{q}_{T}$, by
changing the sign of integration timestep and moving the system through the
learned dynamics. The architecture is illustrated in Figure 1.

Figure 1: NHF architecture

Following the previous notations, let us call $f(.|\mathbf{q}_{T})$ the
density of a normal distribution
$\mathcal{N}(\mu(\mathbf{q}_{T}),\sigma(\mathbf{q}_{T})^{2})$, and
$\mathcal{T}^{-1}$ the backward transformation of phase-space performed by NHF
i.e.
$\mathcal{T}^{-1}(\mathbf{q}_{T},\mathbf{p}_{T})=(\mathbf{q}_{0},\mathbf{p}_{0})$.
Also, we denote $\Pi_{0}$ the joint distribution of $\mathbf{q}_{0}$,
$\mathbf{p}_{0}$. By adding artificial momenta $\mathbf{p}_{T}$ (Toth et al.,
2020), the distribution modeled by our NHF is $m(\mathbf{q}_{T})=\int
M(\mathbf{q}_{T},\mathbf{p}_{T})d\mathbf{p}_{T}=\int\Pi_{0}(\mathcal{T}^{-1}(\mathbf{q}_{T},\mathbf{p}_{T}))d\mathbf{p}_{T}$.
This integral is intractable so instead one can maximize the following ELBO:

$\mathcal{L}(\mathbf{q}_{T})=\mathbf{E}_{f}\left[\log\pi(\mathcal{T}^{-1}(\mathbf{q}_{T},\mathbf{p}_{T}))-\log
f(\mathbf{p}_{T}|\mathbf{q}_{T})\right].$ (4)

This quantity is approximated via Monte Carlo integration.

Once it has learned the transformation, one can reverse the sign of timesteps
and use the same potentials to transform the base distribution into the target
distribution.

The first part of the architecture consists in adding artificial momenta to
simulate a Hamiltonian dynamics. This is done by the Encoder. Here, $\mu$ and
$\sigma$ are approximated by two neural networks. As for the Hamiltonian
transformations, they are made by chaining Leapfrog steps. To do so, one must
design the potential energy $V$ and the kinetic energy $K$ of the system. In
previous work (Toth et al., 2020), authors have proposed to parameterize each
potential by a neural network. We will discuss this choice in the following
section. For now, let us highlight the fact that the choice of neural network
for the energies or within the Encoder is without constraint because the
properties needed for normalizing flows are ensured by the use of a Leapfrog
scheme. This is an illustration of the model flexibility.

##  4 Designing the kinetic energy for NHF

MLP-kinetic NHF. If the kinetic energy is chosen to be a MLP (Toth et al.,
2020), then the model contains two black-boxes that are not easy to interpret
a priori, namely the kinetic and potential energies $K$ and $V$. In
particular, when sampling a multimodal distribution from a unimodal prior, one
may have some troubles understanding where the transfer of multimodality
occurs in the process.

Fixed-kinetic NHF. By fixing the kinetic energy inside NHF, we expect to gain
interpretability on the learned flow by forcing the latter to obey some
Physics principles. Ideally, one would like to find a way of enforcing these
energies to be classical from a classical Physics perspective, i.e. making
learned kinetic to be a quadratic form and learned potential to be the
negative logarithm of the target distribution (or an approximation). This is
an important aspect because we would like to keep track of the transformation
dynamics. This way, we can play with the familiar classical mechanics
framework and interpret the model more easily. By doing so, both
interpretability and maniability (less parameters) are gained through the
Fixed-kinetic NHF. In this model, $K$ is no longer a MLP but a fixed function.
In practice, one can impose the kinetic energy to be written as:

$K(\mathbf{p})=\frac{1}{2}\mathbf{p}^{T}{\cal M}\mathbf{p},$ (5)

with ${\cal M}$ a mass matrix. We impose ${\cal M}$ to be a positive symmetric
definite matrix, as it is the case in classical mechanics. For instance, in
2D, it comes down to optimizing three real parameters, rather than a complete
network. Starting from ($\mathbf{q}_{0},\mathbf{p}_{0})$ drawn from a unimodal
prior distribution and imposing a quadratic kinetic energy significantly
reduces the number of possibilities for the potential energy if one wants to
recover a multimodal $\mathbf{q}_{T}$. Let us examine what happens for a
single Leapfrog transformation. We have $\nabla K(\mathbf{p})={\cal
M}\mathbf{p}$. Using the first step of the Leapfrog algorithm (3), the final
expression for $(\mathbf{q}_{1},\mathbf{p}_{1})$ reads:

$\displaystyle\mathbf{q_{1}}$ $\displaystyle=$
$\displaystyle\mathbf{q_{0}}+{\cal M}\left(\mathbf{p}_{0}-\nabla
V(\mathbf{q}_{0})\frac{dt}{2}\right)dt$ $\displaystyle\mathbf{p}_{1}$
$\displaystyle=$ $\displaystyle\mathbf{p_{0}}-\nabla
V(\mathbf{q_{0}})\frac{dt}{2}$ $\displaystyle-\nabla
V\left(\mathbf{q_{0}}+{\cal M}\left(\mathbf{p}_{0}-\nabla
V(\mathbf{q}_{0})\frac{dt}{2}\right)dt\right)\frac{dt}{2}.$

So, if $V$ is quadratic, then the overall transformation is linear and
multimodality cannot be obtained from the original Gaussian measure. Thus,
multimodality can be achieved only through the potential $V$ in order to
sample from the correct target.

It is possible to create many different versions of NHF by fixing its kinetic
energy, which thus leads to analysis of the learned dynamics in light of a
different Physics framework. One could, for instance, use a relativistic
kinetic energy instead of a classical one. We now numerically show how the
classical choice for kinetic energy yields an interpretable potential $V$ and
such in a robust manner.

##  5 Robustness and optimization of NHF

We present the results of numerical experiments for analyzing the robustness
of NHF. The experiments consist in sampling a 2D Gaussian mixture with 9 modes
(see figure 2). Such toy example, similarly studied in (Toth et al., 2020), is
interesting in the sense that it will allow us to discuss various aspects,
from memory usage to interpretability. Also, traditional generative models
like GAN may suffer from mode-collapse problems even in simple multimodal 2D
settings (Eghbal-zadeh et al., 2019). It should be noted that we never
observed such issue with NHF in all the experiments that have been conducted.

In each model, $\mu$ and $\sigma$ are MLPs with size $(D,H,H,D)$. The energies
are MLPs with size $(D,H,H,1)$, $D$ being the dimension of data, except for
the fixed-kinetic version of NHF where $K$ is a fixed quadratic form. Models
are trained on a 5,000 points dataset with minibatches of size 512. Weights
and biases are optimized with Adam (Kingma & Ba, 2015), setting the learning
rate to $5\times 10^{-4}$.

###  5.1 Impact of Leapfrog-hyperparameters and model complexity

We first discuss the effect of Leapfrog-hyperparameters $L$ (number of
Leapfrog steps) and $T=L\times dt$ (integration time) on the optimization, but
also the impact of the model complexity. The latter is governed by the total
number of neurons in the model, this number being an increasing function of
$H$, the number of neurons per hidden layer in each MLP of the model. If the
model is complex enough, we expect it to learn how to adjust to the number of
Leapfrog steps and choice of integration time. If not, we expect the model to
perform better by increasing the number of computations, i.e. increasing $L$.
As $L$ increases, we do not expect to see any difference as the number of
computations becomes sufficient to solve the problem.

We tested both fixed-kinetic and MLP-kinetic NHF with various choices of $L$,
$T$ and $H$. We use a soft-uniform prior $\propto s(x+3)s(-x+3)$, where $s$ is
the sigmoid function. The corresponding loss decays are illustrated in Figure
3 and additional details can be found in Appendix A.

Figure 2: Density estimation of the 2D multimodal target distribution that is
used in Section 5 with its marginals. It consists in 9 equally-weighted
Gaussians with same covariance matrix $0.5^{2}I_{2}$.  Figure 3: Training loss
as a function of epochs for models with different $H$ (number of neurons per
hidden layer in each neural network of the model), $L$ (number of Leapfrog
steps) and $T$ (integration time). First row: MLP-kinetic NHF; second row:
fixed-kinetic model. From left to right: $H=128,\ 8$. Models were trained for
$15,000$ epochs on a 5,000-point dataset with minibatches of size $512$.

Firstly, the fixed-kinetic NHF model is shown to be more robust than the MLP-
kinetic one to the choices of $L$ and $T$, at fixed $H$, as discrepancy in the
loss decay more clearly appears only with $H=8$. Then, regarding the tuning of
the Leapfrog scheme, at fixed-integration time, models with $L=1$ always reach
higher final value of the loss, this effect being less visible with the fixed-
kinetic model. Increasing the number of leapfrog steps leads to better final
performance of the model even if this effect disappears once the number of
Leapfrog steps gets sufficient and no further expressivity can be achieved.
Finally, as for the effect of integration time $T$, it barely appears for the
fixed-kinetic model, showing that the latter efficiently adjusts to this
parameter. As for the MLP-kinetic model, the effect of the integration time is
clearer, but mostly at $H=8$, where performance improves for $T=1,10$ compared
to $T=0.1$. Overall, as the number of parameters in the model is increasing
($H=32$ or $128$), the impact of the choice of the integration time is very
limited.

Thus, for a given training dataset, there are basically four hyperparameters
that require tuning: three are usual in learning (minibatch size, learning
rate and number of neurons per hidden layer, i.e. number of learning
parameters of the model) and only one is specific to NHF (minimum number of
Leapfrog steps), whose tuning is especially less sensitive when using the
fixed-kinetic variant.

###  5.2 Impact of the prior distribution on the learned dynamics

Now we illustrate the impact of the prior choice on the transfer of
characteristics of the target distribution on the potential $V$, especially
regarding the multimodality nature. All models were trained for 15,000 epochs
using $H=128$, $T=1$ and $L=10$, with a 5,000 points training dataset.

In terms of sampling, all considered schemes recover the nine correct modes
from the target distribution, as illustrated in Figure4. We now consider the
learned potential $V$, as illustrated in Figure5. As the Hamiltonian evolution
only involves its derivative, we also include a shifted version in Figure6.
When choosing a relatively flat soft-uniform prior distribution that covers
the target region, multimodality transfers to the potential energy for both
fixed-kinetic and MLP-kinetic NHF. The potential exhibits indeed local extrema
centered at the modes of the target, which can either be minima or maxima for
the MLP-kinetic NHF but are minima for the fixed-kinetic one. Indeed, with a
MLP-kinetic model, the orientation of the learned energies may change from one
numerical experiment to another, as we do not enforce the positiveness of the
output of $V$ and $K$ making this variant of NHF indifferent to the sign of
the learned energies. Finally, similar results were obtained using a Gaussian
base $\mathcal{N}(0,2.5^{2}I_{2})$ with variance large enough to cover the
support of the target distribution, which stresses the impact of the spatial
expansion rather than the nature of the prior distribution.

On the other hand, when choosing a ”peaked” prior distribution
$\mathcal{N}(0,I_{2})$, we remark that, for the MLP-kinetic NHF, the momenta
$\mathbf{p}_{T}$ generated by the Encoder inherit from the multimodality of
the target distribution, with the same number of modes. As a consequence,
learned energies are different from the classical Physics ones. Learned
potential energy is different from one model to another but in the case of
fixed-kinetic NHF, multimodality is always transferred to the potential
energy, showing the robustness of the model to the choice of prior
distribution.

Thus, using a fixed-kinetic model allows to more robustly transfer important
properties of the target distribution into the learned potential. For both
models though, choosing a wide prior that covers the support of the target
distribution is improving such transfer. Finally, when the learned potential
is not multimodal, it is an indication that multimodality has been transferred
instead to the artificial momenta $\mathbf{p}_{T}$ generated by the Encoder.

Figure 4: Density estimation of the samples produced by the six models
previously defined. First row: MLP-kinetic NHF; second row: Fixed-kinetic
model. From left to right: wide soft-uniform prior; peaked Gaussian prior;
wide Gaussian prior. Models were trained for $15,000$ epochs on a 5,000-point
dataset with minibatches of size $512$, $H=128$, $L=10$ and $T=1$. Figure 5:
Potential energies learned by the six models previously defined. First row:
MLP-kinetic NHF; second row: Fixed-kinetic model. From left to right: wide
soft-uniform prior; peaked Gaussian prior; wide Gaussian prior. Models were
trained for $15,000$ epochs on a 5,000-point dataset with minibatches of size
$512$, $H=128$, $L=10$ and $T=1$. Figure 6: Shifted potential energies learned
by the six models previously defined. First row: MLP-kinetic NHF; second row:
Fixed-kinetic model. From left to right: wide soft-uniform prior; peaked
Gaussian prior; wide Gaussian prior. Models were trained for $15,000$ epochs
on a 5,000-point dataset with minibatches of size $512$, $H=128$, $L=10$ and
$T=1$. Figure 7: Density estimation of the artificial momenta drawn by the
Encoder for the six models previously defined. First row: MLP-kinetic NHF;
second row: Fixed-kinetic model. From left to right: wide soft-uniform prior;
peaked Gaussian prior; wide Gaussian prior. Models were trained for $15,000$
epochs on a 5,000-point dataset with minibatches of size $512$, $H=128$,
$L=10$ and $T=1$.

##  6 Adapting NHF for Bayesian inference

###  6.1 Methodology, derivation of the new loss function

We present how NHF can be used to perform Bayesian inference. It consists in
using Hamiltonian flows to transform the prior distribution, in the sense of
Bayes’ theorem, $\pi_{0}$ of some vector of parameters $\mathbf{q}$ into the
target posterior distribution $\pi(\mathbf{q}|\boldsymbol{d})$ of these
parameters, knowing some data $\boldsymbol{d}$ and likelihood $\ell$. The main
difference with the above described NHF lies in the loss inspired from the KL
phase in Boltzmann Generators (Noé et al., 2019), as well as in the learning
procedure. During training, this NHF takes batches of $\mathbf{q}_{0}$ from
the prior distribution as inputs. For each $\mathbf{q}_{0}$, one
$\mathbf{p}_{0}$ is drawn from a Gaussian distribution whose mean and standard
deviation depend on the $\mathbf{q}_{0}$. The resulting point in phase-space
is then evolved through a series of $L$ Leapfrog steps with integration time
$dt$. The outputs consist in the final positions $\mathbf{q}_{T}$ and momenta
$\mathbf{p}_{T}$, as well as the initial mean $\mu(\mathbf{q}_{0})$, standard
deviation $\sigma(\mathbf{q}_{0})$ and $\mathbf{p}_{0}$. All these outputs, as
well as the data $\boldsymbol{d}$, are used in the loss computation. Once
trained, it is able to transform the prior into the desired posterior
distribution of the parameters. Thus, both training and sampling are now made
following the forward-direction flow from the prior to the posterior.

In order to compute the loss, one needs to have access to the likelihood
distribution $\ell$ of the model, which encapsulates the covariance matrix of
the data as well as the underlying physical mapping between vectors of
parameters and the corresponding data. Then, in the framework of Hamiltonian
dynamics, the full system is made of both positions (the parameters of
interest) and artificial momenta. Let us call $\mathbf{q}_{0},\mathbf{p}_{0}$
the initial position and momentum, respectively, and
$\mathbf{q}_{T},\mathbf{p}_{T}$ the corresponding final position and momentum,
respectively, obtained after $L$ Leapfrog transformations
$\mathcal{T}^{dt}_{1},...,\ \mathcal{T}^{dt}_{L}$ with timestep $dt$, i.e:
$(\mathbf{q}_{T},\mathbf{p}_{T})=\mathcal{T}^{dt}_{L}\circ...\circ\mathcal{T}^{dt}_{1}(\mathbf{q}_{0},\mathbf{p}_{0}):=\mathcal{T}(\mathbf{q}_{0},\mathbf{p}_{0})$.
Also, we introduce the notations for the projections along the final positions
and momenta, i.e.
$\mathbf{q}_{T}:=\mathcal{T}_{q}(\mathbf{q}_{0},\mathbf{p}_{0})$ and
$\mathbf{p}_{T}:=\mathcal{T}_{p}(\mathbf{q}_{0},\mathbf{p}_{0})$. By the
change of variables formula and the fact that the Jacobian determinant of each
Hamiltonian transformation is one, the model joint distribution $M$ is written
as:

$\displaystyle M(\mathbf{q}_{T},\mathbf{p}_{T})$
$\displaystyle=1\times\Pi_{0}(\mathcal{T}^{-dt}_{1}\circ...\circ\mathcal{T}^{-dt}_{L}(\mathbf{q}_{T},\mathbf{p}_{T}))$
$\displaystyle=\Pi_{0}(\mathbf{q}_{0},\mathbf{p}_{0})=\pi_{0}(\mathbf{q}_{0})\times
f(\mathbf{p}_{0}|\mathbf{q}_{0}),$

where $\Pi_{0}$ is the joint prior distribution, $\pi_{0}$ the prior
distribution of the parameters of interest and $f$ the Gaussian distribution
of the Encoder. We fix the target density of the final momenta $g(\mathbf{p})$
(e.g. Gaussian). We now seek to minimize the KL-divergence between the model
joint distribution and the desired target joint distribution conditioned on
data
$\Pi(\mathbf{q},\mathbf{p}|\boldsymbol{d})=\pi(\mathbf{q}|\boldsymbol{d})g(\mathbf{p})$.
We write the latter as the product of a density depending on $\mathbf{q}$ and
an other density depending on $\mathbf{p}$ so that the two are independent.
Using Bayes’ theorem, we have (see Appendix B for more details):

$\begin{split}&D_{KL}(M(\mathbf{q}_{T},\mathbf{p}_{T})\ ||\
\pi(\mathbf{q}_{T}|\boldsymbol{d})g(\mathbf{p}_{T}))=\\\
&\int\pi_{0}(\mathbf{q}_{0},\mathbf{p}_{0})\big{[}\log\pi_{0}(\mathbf{q}_{0})+\log
f(\mathbf{p}_{0}|\mathbf{q}_{0})-\log\pi_{0}(\mathcal{T}_{q}(\mathbf{q}_{0},\mathbf{p}_{0}))\\\
&\hskip
14.22636pt-\log\ell(\boldsymbol{d}|\mathcal{T}_{q}(\mathbf{q}_{0},\mathbf{p}_{0}))-\log
g(\mathcal{T}_{p}(\mathbf{q}_{0},\mathbf{p}_{0})\big{]}d\mathbf{q}_{0}d\mathbf{p}_{0}+\text{cst}.\end{split}$
(6)

###  6.2 Application to cosmology

We apply the herein-above architecture to cosmological analysis. One typical
case is the determination of the cosmic expansion, and more generally of the
cosmological parameters, from the observation of brightness and recession
velocity of Type Ia supernovæ (e.g. Riess et al., 1998; Betoule et al., 2014).
We use it as a toy example, however, while the model used so far has been
simple, it may be expanded in very complicated direction for which sampling
from the probability distribution becomes very complex. New observatory are
presently being built which is expected to deliver tens of thousands of new
supernovæ Ia over the next decade (LSST Science Collaboration, 2009). We now
present the essential feature of the model that we need before applying our
method to the toy example.

According to the $\Lambda$-CDM model, the relation between the distance and
the brightness of Type Ia supernovae is of great interest because it depends
on two cosmological parameters: the matter density parameter $\Omega_{m}$ and
the adimensional Hubble parameter $h$. To be more specific, database of type
Ia supernovae report the distance modulus $\mu$. This quantity is defined as
the difference between the apparent and the absolute magnitude of an
astronomical object, and is directly related to luminosity distance (Weinberg,
1972) and thus a function of the redshift $z$, $\Omega_{m}$ and $h$:

$\mu(z,\Omega_{m},h)=5\log_{10}\left(\frac{D_{L}^{*}(z,\Omega_{m})}{h10\text{pc}}\right)$

where
$D_{L}^{*}(z,\Omega_{m})=\frac{c(1+z)}{H_{0}}\int_{0}^{z}\frac{ds}{\sqrt{1-\Omega_{m}+\Omega_{m}(1+s)^{3}}},$
and $H_{0}=100~{}\text{km s}^{-1}~{}\text{Mpc}^{-1}$, $c$ being the speed of
light. In practice, we avoid computing the integral in $D^{*}_{L}$ by using an
approximation (Pen, 1999) which is only valid for a flat Universe:

$D_{L}^{*}(z,\Omega_{m})=\frac{c(1+z)}{H_{0}}\left[\eta(1,\Omega_{m})-\eta\left(\frac{1}{1+z},\Omega_{m}\right)\right],$

with

$\displaystyle\eta(a,\Omega_{m})$ $\displaystyle=$ $\displaystyle
2\sqrt{1+s^{3}}\left(\frac{1}{a}-0.1540\frac{s}{a^{3}}\right.$
$\displaystyle\left.+0.4304\frac{s^{2}}{a^{2}}+0.19097\frac{s^{3}}{a}+0.066941s^{4}\right).$

The formal definition of these quantities imposes constraints on the possible
values of the parameters, that can only be comprised between zero and one. We
avoid the problem by outputting a sigmoid of $\mathbf{q}_{T}$.

We aim to sample from the posterior distribution
$\pi(\Omega_{m},h|\text{data})$ quantifying the probability that we are living
in a universe whose mean density and expansion is equal to $\Omega_{m}$ and
$h$ given $D$ observations $\text{data}=\\{z_{i},\mu_{i}\\}_{1\leq i\leq D}$
of type Ia supernovae, and the covariance matrix $C$ of the observed distance
moduli. The likelihood of the problem is supposed to be Gaussian, i.e. the
observed data and the simulated output from parameters are different up to a
Gaussian noise. The final momenta distribution $g$ is set to a Normal
distribution. The cumulative plots in figure 8 compare the performance of a
fixed- and MLP-kinetic NHF with a HMC. They are both biased which is expected
since they minimize the KL-divergence between the model distribution and the
target, but less than with an ELBO. We leave for future work a possible
correction using importance sampling methods at the end of training.

Figure 8: Cumulative plots of means and standard deviations of both
$\Omega_{m}$ and $h$ produced by trained NHF models and a HMC on a
20,000-sample dataset, compared to the ground truth. Soft-uniform prior,
30,000 training epochs, $g$ chosen to be a $\mathcal{N}(0,I_{2})$

##  7 Conclusion

In this work, we analyzed how well Hamiltonian Normalizing Flows are suited to
solving some issues in Generative modeling. The main advantage of these
methods is twofold. First, the volume-preservation in phase-space avoids the
costly computation of Jacobian determinants. Then they allow for flexibility
in the type of neural networks that is used since reversibility is ensured by
the symplectic integrator. It is this flexibility that allows us to propose a
NHF variant based on classical kinetic energy. By exploring 2D sampling
problem, we showed that both variants are robust to the choice of
hyperparameters and prior distribution. In particular, we illustrate how the
explicit classical design of the kinetic energy is a way to increase such
robustness and make the model more easily interpretable. Finally, we explained
how to adapt NHF to the context of Bayesian inference to obtain a sampler of
the posterior distribution. Further work could then address how the offsets
could be corrected by importance sampling techniques on a trained model. More
generally, a next step would be testing such architectures on more complex
problems such as image generation to attest their behavior in a high-
dimensional setting.

## Acknowledgments

All the authors are grateful to the support from CNRS through the MITI-Prime
80 project ”CosmoBayes”. All the authors thank the Mésocentre Clermont
Auvergne University (https://mesocentre.uca.fr/) and the HPC cluster at
Fysikum (it.fysik.su.se/hpc/) for providing help, computing and storage
resources. This work was supported by the French ANR under the grant
ANR-20-CE46-0007 (SuSa project) and by the Simons Collaboration on “Learning
the Universe”. JJ acknowledges support by the Swedish Research Council (VR)
under the project 2020-05143 – “Deciphering the Dynamics of Cosmic Structure”.

We acknowledge the use of the following packages in this analysis: Numpy
(Harris et al., 2020), IPython (Pérez & Granger, 2007), Matplotlib (Hunter,
2007), Pytorch (Paszke et al., 2019), SciPy (Jones et al., 2001–).

## References

  * LSST Science Collaboration (2009) Abell, P. A., Allison, J., Anderson, S. F., Andrew, J. R., Angel, J. R. P., Armus, L., Arnett, D., Asztalos, S. J., et al.  LSST Science Book, Version 2.0.  _arXiv e-prints_ , art. arXiv:0912.0201, December 2009.
  * Ancona et al. (2022) Ancona, M., Ceolini, E., Öztireli, C., and Gross, M.  _Gradient-Based Attribution Methods_ , pp. 169–191.  Springer-Verlag, Berlin, Heidelberg, 2022.  ISBN 978-3-030-28953-9.  URL https://doi.org/10.1007/978-3-030-28954-6_9.
  * Arjovsky & Bottou (2017) Arjovsky, M. and Bottou, L.  Towards principled methods for training generative adversarial networks.  In _International Conference on Learning Representations_ , 2017.  URL https://openreview.net/forum?id=Hk4_qw5xe.
  * Berard et al. (2020) Berard, H., Gidel, G., Almahairi, A., Vincent, P., and Lacoste-Julien, S.  A closer look at the optimization landscapes of generative adversarial networks.  In _International Conference on Learning Representations_ , 2020.  URL https://openreview.net/forum?id=HJeVnCEKwH.
  * Betoule et al. (2014) Betoule, M., Kessler, R., Guy, J., Mosher, J., Hardin, D., Biswas, R., Astier, P., El-Hage, P., Konig, M., Kuhlmann, S., Marriner, J., et al.  Improved cosmological constraints from a joint analysis of the SDSS-II and SNLS supernova samples.  _Astronomy & Astrophysics_, 568:A22, aug 2014.  doi: 10.1051/0004-6361/201423413.  URL https://doi.org/10.1051%2F0004-6361%2F201423413.
  * Carpenter et al. (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A.  Stan: A probabilistic programming language.  _Journal of statistical software_ , 76(1), 2017.
  * Dhulipala et al. (2022) Dhulipala, S. L. N., Che, Y., and Shields, M. D.  Bayesian inference with latent hamiltonian neural networks, 2022.  URL https://arxiv.org/abs/2208.06120.
  * Dinh et al. (2014) Dinh, L., Krueger, D., and Bengio, Y.  Nice: Non-linear independent components estimation.  10 2014.
  * Duane et al. (1987) Duane, S., Kennedy, A., Pendleton, B. J., and Roweth, D.  Hybrid monte carlo.  _Physics Letters B_ , 195(2):216–222, 1987.  ISSN 0370-2693.  doi: https://doi.org/10.1016/0370-2693(87)91197-X.  URL https://www.sciencedirect.com/science/article/pii/037026938791197X.
  * Eghbal-zadeh et al. (2019) Eghbal-zadeh, H., Zellinger, W., and Widmer, G.  Mixture density generative adversarial networks.  _Proceedings / CVPR, IEEE Computer Society Conference on Computer Vision and Pattern Recognition. IEEE Computer Society Conference on Computer Vision and Pattern Recognition_ , 06 2019.  doi: 10.1109/CVPR.2019.00597.
  * Frazer et al. (2021) Frazer, J., Notin, P., Dias, M., Gomez, A., Min, J. K., Brock, K., Gal, Y., and Marks, D. S.  Disease variant prediction with deep generative models of evolutionary data.  _Nature_ , 599(7883):91–95, 2021.
  * Gilpin et al. (2018) Gilpin, L. H., Bau, D., Yuan, B. Z., Bajwa, A., Specter, M., and Kagal, L.  Explaining explanations: An overview of interpretability of machine learning.  In _2018 IEEE 5th International Conference on data science and advanced analytics (DSAA)_ , pp. 80–89. IEEE, 2018.
  * Goodfellow et al. (2014) Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y.  Generative adversarial nets.  In Ghahramani, Z., Welling, M., Cortes, C., Lawrence, N., and Weinberger, K. (eds.), _Advances in Neural Information Processing Systems_ , volume 27. Curran Associates, Inc., 2014.  URL https://proceedings.neurips.cc/paper/2014/file/5ca3e9b122f61f8f06494c97b1afccf3-Paper.pdf.
  * Greydanus et al. (2019) Greydanus, S., Dzamba, M., and Yosinski, J.  _Hamiltonian Neural Networks_.  Curran Associates Inc., Red Hook, NY, USA, 2019.
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del R’ıo, J. F., Wiebe, M., Peterson, P., G’erard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., and Oliphant, T. E.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, September 2020\.  doi: 10.1038/s41586-020-2649-2.  URL https://doi.org/10.1038/s41586-020-2649-2.
  * Homan & Gelman (2014) Homan, M. D. and Gelman, A.  The no-u-turn sampler: Adaptively setting path lengths in hamiltonian monte carlo.  _J. Mach. Learn. Res._ , 15(1):1593–1623, jan 2014.  ISSN 1532-4435.
  * Hunter (2007) Hunter, J. D.  Matplotlib: A 2d graphics environment.  _Computing In Science & Engineering_, 9(3):90–95, 2007.  doi: 10.1109/MCSE.2007.55.
  * Jimenez Rezende et al. (2019) Jimenez Rezende, D., Racanière, S., Higgins, I., and Toth, P.  Equivariant Hamiltonian Flows.  _arXiv e-prints_ , art. arXiv:1909.13739, September 2019.
  * Jones et al. (2001–) Jones, E., Oliphant, T., Peterson, P., et al.  SciPy: Open source scientific tools for Python, 2001–.  URL http://www.scipy.org/.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  In _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.  URL http://arxiv.org/abs/1412.6980.
  * Landau & Lifshitz (1982) Landau, L. and Lifshitz, E.  _Mechanics: Volume 1_.  Number vol. 1. Elsevier Science, 1982.  ISBN 9780080503479.  URL https://books.google.fr/books?id=bE-9tUH2J2wC.
  * Lapuschkin et al. (2019) Lapuschkin, S., Wäldchen, S., Binder, A., Montavon, G., Samek, W., and Müller, K.-R.  Unmasking clever hans predictors and assessing what machines really learn.  _Nature communications_ , 10(1):1096–1096, 2019\.  ISSN 2041-1723.
  * Lin et al. (2018) Lin, Z., Khetan, A., Fanti, G., and Oh, S.  Pacgan: The power of two samples in generative adversarial networks.  In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.  URL https://proceedings.neurips.cc/paper/2018/file/288cc0ff022877bd3df94bc9360b9c5d-Paper.pdf.
  * Lopez et al. (2020) Lopez, R., Gayoso, A., and Yosef, N.  Enhancing scientific discoveries in molecular biology with deep generative models.  _Molecular Systems Biology_ , 16(9):e9198, 2020\.
  * Menier et al. (2022) Menier, E., Bucci, M. A., Yagoubi, M., Mathelin, L., and Schoenauer, M.  Continuous Methods : Hamiltonian Domain Translation.  _arXiv e-prints_ , art. arXiv:2207.03843, July 2022.
  * Noé et al. (2019) Noé, F., Olsson, S., Köhler, J., and Wu, H.  Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning.  _Science_ , 365(6457):eaaw1147, 2019.  doi: 10.1126/science.aaw1147.  URL https://www.science.org/doi/abs/10.1126/science.aaw1147.
  * Papamakarios et al. (2022) Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S., and Lakshminarayanan, B.  Normalizing flows for probabilistic modeling and inference.  _J. Mach. Learn. Res._ , 22(1), jul 2022.  ISSN 1532-4435.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. Curran Associates, Inc., 2019.  URL http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * Pen (1999) Pen, U.-L.  Analytical fit to the luminosity distance for flat cosmologies with a cosmological constant.  _The Astrophysical Journal Supplement Series_ , 120(1):49–50, jan 1999.  doi: 10.1086/313167.  URL https://doi.org/10.1086%2F313167.
  * Pérez & Granger (2007) Pérez, F. and Granger, B. E.  IPython: a system for interactive scientific computing.  _Computing in Science and Engineering_ , 9(3):21–29, May 2007.  ISSN 1521-9615.  doi: 10.1109/MCSE.2007.53.  URL https://ipython.org.
  * Raissi et al. (2019) Raissi, M., Perdikaris, P., and Karniadakis, G. E.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational physics_ , 378:686–707, 2019\.
  * Rezende & Mohamed (2015) Rezende, D. and Mohamed, S.  Variational inference with normalizing flows.  In Bach, F. and Blei, D. (eds.), _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pp. 1530–1538, Lille, France, 07–09 Jul 2015\. PMLR.  URL https://proceedings.mlr.press/v37/rezende15.html.
  * Ribeiro et al. (2016) Ribeiro, M. T., Singh, S., and Guestrin, C.  ”why should i trust you?”: Explaining the predictions of any classifier.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pp. 1135–1144, New York, NY, USA, 2016. Association for Computing Machinery.  ISBN 9781450342322.  doi: 10.1145/2939672.2939778.  URL https://doi.org/10.1145/2939672.2939778.
  * Riess et al. (1998) Riess, A. G., Filippenko, A. V., Challis, P., Clocchiatti, A., Diercks, A., Garnavich, P. M., Gilliland, R. L., Hogan, C. J., Jha, S., Kirshner, R. P., Leibundgut, B., Phillips, M. M., Reiss, D., Schmidt, B. P., Schommer, R. A., Smith, R. C., Spyromilio, J., Stubbs, C., Suntzeff, N. B., and Tonry, J.  Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant.  _The Astronomical Journal_ , 116(3):1009–1038, September 1998.  doi: 10.1086/30049910.48550/arXiv.astro-ph/9805201.
  * Robert & Casella (2004) Robert, C. P. and Casella, G.  _Monte Carlo Statistical Methods_.  Springer, 2nd edition, 2004.  ISBN 0387212396.
  * Rodriguez et al. (2018) Rodriguez, A. C., Kacprzak, T., Lucchi, A., Amara, A., Sgier, R., Fluri, J., Hofmann, T., and Réfrégier, A.  Fast cosmic web simulations with generative adversarial networks.  _Computational Astrophysics and Cosmology_ , 5(1):1–11, 2018.
  * Samek & Müller (2019) Samek, W. and Müller, K.-R.  Towards explainable artificial intelligence.  In _Explainable AI: interpreting, explaining and visualizing deep learning_ , pp. 5–22. Springer, 2019.
  * Shen et al. (2023) Shen, J., Chowdhury, J., Banerjee, S., and Terejanu, G.  Machine fault classification using hamiltonian neural networks, 2023.  URL https://arxiv.org/abs/2301.02243.
  * Sohl-Dickstein et al. (2015) Sohl-Dickstein, J., Weiss, E., Maheswaranathan, N., and Ganguli, S.  Deep unsupervised learning using nonequilibrium thermodynamics.  In Bach, F. and Blei, D. (eds.), _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pp. 2256–2265, Lille, France, 07–09 Jul 2015\. PMLR.  URL https://proceedings.mlr.press/v37/sohl-dickstein15.html.
  * Tabak & Vanden-Eijnden (2010) Tabak, E. and Vanden-Eijnden, E.  Density estimation by dual ascent of the log-likelihood.  _Communications in Mathematical Sciences_ , 8(1):217–233, 2010.  ISSN 1539-6746.  doi: 10.4310/CMS.2010.v8.n1.a11.
  * Toth et al. (2020) Toth, P., Rezende, D. J., Jaegle, A., Racanière, S., Botev, A., and Higgins, I.  Hamiltonian generative networks.  In _International Conference on Learning Representations_ , 2020.  URL https://openreview.net/forum?id=HJenn6VFvB.
  * Weinberg (1972) Weinberg, S.  _Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity_.  1972\.
  * Winkler et al. (2019) Winkler, C., Worrall, D., Hoogeboom, E., and Welling, M.  Learning likelihoods with conditional normalizing flows, 2019.  URL https://arxiv.org/abs/1912.00042.

##  Appendix A Additional plots about robustness and optimization of NHF

We present the results of additional experiments with $H=32$ in figure 9.
Removing more than 90% of parameters (passage from $H=128$ to $H=32$), the
final values are always higher by less than 4%, for both models. The different
final values of the loss function can be represented on scatter plots, see
Figure 10. The latter clearly illustrate the robustness of the fixed-kinetic
model, see Figure 9. It also shows that models with $L=1$ perform poorer than
with $L=2,\ 10,\ 50$.

Figure 9: Training loss as a function of epochs for models with different $H$
(number of neurons per hidden layer in each neural network of the model), $L$
(number of Leapfrog steps) and $T$ (integration time). First row: MLP-kinetic
NHF; second row: fixed-kinetic model. From left to right: $H=128,\ 32,\ 8$.
Models were trained for $15,000$ epochs on a 5,000-point dataset with
minibatches of size $512$. Figure 10: Scatter plots $L$ vs. $T$ of final
values of the loss averaged on last 500 epochs. First row: MLP-kinetic NHF;
second row: fixed-kinetic model. From left to right: $H=128,\ 32,\ 8$. Models
were trained for $15,000$ epochs on a 5,000-point dataset with minibatches of
size $512$.

##  Appendix B Derivation of the KL-divergence and ELBO for the inference
problem

The KL-divergence suited to the inference problem is derived as follows:

$\begin{split}&D_{KL}(M(\mathbf{q}_{T},\mathbf{p}_{T})\ ||\
\pi(\mathbf{q}_{T}|\boldsymbol{d})g(\mathbf{p}_{T}))=\int
M(\mathbf{q}_{T},\mathbf{p}_{T})\log
M(\mathbf{q}_{T},\mathbf{p}_{T})d\mathbf{q}_{T}d\mathbf{p}_{T}-\int
M(\mathbf{q}_{T},\mathbf{p}_{T})\left[\log\pi(\mathbf{q}_{T}|\boldsymbol{d})+\log
g(\mathbf{p}_{T})\right]d\mathbf{q}_{T}d\mathbf{p}_{T}\\\
&=\int\Pi_{0}(\mathcal{T}^{-1}(\mathbf{q}_{T},\mathbf{p}_{T}))\log\Pi_{0}(\mathcal{T}^{-1}(\mathbf{q}_{T},\mathbf{p}_{T}))d\mathbf{q}_{T}d\mathbf{p}_{T}-\int
M(\mathbf{q}_{T},\mathbf{p}_{T})\left[\log\pi_{0}(\mathbf{q}_{T})+\log\ell(\boldsymbol{d}|\mathbf{q}_{T})-\log
p(\boldsymbol{d})+\log
g(\mathbf{p}_{T})\right]d\mathbf{q}_{T}d\mathbf{p}_{T}\\\
&=\int\Pi_{0}(\mathbf{q}_{0},\mathbf{p}_{0})\left[\log\pi_{0}(\mathbf{q}_{0})+\log
f(\mathbf{p}_{0}|\mathbf{q}_{0})\right]d\mathbf{q}_{0}d\mathbf{p}_{0}-\int
M(\mathbf{q}_{T},\mathbf{p}_{T})\left[\log\pi_{0}(\mathbf{q}_{T})+\log\ell(\boldsymbol{d}|\mathbf{q}_{T})+\log
g(\mathbf{p})\right]d\mathbf{q}_{T}d\mathbf{p}_{T}+\text{cst}\\\
&=\int\Pi_{0}(\mathbf{q}_{0},\mathbf{p}_{0})\left[\log\pi_{0}(\mathbf{q}_{0})+\log
f(\mathbf{p}_{0}|\mathbf{q}_{0})-\log\pi_{0}(\mathcal{T}_{q}(\mathbf{q}_{0},\mathbf{p}_{0}))-\log\ell(\boldsymbol{d}|\mathcal{T}_{q}(\mathbf{q}_{0},\mathbf{p}_{0}))-\log
g(\mathcal{T}_{p}(\mathbf{q}_{0},\mathbf{p}_{0})\right]d\mathbf{q}_{0}d\mathbf{p}_{0}+\text{cst}.\end{split}$
(7)

We can also adapt the ELBO from (Toth et al., 2020) to our inference
framework:

$\displaystyle\ln\pi_{0}(\mathbf{q}_{0})$
$\displaystyle=\ln\int\Pi_{0}(\mathbf{q}_{0},\mathbf{p}_{0})d\mathbf{p}_{0}$
$\displaystyle=\ln\int\frac{\Pi_{0}(\mathbf{q}_{0},\mathbf{p}_{0})}{f(\mathbf{p}_{0}|\mathbf{q}_{0})}f(\mathbf{p}_{0}|\mathbf{q}_{0})d\mathbf{p}_{0}$
$\displaystyle=\ln
E_{f}\left[\frac{\Pi_{0}(\mathbf{q}_{0},\mathbf{p}_{0})}{f(\mathbf{p}_{0}|\mathbf{q}_{0})}\right]$
$\displaystyle\geq E_{f}\left[\ln\Pi_{0}(\mathbf{q}_{0},\mathbf{p}_{0})-\ln
f(\mathbf{p}_{0}|\mathbf{q}_{0})\right]$ $\displaystyle=E_{f}\left[\ln
M(\mathcal{T}(\mathbf{q}_{0},\mathbf{p}_{0}))-\ln
f(\mathbf{p}_{0}|\mathbf{q}_{0})\right]$

Then, expliciting $M(q,p)=\pi_{0}(q)\ell(d|q)g(p)$:

$\displaystyle\text{ELBO}(\mathbf{q}_{0})$
$\displaystyle=E_{f}\left[\ln\left[\pi_{0}(T_{q}(\mathcal{T}(\mathbf{q}_{0},\mathbf{p}_{0})))\ell(d|T_{q}(\mathcal{T}(\mathbf{q}_{0},\mathbf{p}_{0})))g(T_{p}(\mathcal{T}(\mathbf{q}_{0},\mathbf{p}_{0})))\right]-\ln
f(\mathbf{p}_{0}|\mathbf{q}_{0})\right]$ (8)


#  Bäcklund Transformations of Kaup–Kupershmidt equations with Multi-soliton
solutions in Darboux Framework

Irfan Mahmood
Centre for High Energy Physics, University of the Punjab, 54590 Lahore,
Pakistan
Lahore
<EMAIL_ADDRESS>

###### Abstract

This article encloses the derivation of Darboux solutions for Kaup–Kupershmidt
equations with their generalization in determinantal form. One of the main
focuses of this work is to construct the Bäcklund transformation for the
different solutions of that equation through its associated Riccati equation
and then that transformations further reduces to its algebraic analogue with
the help of One-fold Darboux solution. Finally, its exact solutions upto three
solitons are calculated with their graphical representations which reveal
dynamical profiles of these solutions.

_K_ eywords Kaup–Kupershmidt equations, Darboux Transformation, Algebraic
Bäcklund Transformation, Solitons

##  1 Introduction

The nonlinear evolution equations which admit solitonic solutions have
achieved substantial attention in theory of integrable systems due to their
wide applications in various domains of physics and applied mathematics. The
importance to study such systems can be acknowledged by exploring the
applications of the One of most earliest paradigm equation, the Korteweg–De
Vries (KdV) equation, which plays a crucial role in the study of
hydrodynamics, describing the geometrical properties of wave propagations in
shallow water [1]. This equation with its higher order forms also applied
successfully as an integrable models in the analysis of electron plasma waves
associated with cylindrical plasma systems [2].

On a generalized fifth order KdV equationsIn this article, we discuss the
integrability of fifth order KdV like equations as the Kaup–Kupershmidt
equations [3] with their multi-soliton solutions. These equations are
efficiently applied as integrable models in control theory, to study nuclear
dynamics, in probe of fluid flow in elastic viscose medium [4, 5]. Moreover
the symmetric cases of these equations discussed in [6, 7] which consolidate
that equation from mathematical and physical poit of views. In literature, two
versions of Kaup–Kupershmidt equation exist which differ by coefficients, one
is called the Kaup–Kupershmidt type-I (KK-1) equation and its partner equation
is named as the Kaup–Kupershmidt type-II (KK-II) equation [8, 9]. These
computations include the investigation of the various integrable aspects of
only KK-I equation

$\displaystyle
u_{t}+4u^{3}u_{x}-\frac{75}{2}u_{x}u_{xx}-15uu_{xxx}+u_{xxxxx}=0$ (1)

rather then for the both versions because the second version of
Kaup–Kupershmidt equation

$\displaystyle v_{t}+20v^{3}v_{x}-25v_{x}v_{xx}+10vv_{xxx}+v_{xxxxx}=0$ (2)

is connected to KK-I equation through the scale transformation
$u=-\frac{2}{3}v$. These two equations also have been implemented successfully
in various physical systems as in nonlinear optics [10], the description of
the propagation of magneto-ion acoustic waves in plasma [11] and the analysis
of waves in oceans.

This work is pursued with several motivations one of these is that as KK-1
equation involves less mathematical complexity compared to higher oder KdV
systems, as Kaup [8] attempted to investigate its solitonic solutions but only
possibly obtained a solitary wave solution. Subsequently, its two- and three-
soliton solutions were derived in [9] using its Hirota form, without
presenting higher-order solitons. Another aspect which distinguishes the KK
equations from the KdV higher system is the absence of the $N$ generalization
of its multi-solitonic solutions in terms of Wronskian. The derivation of the
Hamiltonian structures of KK equations by C. T. Lee and C. C. Lee [8], with
their canonical forms and Lax representations, is acknowledged as a
substantial contribution towards exploring their more integrable aspects.

The main emphasis of this work is to investigate the multi-soliton solutions
of the KK-1 equation with their generalization in terms of Wronskian in
Darboux approach. The Darboux transformation [12, 13, 15] has been adopted as
an efficient integrable tool in soliton theory in calculating the exact
solutions, to construct the Bäcklund transformations deriving discrete
analogues of integrable systems which enhance its significance from a physical
point of view. Among the various remarkable applications of this method one of
these is its implementation to analyze the electrodynamical features [16]
associated with quantum cavity problems and also to elaborate the geometrical
properties of the graphene system [17] with its exact solitonic solutions. One
of the main focuses of this work is to derive the Bäcklund transformations of
KK-1 equation via its Riccati form in drivative forms then by make use one-
fold Darboux solution.

##  2 Linear representations and The Darboux solutions

This section encloses the Lax representation of KK-1 equation (1) and the
derivation of its one-fold, two-fold, and three-fold Darboux solutions by
using the Darboux transformation [12, 13] through the arbitrary function,
subsequently we generalize their $N$th-fold Darboux solutions in determinantal
form.

###  2.1 One Fold-Darboux solution

This can be manifested that KK-I equation (1) aries from the compatibility
condition of following linear system

$\displaystyle L\psi=\lambda\psi$ (3) $\displaystyle\psi_{t}=M\psi$ (4)

with Lax pair $L=\partial_{x}^{3}-3u\partial_{x}-\frac{3}{2}u_{x}$ and
$M=4\partial_{x}^{5}-45u\partial_{x}^{3}-\frac{135}{2}u_{x}\partial_{x}^{2}+\left(45u^{2}-\frac{105}{2}u_{x}x\right)\partial_{x}u+45uu_{x}-15u_{xxx}$
and $\lambda$ is a spectral parameter. The system (3) under the Draboux
transformation on arbitrary function $\psi$ as

$\displaystyle\psi\longrightarrow\tilde{\psi}=\psi_{x}-\sigma\psi$ (5)

generates one-fold transformation on $u$ in following form

$\displaystyle\tilde{u}=u-\sigma_{x}$ (6)

with inhomogeneous Riccati equation for $\sigma$

$\displaystyle\sigma^{2}+\sigma_{x}-2u=\lambda$ (7)

here $\sigma=\frac{\psi_{1x}}{\psi_{1}}$ and equation (6) connects the old
solution $u$ of KK-I equation with its new solution $u_{[1]}$ through the
$\sigma$ that satisfies the Riccati equation (7).

## Remark

The KK-II equations can be obtained from the compatibility of linear system
(3) and (4) with Lax pair $L=\partial_{x}^{3}+2v\partial_{x}+v_{x}$ and
$M=9\partial_{x}^{5}+30v\partial_{x}^{3}+45v_{x}\partial_{x}^{2}+\left(20v^{2}+35v_{x}x\right)\partial_{x}u+20vv_{x}-10v_{xxx}$.
Uner the Darboux transformation (5) we get $v_{[1]}=v+\frac{3}{2}\sigma_{x}$
that can also be obtained directly by replacing $u$ with $-\frac{2}{3}v$ in
KK-I one-fold Darbpux tranformation (6) therefore only integrability of KK-I
is studided in this this work.

##  3 Bäcklund Transformation

This section is devoted to construct the Bäcklund Transformation for KK-1
equation and then it is attempted to derive the algebraic analogue for
Bäcklund Transformation first order in space derivative. For this purpose let
start from Riccati equation (7) with one fold-Darboux tranformation. Now let
introduce $P_{x}=u$ in these equations and then combing the resulting
expression we get an equation relating $P$ and $\tilde{P}$ as below

$(\tilde{P}-P)^{2}-(\tilde{P}+P)_{x}=\lambda.$ (8)

Above last expression represents the parametric Bäcklund Transformation
connecting old solution $P_{x}=u$ with new solution $\tilde{P}=\tilde{u}$
which also can be acknowledged as nonlinear principle of superposition. In
order to construct the Bäcklund Transformation for the temporal part, first
write equation (1) for $P$ and $\tilde{P}$ then add the resulting expression
which involves $(\tilde{P}+P)_{xxxxx}$ that can be calculated in reduced order
by taking the fourth derivation of equation (8), the added expression with
this value will take the following form after simplification

$(\tilde{P}+P)_{t}=2(P-\tilde{P})(\tilde{P}_{xxxx}-P_{xxxx})+(15P-8P_{x})(\tilde{P}_{xxx}+P_{xxx})+(\tilde{P}_{x}+\tilde{P}_{xx})(4P_{xx}-6\tilde{P}_{x})-P_{x}Q_{1}+\tilde{P}_{x}Q_{2}+2\tilde{P}P_{xx}$
(9)

here $Q_{1}=6P_{x}+4P^{3}-\frac{75}{2}P_{xx}$ and
$Q_{2}=2P_{xxx}-4\tilde{P}^{3}+\frac{75}{2}\tilde{P}_{xx}$ and last result
represents the time derivative part of Bäcklund transformation. In following
proposition,

### Proposition3.1

In this proposition an algebraic analogue of Bäcklund Transformation (8) is
derived in following form

$(\hat{P}-\tilde{P})(\tilde{\hat{P}}+P)=\lambda-\eta$ (10)

Proof:
Let define transformation in two different directions with distinct
transformation parameters as $BT_{\lambda}:P\longrightarrow\tilde{P}$ and
$BT_{\eta}:P\longrightarrow\hat{P}$, the for parameter $\eta$, the Bäcklund
Transformation becomes

$(P-\hat{P})^{2}-(\hat{P}+P)_{x}=\eta.$ (11)

and now the above transormations (8) and (11) with two-fold shift in different
direction can be written as

$(\tilde{P}-\hat{\tilde{P}})^{2}-(\tilde{P}+\hat{\tilde{P}})_{x}=\eta$ (12)
$(\hat{P}-\tilde{\hat{P}})^{2}-(\hat{P}+\tilde{\hat{P}})_{x}=\lambda.$ (13)

Now after eliminating all the derivatives from one-fold and two-fold Bäcklund
transformations with commutative condition $\tilde{\hat{P}}=\hat{\tilde{P}}$
and after simplification we obtain expression (10) which does not involve any
derivative and connecting new solution solutions with initial solution $P$
through the transformation parameters $\lambda$ and $\eta$ that can be
regarded as the algebraic Bäcklund Transformation. We can also calculate
$\tilde{\hat{P}}$ direcly from that expression with initial solution $P$ and
$\tilde{P}$, $\hat{P}$ from Riccati equation.

###  3.1 Two and Three Fold Darboux transformations

In order to calculate the multi-fold Darboux tranformation let denote the
$\tilde{u}=u[1]$ as for one-fold transformation and so on. The two-fold
Darboux transformation for arbitrary function $\psi$ can be expressed as below

$\psi_{[2]}=\left(\frac{d}{dx}-a_{2}\right)\left(\frac{d}{dx}-a_{1}\right)\psi,$
(14)

here $a_{2}=\frac{\psi_{2x}[1}{\psi_{2}[1]}$,
$a_{1}=\frac{\psi_{1x}}{\psi_{1}}$ and $\psi_{2}[1]$ can be calculated form
following expression

$\psi_{2}[1]=\psi_{2x}-a_{1}\psi_{2},$ (15)

where $\psi_{2}$ is the particular solution of linear systems (3) and (4) at
$\lambda=\lambda_{2}$, simply we can write two-fold transformation (14) as
ratio of Wronskians

$\psi_{[2]}=\frac{\textit{Wron}(\psi_{1},\psi_{2},\psi)}{\textit{Wron}(\psi_{1},\psi_{2})},$
(16)

with

$\displaystyle\textit{Wron}(\psi_{1},\psi_{2},\psi)$ $\displaystyle=$
$\displaystyle\begin{vmatrix}\psi_{1}&\psi_{2}&\psi\\\
\psi_{1x}&\psi_{2x}&\psi_{x}\\\ \psi_{1xx}&\psi_{2xx}&\psi_{xx}\\\
\end{vmatrix},$ $\displaystyle\textit{Wron}(\psi_{1},\psi_{2})$
$\displaystyle=$ $\displaystyle\begin{vmatrix}\psi_{1}&\psi_{2}\\\
\psi_{1x}&\psi_{2x}\\\ \end{vmatrix},$ (17)

here $\psi_{i}$ are the particular solutions of linear systems at
$\lambda=\lambda_{i}$. Similarly, we can construct the three-fold
transformation for arbitrary function $\psi$ as
$\psi[3]=\left(\frac{d}{dx}-\frac{\psi_{3x}[2]}{\psi_{3}[2]}\right)\psi[2]$ in
terms of the Wronskians as below

$\psi[3]=\frac{\textit{Wron}(\psi_{1},\psi_{2},\psi_{3},\psi)}{\textit{Wron}(\psi_{1},\psi_{2},\psi_{3})},$
(18)

where $\textit{Wron}(\psi_{1},\psi_{2},\psi_{3},\psi)$ and
$\textit{Wron}(\psi_{1},\psi_{2},\psi_{3})$ are the Wronskian of order $4$ and
order $3$ respectively. Upto here one-fold, two-fold and three-fold Darboux
transformations for arbitrary function $\psi$ are derived but following
proposition $2.2$, we elaborate the procedure [12, 13, 14] to generalize the
$N$th-fold transformation for $u$ in compact form as the logarithmic
derivative of $N$-th order Wronskians.

### Proposition3.2

This can be shown that by using one-fold DT, we may construct $N$th-fold
Darboux transformations for the field variable $u$ in the following compact
form

$u[N]=u-h_{1x},$ (19)

here
$h_{1}=\frac{d}{dx}\log{\textit{W}\left(\psi_{1},\psi_{2},...,\psi_{N}\right)}$
and $\textit{W}(\psi_{1},\psi_{2},......,\psi_{N})$ is Wronskian of order $N$.
Proof:
As the second step Darboux transformation on $\psi$ can be expressed as below

$\psi[2]=\left(\frac{d}{dx}-a_{2}\right)\psi[1]=\left(\frac{d}{dx}-a_{2}\right)\left(\frac{d}{dx}-a_{1}\right)\psi,$
(20)

and similarly the two-fold Darboux transformation on $u$ immediately

$u[2]=u+a_{1x}+a_{2x}.$ (21)

Now after $N$ iteration we obtain the $N$th-fold Darboux transformations as
below

$\psi[N]=\left(\frac{d}{dx}-a_{N}\right)\left(\frac{d}{dx}-a_{N-1}\right)\left(\frac{d}{dx}-a_{N-2}\right)...\left(\frac{d}{dx}-a_{2}\right)\left(\frac{d}{dx}-a_{1}\right)\psi.$
(22) $u[N]=u-\sum_{i=1}^{N}a_{ix}.$ (23)

with $a_{i}=\psi_{ix}[i-1]\psi^{-1}_{i}[i-1]$. Now in order to construct
$N$th-fold Darboux transformations in terms of Wronskian, we start with $N$th-
fold Darboux transformation for arbitrary function $\psi$ in following form

$\psi[N]=D[N]\psi=\psi^{(N)}+h_{1}\psi^{(N-1)}+h_{2}\psi^{(N-2)}+.....+h_{N-1}\psi^{(1)}+h_{N}\psi,$
(24)

which is an equivalent representation of (22), now we can easily construct
$N$th-fold transformation for $u$ in following form

$u[N]=u-h_{1x}$ (25)

and $\eta_{1}$ may be determine from the following system of $N$ linear
algebraic equations

$\sum_{i=1}^{N}\eta_{i}\psi^{N-i}=-\psi^{N}_{i}.$ (26)

by applying the Kramer rule and explicitly can be written as the logarithmic
derivation of Wronskian as below

$h_{1}=\frac{d}{dx}\log{Wro\left(\psi_{1},\psi_{2},...,\psi_{N}\right)}.$ (27)

Now we can directly write the $N$th-fold expression for arbitrary function
$\psi$ as the ratio of Wronskians in following form

$\psi[N]=\frac{Wron(\psi_{1},\psi_{2},......,\psi_{N},\psi)}{Wron(\psi_{1},\psi_{2},......,\psi_{N})},$
(28)

here

$Wron(\psi_{1},\psi_{2},....\psi_{N})=\begin{vmatrix}\psi_{1}&\psi_{2}&\cdots&\psi_{N}\\\
\psi_{1}^{(1)}&\psi_{2}^{(1)}&\cdots&\psi_{N}^{(1)}\\\
\vdots&\vdots&\ddots&\vdots\\\
\psi_{1}^{(n-1)}&\psi_{2}^{(n-1)}&\cdots&\psi_{N}^{(n-1)}\end{vmatrix},$ (29)

and

$Wron(\psi_{1},\psi_{2},....\psi_{N},\psi)=\begin{vmatrix}\psi_{1}&\psi_{2}&\cdots&\psi_{N}&\psi\\\
\psi_{1}^{(1)}&\psi_{2}^{(1)}&\cdots&\psi_{N}^{(1)}&\psi^{(1)}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
\psi_{1}^{(n-1)}&\psi_{2}^{(n-1)}&\cdots&\psi_{N}^{(n-1)}&\psi^{(n-1)}\end{vmatrix}$
(30)

in above determinants $\psi_{j}^{(i)}$ stands for $i$-th derivative of
$\psi_{j}$ with respect to $x$ as
$\psi_{j}^{(i)}=\frac{d^{i}\psi_{j}}{dx^{i}}$.

##  4 Exact soliton solutions

This section encloses the derivation of exact soliton solutions of KK-I
equation (1) upto three-soliton solution through its $N$-fold Darboux
transformation taking intial solution $u=0$ with their graphical
representations that reveal their dynamical profiles beforle and after
interactions.

###  4.1 One soliton solutions

Staring with simplest trivial solution of equation (1) as $u=0$ the particular
solution $\psi_{1}$ at $\lambda=\lambda_{1}$ from the linear equations (3) and
(4) can be calculated as below

$\displaystyle\psi_{1}=2\cosh\left(4\text{$\lambda_{1}$}^{5/3}t+\sqrt[3]{\text{$\lambda_{1}$}}x\right).$
(31)

Now substituting the value of $\psi_{1}$ into one-fold transformation (6) with
$u=0$ and the after simplification we get

$\displaystyle
u[1]=-2\text{$\lambda_{1}$}^{2/3}\text{sech}^{2}\left(4\text{$\lambda_{1}$}^{5/3}t+\sqrt[3]{\text{$\lambda_{1}$}}x\right)$
(32)

is the one soliton solution of equation (1) and its graphical representations
in different frames are manifested as below

(a) 1-D

(b) 3-D

Figure 1: (a) represents the one dimensional dynamics of one soliton solution
where as in (b) its three dimensional profile has been shown.

###  4.2 Two soliton solutions

The two soliton solution can be calculated by taking $N=2$ in its $N$-fold
Darboux transformations or from second iteration of its one-fold Darboux
transformation as

$u[2]=u[1]-\psi_{2}^{\prime}[1]\psi_{2}^{-1}[1].$ (33)

and the value for $\psi_{2}[1]$ can be obtained from

$\psi_{2}[1]=\left(\frac{d}{dx}-\sigma_{1}\right)\psi_{2}.$ (34)

here $\psi_{2}$ is the particular solution of system (3) and (3) and
explicitly can be written as

$\displaystyle\psi_{2}=2\sinh\left(4\text{$\lambda$2}^{5/3}t+\sqrt[3]{\text{$\lambda$2}}x\right).$
(35)

Now the value of $\psi_{2}[1]$ can be obtained from expression (34) with the
help of above calculated functions $\sigma_{1}$ and $\psi_{2}$ in subsequent
form and the two soliton solution is calculated as follows

$u[2]=\frac{2\left(\text{$\lambda_{1}$}^{2/3}-\text{$\lambda_{2}$}^{2/3}\right)\left(-\text{$\lambda_{1}$}^{2/3}+\text{$\lambda_{2}$}^{2/3}+\text{$\lambda_{1}$}^{2/3}\cosh\left(2\theta_{2}\right)+\text{$\lambda_{2}$}^{2/3}\cosh\left(2\theta_{1}\right)\right)}{\left(\sqrt[3]{\text{$\lambda_{2}$}}\cosh\left(\theta_{1}\right)\cosh\left(\theta_{2}\right)-\sqrt[3]{\text{$\lambda_{1}$}}\sinh\left(\theta_{1}\right)\sinh\left(\theta_{2}\right)\right)^{2}}$
(36)

where
$\theta_{1}=\sqrt[3]{\text{$\lambda_{1}$}}\left(4\text{$\lambda_{1}$}^{4/3}t+x\right)$
and
$\theta_{2}=\sqrt[3]{\text{$\lambda_{2}$}}\left(4\text{$\lambda_{2}$}^{4/3}t+x\right)$

(a) 1-D

(b) 1-D

Figure 2: (a) represents the one dimensional dynamics of two soliton solution
before interaction where as in (b) one dimensional dynamics at the time of
interaction has been shown.

(a) 1-D

(b) 3-D

Figure 3: (a) represents the one dimensional dynamics of two soliton solution
after interaction where as in (b) the three dimensional dynamics of two
soliton solution has been shown

###  4.3 Three soliton solution

The one step forward iteration on two fold Darboux transformation (33) yields
three soliton solution as below

$u[3]=u[2]-\psi_{3}^{\prime}[2]\psi_{3}^{-1}[2].$ (37)

and the value for $\psi_{3}[2]$ can be obtained from

$\psi_{3}[2]=\left(\frac{d}{dx}-\sigma_{2}\right)\psi_{3}.$ (38)

here $\psi_{3}$ is the particular solution of system (3) and (3)
at$\sigma_{1}$, explicitly can be written as

$\displaystyle\psi_{3}=2\cosh\left(4\text{$\lambda_{3}$}^{5/3}t+\sqrt[3]{\text{$\lambda_{3}$}}x\right)$
(39)

Now the value of $\psi_{2}[1]$ can be obtained from expression (34) with the
help of above calculated functions $\sigma_{1}$ and $\psi_{2}$ and after
putting all the values in Eq.(37) ,we get the three soliton solution as follow

$\displaystyle u[3]$ $\displaystyle=$
$\displaystyle-\Bigg{[}\Bigg{(}-(\lambda_{1}^{1/3}(-\lambda_{1}^{2/3}+\lambda_{3}^{2/3})\lambda_{3}^{1/3}\sinh{(\theta_{1})}\sinh{(\theta_{2})}\sinh{(\theta_{3})}+\lambda_{2}^{1/3}\cosh{(\theta_{2})}(\lambda_{1}^{1/3}(\lambda_{1}^{2/3}$
(40) $\displaystyle-$
$\displaystyle\lambda_{2}^{2/3})\cosh{\theta_{3}}\sinh{\theta_{1}}+(\lambda_{2}^{2/3}-\lambda_{3}^{2/3})\lambda_{3}^{1/3}\cosh{(\theta_{1})}\sinh{(\theta_{3})}))^{2}+(\lambda_{1}^{1/3})(-\lambda_{2}^{2/3}$
(41) $\displaystyle+$
$\displaystyle\lambda_{3}^{2/3})\cosh{(\theta_{3})}\sinh{(\theta_{1})}\sinh{(\theta_{2})}+\cosh{(\theta_{1})}(\lambda_{2}^{1/3}(\lambda_{1}^{2/3}-\lambda_{3}^{2/3})\cosh{(\theta_{2})}\cosh{(\theta_{3})}$
(42) $\displaystyle+$
$\displaystyle(-\lambda_{1}^{2/3}+\lambda_{2}^{2/3})\lambda_{3}^{1/3}\sinh{(\theta_{2})}\sinh{(\theta_{3})}))(\lambda_{1}^{1/3}(\lambda_{1}^{2/3}-\lambda_{2}^{2/3}-\lambda_{3}^{2/3})(\lambda_{2}^{2/3}$
(43) $\displaystyle-$
$\displaystyle\lambda_{3}^{2/3})\cosh{(\theta_{3})}\sinh{(\theta_{1})}\sinh{(\theta_{2})}+\cosh{(\theta_{1})}(-\lambda_{2}^{1/3}(\lambda_{1}^{2/3}-\lambda_{3}^{2/3})(-\lambda_{1}^{2/3}+\lambda_{2}^{2/3}$
(44) $\displaystyle-$
$\displaystyle\lambda_{3}^{2/3})\cosh{(\theta_{2})}\cosh{(\theta_{3})}-(\lambda_{1}^{2/3}-\lambda_{2}^{2/3})(\lambda_{1}^{2/3}+\lambda_{2}^{2/3}-\lambda_{3}^{2/3})\lambda_{3}^{1/3}\sinh{(\theta_{2})}\sinh{\theta_{3}}))\Bigg{)}$
(45) $\displaystyle/$
$\displaystyle\Bigg{(}\lambda_{1}^{1/3}(-\lambda_{2}^{2/3}+\lambda_{3}^{2/3})\cosh{(\theta_{3})}\sinh{(\theta_{1})}\sinh{(\theta_{2})}+\cosh{(\theta_{1})}(\lambda_{2}^{1/3}(\lambda_{1}^{2/3}$
(46) $\displaystyle-$
$\displaystyle\lambda_{3}^{2/3})\cosh{(\theta_{2})}\cosh{(\theta_{3})}+(-\lambda_{1}^{2/3}+\lambda_{2}^{2/3})\lambda_{3}^{1/3}\sinh{(\theta_{2})}\sinh{(\theta_{3})})\Bigg{)}^{2}\Bigg{]}$
(47)

where
$\theta_{1}=\sqrt[3]{\text{$\lambda_{1}$}}\left(4\text{$\lambda_{1}$}^{4/3}t+x\right)$,
$\theta_{2}=\sqrt[3]{\text{$\lambda_{2}$}}\left(4\text{$\lambda_{2}$}^{4/3}t+x\right)$
and
$\theta_{3}=\sqrt[3]{\text{$\lambda_{3}$}}\left(4\text{$\lambda_{3}$}^{4/3}t+x\right)$

(a) 1-D

(b) 1-D

Figure 4: (a) represents the one dimensional dynamics of three soliton
solution before interaction where as in (b) one dimensional dynamics at the
time of interaction has been shown.

(a) 1-D

(b) 3-D

Figure 5: (a) represents the one dimensional dynamics of three soliton
solution after interaction where as in (b) the three dimensional view of three
soliton solution has been shown (a) 1-D Figure 6: (a) represents the contour
plot of three soliton solution

## Conclusion

This work enclosed the derivation of exact solutions of Kaup–Kupershmidt
equation with its Bäcklund Transformation in aglebraic form through Draboux
transformations. Further motivation is that through the Draboux transfirmation
presented in this paper one may construct its discrete and quadrilaterla form.
More interestingly these computations may be extended to noncommutative spaces
to study the $N$-fold Darboux solutions in quasideterminat form.

## Acknowledgement

I am very thankful to the Punjab University 54590 on providing me facilities
to complete this research project

## References

  * [1] N. Zabusky and C. Galvin. Shallow-water waves, the Korteweg-deVries equation and solitons. Journal of Fluid Mechanics, 47(4):811–824, 1971.
  * [2] A. Khater, D. Callebaut, and A. Seadawy. General soliton solutions for nonlinear dispersive waves in convective type instabilities. Physica Scripta, 74(3):384, 2006.
  * [3] D. Kaya and S. M. El-Sayed. On a generalized fifth order KdV equations. physics Letters A, 310(1):44–51, 2003.
  * [4] J. J. T. A. A. Kilbas, H. M. Srivastava. Theory and Applications of Fractional Differential Equations. Elsevier, New York, 204, 2006.
  * [5] O. M. S.G. Samko, A.A. Kilbas. Fractional Integrals and Derivatives: Theory and Applications. 2002.
  * [6] E. G. Reyes. Nonlocal symmetries and the Kaup–Kupershmidt equation. Journal of mathematical physics, 48, 2005\.
  * [7] G. S. EG Reyes. Explicit solutions to the Kaup–Kupershmidt equation via nonlocal symmetries. International Journal of Bifurcation and Chaos, 08, 2007.
  * [8] C. Lee and C. Lee. Lax pairs and Hamiltonians for the Kaup–Kupershmidt-type equation. Physica Scripta, 85(3):035004, 2012.
  * [9] W. Hereman and A. Nuseir. Symbolic methods to construct exact solutions of nonlinear partial differential equations. Mathematics and Computers in Simulation, 43(1):13–27, 1997.
  * [10] M. Gedalin and T. Scott. Band YB. J. Phys. Rev. Lett., 1997:78, 1997.
  * [11] G. B. Whitham. Linear and nonlinear waves. John Wiley $\And$ Sons, 2011.
  * [12] V. B. Matveev and M. A. Salle. Darboux transformations and solitons. Springer series in nonlinear dynamics, 1991\.
  * [13] C. Gu, H. Hu, and Z. Zhou. Darboux transformations in integrable systems: theory and their applications to geometry. Springer Science $\And$ Business Media, 2004.
  * [14] Q. H. Park and H. J. Shin. Darboux transformation and Crum’s formula for multi-component integrable equations. Phys. D Nonlinear Phenom, 2001.
  * [15] A. Trisetyarso. Application of Darboux Transformation to solve Multisoliton Solution on Non-linear Schrödinger Equation. arXiv preprint arXiv:0910.0901, 2009.
  * [16] A. Trisetyarso. Correlation of Dirac potentials and atomic inversion in cavity quantum electrodynamics. Journal of Mathematical Physics, 51(7), 2010.
  * [17] A. Trisetyarso. Dirac four-potential tunings-based quantum transistor utilizing the Lorentz force. arXiv preprint arXiv:1003.4590, 2010.