# Average dissipation for stochastic transport equations with Lévy noise

Franco Flandoli 111Faculty of Sciences, Scuola Normale Superiore Pisa, Italy.
E-mail<EMAIL_ADDRESS>Andrea Papini 222Department of Mathematical
Sciences, Chalmers University of Technology & University of Gothenburg,
Gothenburg, Sweden. E-mail<EMAIL_ADDRESS>Marco Rehmeier333Faculty of
Sciences, Scuola Normale Superiore Pisa, Italy 444Faculty of Mathematics,
Bielefeld University, Bielefeld, Germany. E-mail<EMAIL_ADDRESS>bielefeld.de

###### Abstract

We show that, in one spatial and arbitrary jump dimension, the averaged
solution of a Marcus-type SPDE with pure jump Lévy transport noise satisfies a
dissipative deterministic equation involving a fractional Laplace-type
operator. To this end, we identify the correct associated Lévy measure for the
driving noise. We consider this a first step in the direction of a non-local
version of enhanced dissipation, a phenomenon recently proven to occur for
Brownian transport noise and the associated local parabolic PDE by the first
author. Moreover, we present numerical simulations, supporting the fact that
dissipation occurs for the averaged solution, with a behavior akin to the
diffusion due to a fractional Laplacian, but not in a pathwise sense.

Keywords: Marcus canonical integral; stochastic PDEs; Lévy process; transport
equation; dissipation
MSC2020 subject classifications: 60H15; 76F55; 76F25; 35Q49; 35K10; 60J76

##  1 Introduction

Stochastic transport (and advection, not discussed in this paper) attracts
more and more attention for its potentialities to describe small scale
turbulence in several models and applications, see for instance the volumes
[8, 9, 10] or the application to raindrop formation [7] and turbulence in
pipes [2], among many others. Small scale turbulence is described by
stochastic processes, space-dependent, either given a priori or inferred from
data. Most models deal with white noise or Ornstein–Uhlenbeck processes, the
first basic paradigms for any investigation of this kind. However, turbulent
signals may be more complex. Two classes of processes seem to be the first
ones to be considered after Gaussian noise: Fractional Gaussian noise, and
$\alpha$-stable processes. The first one has been considered in [11],
reporting some preliminary results. In this paper we give some preliminary
results on $\alpha$-stable processes, which seem to be the first of its kind.
We also include a list of questions for future work in this direction.

The property of turbulent fluids we want to emphasize is the additional
dissipation produced by turbulent eddies. In the white noise case this has
been widely investigated, see for instance [9] and the references therein.
Here, we consider this property in the case of $\alpha$-stable transport
noise. The most interesting case is undoubtedly when this transport noise
models the small scales acting on the large ones in nonlinear models, but such
a case is still too difficult in the $\alpha$-stable case. We limit ourselves
to linear transport of a scalar quantity, for instance heat, by an
$\alpha$-stable noise. Even in this simplified setting, many difficult
technical questions emerge. The first one is which notion of integral should
be used. In the Gaussian white noise case the basic rule is to choose a noise
providing the correct invariance (conservation) properties. This leads to
Stratonovich noise. In the $\alpha$-stable case the same invariance is given
by the so-called Markus noise. This type of stochastic integral has already
been developed in the literature, both in the finite- and infinite-dimensional
case [13, 3, 12, 1]. Its main advantage is that it preserves the ordinary
rules of calculus for general Lévy processes as integrators. In this sense,
the Marcus integral can be considered a natural extension of the Stratonovich
integral. Indeed, for diffusion integrators, both integral notions coincide.
We briefly review its definition in Section 2.

Given the model, the noise and the meaning of stochastic calculus, we limit
ourselves here to investigate the following question: whether the expected
value of the solution is dissipated, and whether this expected value satisfies
a closed equation, similarly to the Gaussian noise case. We give affirmative
answers to these questions: the expected value of the solution to the Marcus-
type SPDE (2.1), for a suitable choice of an essentially $\alpha$-stable
symmetric Lévy measure, satisfies the second-order parabolic deterministic
equation (3.4), see Proposition 3.2. The precise shape of the operator in
(3.4) depends, of course, on the chosen Lévy measure for the pure jump Lévy
process $Z$ from (2.1). With the aforementioned choice of an $\alpha$-stable
symmetric Lévy measure, the resulting operator is close (but not identical) to
the $\alpha$-fractional Laplacian.

In order to support our claims, we include numerical simulations, which
underline the dissipative character of the expected value of the solution to
(2.1). We limit ourselves to dimension $d=1$, in which an explicit solution
can be computed via the method of characteristics. We simulate the evolution
of the averaged solution profile, obtained with Monte Carlo method, showing a
decay in both time and space, as expected from our theoretical results. More
so, starting with a compact supported initial condition, we analyze, in
Figures 2(b), 3(a) and 3(b), the decay in time in the origin $x=0$, obtaining
an asymptotic behavior $\sim\beta(\sigma,m,\alpha)t^{-1/\alpha}$, with $\beta$
depending on the velocity field norm, the dimension of the Lévy process, and
the parameter $\alpha$ of the Lévy measure $\nu$. This decay corresponds to
the one of nonlocal PDEs. We remark that the pathwise solution profile shows
no dissipativity when $\sigma$ is constant and $d=1$. It will be interesting
to devote future work to the question whether pathwise dissipative behavior
can be observed as well, for instance by improving the mixing property of the
velocity field.

The organization of this paper is as follows. In Section 2, we introduce our
model and recall the notion and basic properties of Marcus stochastic integral
equations. In Section 3, we state our theoretical results, see in particular
Proposition 3.2. We present and discuss numerical simulations in Section 4
and, finally, pose some open questions for future research in Section 5.

##  2 Stochastic transport equations of Marcus-type

We consider the following transport Marcus-type SPDE on
$\mathbb{R}_{+}\times\mathbb{R}^{d}$

$du(t,x)=(\sigma(x)\nabla u(t,x))\diamond dZ_{t},\quad u(0,x)=u_{0}(x),$ (2.1)

where $\sigma:\mathbb{R}^{d}\to\mathbb{R}^{d\times m}$ and
$u_{0}:\mathbb{R}^{d}\to\mathbb{R}$. $Z$ is an $m$-dimensional pure jump Lévy
process on a filtered probability space
$(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0},\mathbb{P})$,

$Z_{t}=\int_{0}^{t}\int_{\overline{B_{1}(0)}}z\tilde{N}(dz,ds)+\int_{0}^{t}\int_{\overline{B_{1}(0)}^{c}}zN(dz,ds),\quad
t\geq 0,$

with Poisson random measure $N$, Lévy measure $\nu$ (i.e. $\nu$ is a Borel
probability measure on $\mathbb{R}^{m}$ with $\nu(\\{0\\})=0$ and
$\int_{\mathbb{R}^{m}}\min(1,z^{2})d\nu(z)<\infty$), and
$\tilde{N}(dz,dt)=N(dz,dt)-\mathds{1}_{\overline{B_{1}(0)}}\nu(dz)dt$. Here
$B_{1}(0)$ denotes the Euclidean ball in $\mathbb{R}^{m}$ with radius $1$
centered at $0$, and $\overline{B_{1}(0)}$ its closure. We make specific
choices for $\sigma$ and $\nu$ below. The symbol $\diamond$ denotes the Marcus
stochastic integral, i.e. (2.1) is understood in the following integral sense:

$\displaystyle u(t,x)=u_{0}(x)$
$\displaystyle+\int_{0}^{t}\int_{\overline{B_{1}(0)}}e^{\sigma
z}u(s-,x)-u(s-,x)\,\tilde{N}(dz,ds)$ (2.2)
$\displaystyle+\int_{0}^{t}\int_{\mathbb{R}^{m}\backslash\overline{B_{1}(0)}}e^{\sigma
z}u(s-,x)-u(s-,x)\,N(dz,ds)$
$\displaystyle+\int_{0}^{t}\int_{\overline{B_{1}(0)}}e^{\sigma
z}u(s-,x)-u(s-,x)-\nabla
u(s-,x)\cdot(\sigma(x)z)\,d\nu(z)ds,\quad(t,x)\in\mathbb{R}_{+}\times\mathbb{R}^{d},$

where for $z\in\mathbb{R}^{m}$ and $f:\mathbb{R}^{d}\to\mathbb{R}$, $e^{\sigma
z}f$ denotes the solution $g$ of

$\displaystyle\partial_{t}g(t,x)=\nabla g(t,x)\cdot(\sigma(x)z),\quad
g(0,x)=f(x),$ (2.3)

evaluated at $t=1$. Here $s-$ denotes the left limit of $s\in\mathbb{R}$. When
$\sigma$ and $f$ are sufficiently regular, the solution to this first-order
linear transport PDE is unique and given by $g(t,x)=f(\phi_{t,0}(x))$, where
$(t,x)\mapsto\phi_{t,0}(x)$ is the inverse of the unique solution flow
$(t,x)\mapsto\phi_{0,t}(x)$ for the ODE

$\partial_{t}\phi_{0,t}(x)=-\sigma(\phi_{0,t}(x))z,\quad\phi_{0,0}(x)=x$

on $\mathbb{R}\times\mathbb{R}^{d}$. We recall the following definition and
result from [12].

######  Definition 2.1.

An $(\mathcal{F}_{t})$-adapted random field
$u:\mathbb{R}_{+}\times\mathbb{R}^{d}\times\Omega\to\mathbb{R}$ is a solution
to (2.1), if it is a càdlàg $C^{2}$-semimartingale and (2.2) is satisfied for
$\mathbb{P}$-a.e. $\omega\in\Omega$.

######  Proposition 2.2.

If $u_{0}\in C^{2}_{b}(\mathbb{R}^{d})$ and $\sigma\in
C^{4}_{b}(\mathbb{R}^{d},\mathbb{R}^{d\times m})$, then there is a unique
solution to (2.1), and it is given by

$u(t,x)=u_{0}(\varphi_{t,0}(x)),$

where $(t,x)\mapsto\varphi_{t,0}(x)$ denotes the inverse of the stochastic
flow of the Marcus-SDE

$\varphi_{0,t}(x)=x-\int_{0}^{t}\sigma(\varphi_{0,s-}(x))\diamond dZ_{t},\quad
t\geq 0,x\in\mathbb{R}^{d}.$

The definition of solution to this SDE is similar to the infinite-dimensional
case, precisely it is given by

$\displaystyle\varphi_{0,t}(x)=x$
$\displaystyle-\int_{0}^{t}\int_{\overline{B_{1}(0)}}\Psi_{z}(\varphi_{0,s-}(x))-\varphi_{0,s-}(x)\,d\tilde{N}(dz,ds)$
$\displaystyle-\int_{0}^{t}\int_{\mathbb{R}^{m}\backslash\overline{B_{1}(0)}}\Psi_{z}(\varphi_{0,s-}(x))-\varphi_{0,s-}(x)\,dN(dz,ds)$
$\displaystyle-\int_{0}^{t}\int_{\overline{B_{1}(0)}}\Psi_{z}(\varphi_{0,s-}(x))-\varphi_{0,s-}(x)-z\sigma(\varphi_{0,s-}(x))\,\nu(z)ds,$

where $\Psi_{z}(y)$ denotes the solution to

$\partial_{t}f(t)=z\sigma(f(t)),\quad t\in\mathbb{R},\quad\quad
f(0)=y\in\mathbb{R}^{d},$

evaluated at $t=1$, see [13] and [1, Ch.4,6].

###  2.1 Special cases

We are specifically interested in the case $d=1$,
$\sigma(x)=\sigma\in\mathbb{R}^{m}$ constant and
$\nu=\frac{C}{|z|^{m+\alpha}}dz$, where $C$ is either a constant depending on
$m$ and $\alpha\in(0,2)$, or a function of $z$. In this case, the solution to
(2.3) is given by $g(t,x)=f(x+\sigma\cdot zt)$, and the last integral term in
(2.2) simplifies to

$\int_{0}^{t}\int_{\overline{B_{1}(0)}}u(s-,x+\sigma\cdot z)-u(s-,x)-\nabla
u(s-,x)\sigma z\,d\nu(z)ds.$

Moreover, in this case we have $\varphi_{0,t}(x)=x-\sigma\cdot Z_{t}$. Since
for any choice of $\sigma$ and $\nu$ both stochastic integrals from (2.2) are
martingales, taking expectation yields

$\mathbb{E}[u(t,x)]-\mathbb{E}[u_{0}(x)]=\mathbb{E}\bigg{[}\int_{0}^{t}\int_{\overline{B_{1}(0)}}u(s-,x+\sigma\cdot
z)-u(s-,x)-\nabla u(s-,x)\sigma z\,d\nu(z)ds\bigg{]},$ (2.4)

where we write $\mathbb{E}[X]=\int_{\Omega}X\,d\mathbb{P}$ for a random
variable $X:\Omega\to\mathbb{R}$, provided the integral is defined. Also note
that for $d=1$, every divergence-free vector field is constant.

##  3 Averaged enhanced dissipation

Let $d=1$, $\sigma$ be constant, and $\nu$ have a radially symmetric density
(for instance, the classical symmetric $\alpha$-stable density
$\frac{1}{|z|^{m+\alpha}}$), and set $\theta:=|\sigma|$. Then, due to the
radial symmetry, the RHS of (2.4) without expectation, i.e. for each fixed
$\omega\in\Omega$, equals

$\displaystyle\int_{0}^{t}\int_{\overline{B_{1}(0)}}u(s-,x+\theta
z_{1})-u(s-,x)-\nabla u(s-,x)\theta z_{1}\,d\nu(z)ds,$ (3.1)

where we denote by $z_{i}$ the $i$-th component of
$z=(z_{1},\dots,z_{m})\in\mathbb{R}^{m}$. In order to further calculate this
integral, we need the following lemma. Below, we denote by $\pi_{1}$ the
canonical projection $\pi_{1}:\mathbb{R}^{m}\to\mathbb{R}$,
$\pi_{1}(z)=z_{1}$.

######  Lemma 3.1.

Let $\alpha\in(0,2)$.

  1. (i)

Set $\nu_{\alpha}:=\frac{1}{|z|^{m+\alpha}}dz$. Then

$\nu_{\alpha,1}:=\nu\circ\pi_{1}^{-1}=\frac{C(m,\alpha)}{|y|^{1+\alpha}}dy,$

with

$C(m,\alpha)=|\mathbb{S}^{m-2}|\int_{0}^{\infty}(1+r^{2})^{-\frac{m+\alpha}{2}}r^{m-2}dr<\infty,$

where $|\mathbb{S}^{m-2}|$ denotes the surface area of $\mathbb{S}^{m-2}$, the
unit sphere in $\mathbb{R}^{m-1}$.

  2. (ii)

Let
$\nu_{\mathds{1},\alpha}:=\mathds{1}_{\overline{B_{1}(0)}}(z)\frac{1}{|z|^{m+\alpha}}dz$.
Then

$\nu_{\mathds{1},\alpha,1}:=\nu_{\mathds{1},\alpha}\circ\pi_{1}^{-1}=\frac{C(y,m,\alpha)}{|y|^{1+\alpha}}dy,$

with

$C(y,m,\alpha):=\mathds{1}_{[-1,1]}(y)|\mathbb{S}^{m-2}|\int_{0}^{\frac{\sqrt{1-y^{2}}}{|y|}}(1+r^{2})^{-\frac{m+\alpha}{2}}r^{m-2}dr,\quad
y\in\mathbb{R}.$

###### Proof.

  1. (i)

Let $A\in\mathcal{B}(\mathbb{R})$, and for
$z=(z_{1},\dots,z_{m})\in\mathbb{R}^{m}$, write
$z^{\prime}=(z_{2},\dots,z_{m})$. Then

$\displaystyle\nu_{\alpha,1}(A)=\nu_{\alpha}(A\times\mathbb{R}^{m-1})$
$\displaystyle=\int_{A}\int_{\mathbb{R}^{m-1}}(z_{1}^{2}+|z^{\prime}|^{2})^{-\frac{m+\alpha}{2}}dz_{1}dz^{\prime}$
$\displaystyle=\int_{A}|z_{1}|^{-m-\alpha}\int_{\mathbb{R}^{m-1}}\bigg{(}1+\frac{|z^{\prime}|^{2}}{z_{1}^{2}}\bigg{)}^{-\frac{m+\alpha}{2}}dz^{\prime}dz_{1}$
$\displaystyle=\int_{A}|z_{1}|^{-1-\alpha}\int_{\mathbb{R}^{m-1}}(1+|z^{\prime}|)^{-\frac{m+\alpha}{2}}dz^{\prime}dz_{1}$
$\displaystyle=\int_{A}C(m,\alpha)|z_{1}|^{-1-\alpha}dz_{1},$

where the third equality follows from the transformation rule and the final
one by calculating the inner integral by spherical coordinates.

  2. (ii)

Due to the definition of $\nu_{\mathds{1},\alpha}$, the proof is similar to
the first part.

∎

To convey our further procedure, first consider (3.1) with domain of
integration $\mathbb{R}^{m}$ instead of $\overline{B_{1}(0)}$, and choose
$\nu=C(m,\alpha)^{-1}\nu_{\alpha}$. Then, by Lemma 3.1 (i), (3.1) equals

$\int_{0}^{t}\int_{\mathbb{R}}\frac{u(s-,x+\theta y)-u(s-,x)-\nabla
u(s-,x)\theta y}{|y|^{1+\alpha}}\,dyds.$

Inserting in (2.4) and interchanging the expectation with the temporal and
spatial integral as well as with the gradient shows that
$U(t,x):=\mathbb{E}[u(t,x)]$ solves

$\partial_{t}U(t,x)=\mathcal{L}_{\alpha}U(t-,x),$ (3.2)

where the operator $\mathcal{L}_{\alpha}$ is defined by

$\mathcal{L}_{\alpha}f(x):=\int_{\mathbb{R}}\frac{f(x+\theta y)-f(x)-\nabla
f(x)\theta y}{|y|^{1+\alpha}}dy.$

Note that (3.2) is a deterministic, nonlocal second-order parabolic equation
for the expected value of $u(t,x)$, which itself solves (pathwise) (2.1). We
point out the similarity of $\mathcal{L}_{\alpha}$ with the fractional
Laplacian $(-\Delta)^{\alpha}$ on $\mathbb{R}$,

$(-\Delta)^{\alpha}f(x)=\int_{\mathbb{R}}\frac{f(x+y)-f(x)-\nabla
f(x)y\mathds{1}_{(-1,1)}(y)}{|y|^{1+\alpha}}dy.$

In fact, for $\theta=1$ and $\alpha>1$,
$\mathcal{L}_{\alpha}=(-\Delta)^{\alpha}$, since in this case
$\int_{B_{1}(0)^{c}}\frac{y}{|y|^{1+\alpha}}dy=0$.

Now, in order to take into account the proper domain of integration
$\overline{B_{1}(0)}$ in (3.1), we repeat the previous lines with the choice

$\nu=|\mathbb{S}^{m-2}|^{-1}\nu_{\mathds{1},\alpha}.$ (3.3)

Then, similarly to (3.2), we arrive at

$\partial_{t}U(t,x)=\mathcal{L}_{\mathds{1},\alpha}U(t-,x),\quad(t,x)\in\mathbb{R}_{+}\times\mathbb{R},$
(3.4)

satisfied by $U(t,x)=\mathbb{E}[u(t,x)]$, where we set

$\mathcal{L}_{\mathds{1},\alpha}f(x):=\int_{-1}^{1}f(x+\theta y)-f(x)-\nabla
f(x)\theta
y\,d\nu_{\mathds{1},\alpha,1}(y)=\int_{-1}^{1}c(m,\alpha,y)\frac{f(x+\theta
y)-f(x)-\nabla f(x)\theta y}{|y|^{1+\alpha}}\,dy$

with

$c(y,m,\alpha):=\int_{0}^{\frac{\sqrt{1-y^{2}}}{|y|}}(1+r^{2})^{-\frac{m+\alpha}{2}}r^{m-2}dr.$

For the finite positive weight $c(y,m,\alpha)$, we note
$c(y,m,\alpha)\xrightarrow{|y|\to 1}0$, and that $c(y,m,\alpha)$ is symmetric
around $y=0$. Since $c(y,m,\alpha)$ is bounded on $(-1,1)$, by Taylor formula
$\mathcal{L}_{\mathds{1},\alpha}f(x)$ is well-defined and finite for any
$\alpha\in(0,2)$, $x\in\mathbb{R}$, and $f\in C^{2}(\mathbb{R})$. The latter
is satisfied for $x\mapsto u(t-,x)$, for every $t\geq 0$ and $\mathbb{P}$-a.e.
$\omega\in\Omega$, as well as for $x\mapsto U(t-,x)$, for every $t\geq 0$.
Interchanging the expectation with the temporal and spatial integral in (3.1)
is justified, since $u$ is given as in Proposition (2.2) and since $\varphi$
is a stochastic flow of smooth diffeomorphisms, see [1, Thm.6.10.10].
Therefore, we have arrived at the following result.

######  Proposition 3.2.

Consider (2.1) for $d=1$, $u_{0}\in C^{2}_{b}(\mathbb{R})$,
$\sigma\in\mathbb{R}^{m}$ constant, $Z$ with Lévy measure
$\nu_{\mathds{1},\alpha}$ as in 3.3, and let $u$ be the unique solution in the
sense of Definition 2.1. Then, $U(t,x):=\mathbb{E}[u(t,x)]$ solves (3.4).

######  Remark 3.3.

Note that equation (3.4) is dissipative since it involves a fractional
Laplace-like operator. For such operators, there are results for the decay of
solutions, yielding a link with our numerical results (see Section 4) and the
operator $\mathcal{L}_{\mathds{1},\alpha}$. More precisely, consider

$\partial_{t}u(t,x)=(-\Delta)^{\alpha/2}u(t,x),\quad
u(0,x)=u_{0}(x),\quad(t,x)\in\mathbb{R}_{+}\times\mathbb{R}.$

There exists a $C_{0}-$semigroup $(\mathcal{S}_{\alpha}(t))_{t\geq 0}$ such
that $u(t,x)=\mathcal{S}_{\alpha}(t)*u_{0}(x)$. The following estimate in
dimension $d=1$ was obtained in [4].

$\|u(t)\|_{\infty}=\|\mathcal{S}_{\alpha}(t)*u_{0}\|_{\infty}\leq
t^{-\alpha/2}\|u_{0}\|_{1},$

where $\|\cdot\|_{\infty}$ and $||\cdot||_{1}$ denote the usual $L^{\infty}$\-
and $L^{1}$-norms, respectively, thus showing the dissipative behavior of the
solution. This behavior is retrieved also numerically (see below), suggesting
a similar behavior for $\mathcal{L}_{\mathds{1},\alpha}$ due to is similarity
with the fractional Laplacian.

So far, we were unable to prove similar results for individual paths of the
solution to (2.1). In fact, the pathwise profile shows no sign of
dissipativity (see Figure 1) in our simple special case, where $\sigma$ is
constant and $d=1$. So, a conclusion on the pathwise dissipative behavior and
an Ito-Stratonovich diffusion limit-type result is yet to be reached. A main
reason for this appears to be the absence of any mixing property of
divergence-free vector fields in dimension $d=1$, yielding the pathwise
profile a simple translation in time. We expect that in dimension $d>1$ and
for suitable vector fields as in [6, 5], a pathwise dissipative behavior of
the solution to (2.1) solution is possible.

##  4 Numerical results

Figure 1: Solution trajectory at several times

In order to support the results of the previous section, here we present
numerical simulations of the SPDE (2.1) in our special case, thereby
underlining the dissipative character of the expected value of its solution.
As in Section 3, we limit ourselves to dimension $d=1$, we consider a constant
vector field $\sigma\in\mathbb{R}^{m}$, and we choose the Lévy measure of the
driving pure jump Lévy process to be $\nu_{\mathds{1},\alpha}$, as defined in
Lemma 3.1 (ii). The jump dimension $m\in\mathbb{N}$ of the Lévy process is
arbitrary. Under these assumptions, an explicit solution to (2.1) can be
computed via the method of characteristics, which can be exploited for simple
numerical simulations, namely

$\displaystyle u(t,x;\omega):=u_{0}(x+\sigma\cdot Z_{t}(\omega)).$

For all simulations below, we fixed a smooth bump function
$u_{0}(x):=\exp(-\frac{0.01}{0.01-\min(0.01,x^{2})})$ as initial condition. To
simulate our $\alpha$-stable Lévy process trajectory, we need to take into
account the fact that our choice of Lévy measure $\nu_{\mathds{1},\alpha}$
neglects large jumps. To this end, using the independent increment and self-
similarity property, we compute the next step of the trajectory by cutting
away jumps larger than one and generating a new realization until the jump
size is sufficiently small. To implement the $\alpha-$stable distribution we
used the $R$-package $Stabledif$ [14]. The time step is $dt=10^{-4}$, in the
range $[0,2]$, while the space domain is selected to be the interval $[-1,1]$.
Note that the equation 2.1 is posed on $\mathbb{R}^{d}$, i.e. there is a
slight discrepancy between the equation and our numerical simulations that
needs to be taken into account when we interpreting the results. The support
of $\operatorname*{supp}u_{0}$ is sufficiently small so that almost no mass
escapes the system throughout the time evolution. The space discretization
size is selected at $dx=10^{-3}$.

(a) Decay in space of averaged solution at several times

(b) Decay in time of averaged solution at $x=0$

Figure 2: Averaged solution to equation 2.1

In our numerical analysis, we set our parameters in the range
$\theta=\|\sigma\|\in[0,1]$, $\alpha\in(0,2)\setminus\\{1\\}$ and
$m=1,2,5,10$. Our results are qualitatively consistent, and we are here
analyzing and presenting figures for $\theta=0.5,\ \alpha=1.5$ and $m=2$. In
future works, we expect to give a precise quantification of the dependence on
these parameters or the decay rate in time of the averaged solution. Our
results, as of now, focus on the qualitative behavior in time for the averaged
solution. This behavior depends only on the parameter $\alpha$ of the stable
distribution.

(a) Nonlinear regression $\sim t^{-1/alpha}$, $R.E.=0.003$

(b) Log-Log Plot and tail behavior. $z:=t\ (time),\ y:=\mathbb{E}[u(t,x=0)]$.
Left: full decay. Right: zoom on tail for larger times.

Figure 3: Nonlinear regression on 3.4

Figure 1 presents a solution trajectory to (2.1) at different times.
Precisely, the black curve represents the initial condition, while the blue
ones show the solution profile at times $t=1/2,\ t=1$. Similarly, the grey
curves represent the times $t=3/2,\ t=2$. As expected, the solution
trajectories are translations, with no mixing property arising from the
constant vector field $\sigma$, thus not showing any dissipation. The $L^{2}$
norm in space is preserved, and even for large times the solution preserves
energy with an error of $10^{-14}$, which is only due to the space domain
being finite in the simulations. Nonetheless it is clear that the energy is
preserved, as expected, due to the transport character of (2.1).

To show consistency with the theoretical results, we simulate the evolution of
the averaged solution profile, obtained with Monte Carlo methods over
averaging $5000$ samples, as pictured in Figures 2(a), 2(b), showing a decay
in both time and space. In the first figure, analogously to the pathwise
result, we present the averaged solution at different times. The black graph
represents the initial condition, while the green and purple ones show the
average solution at times $t=1/2,\ t=1$, respectively, and at times $t=3/2,\
t=2$ for the red and blue ones. In this case, a diffusive behavior is present
and a decay in space and time is observed. Concerning the space behavior, the
profiles are still not smooth. The reason is twofold, one being the averaging
procedure, the second one arising from the fractional operator obtained for
the equation modeling the averaged system (i.e. (3.4)) and the corresponding
small jumps. Concerning the decay, particularly care for $x=0$, in which the
initial condition has its maximum, and note the decay of the averaged solution
in time with a power law-like asymptotic behavior.

More so, starting with a compact supported initial condition, we show, in
Figures 3(a) and 3(b), the time decay in the origin $x=0$, with a nonlinear
regression to estimate the asymptotic and power law-like behavior, linking it
to the operator $\mathcal{L}_{\mathds{1},\alpha}$ proposed in the theoretical
section. In particular, in Figure 3(a), we performed a regression using the
insight of Remark 3.3 numerically to show that, with a residual error of less
than $0.003$, we have a decay of the averaged solution in time of the
following form:

$\mathbb{E}[u(t,x=0)]\sim\beta(\sigma,m,\alpha)t^{-1/\alpha},$

with $\beta$ depending on $\theta:=\|\sigma\|$, the dimension of the Lévy
process $m\in\mathbb{N}$, and the parameter $\alpha$ of the Lévy measure
$\nu$. Here, we have not delved into an analysis of the behavior of $\beta$,
which could be theoretically examined, as discussed in Remark 3.3. Therefore,
it is crucial for future research to explore how the strength of $\beta$, and
consequently the velocity field, interact with the decay of the profile.

In Figure 3(b), we plot the log-log version of the decay in time, showing the
inverse asymptotic behavior of the averaged solution profile in space and time
and its concordance with the results of the theoretical section and the
nonlinear regression. More so, on the right in Figure 3(b), the tail was
analyzed in the time frame $z[500:1001]:=[0.5,1]$, showing a rough behavior,
but with respect to the regression line, the error is in the range
$0.003-0.006$, which validates our results.

##  5 Open questions

We conclude this work with a few questions for future work in this direction.

  1. (i)

Can similar results be obtained for non-constant $\sigma$ (in dimension $d\geq
2$)?

  2. (ii)

Related to (i), do vector fields with suitable mixing properties lead to a
pathwise dissipation result for solutions to (2.1)?

  3. (iii)

Is it possible to obtain precisely the fractional Laplacian as the operator of
the deterministic equation (3.4) instead of $\mathcal{L}_{\mathds{1},\alpha}$,
and what is the corresponding Lévy measure?

  4. (iv)

Which operators does one obtain by choosing Lévy measures in (2.1) which are
not $\alpha$-stable?

#### Acknowledgements.

The research of F.F. is funded by the European Union (ERC, NoisyFluid, No.
101053472. A.P. is supported by the European Union (ERC, StochMan, No.
101088589). M.R. is funded by the German Research Foundation (DFG) - Project
number 517982119. The authors would also like to thank Marvin Weidner for a
very valuable comment regarding Lemma 3.1 and Gaia Tramonte for helpful
suggestions.

## References

  * [1] D. Applebaum.  Lévy Processes and Stochastic Calculus.  Cambridge Studies in Advanced Mathematics. Cambridge University Press, 2nd edition, 2009.
  * [2] A. Cavalieri, É. Mémin, and G. Tissot.  Input-output analysis of the stochastic Navier–Stokes equations: Application to turbulent channel flow.  Phys. Rev. Fluids, 8:033904, 2023.
  * [3] A. Chechkin and I. Pavlyukevich.  Marcus versus Stratonovich for systems with jump noise.  Journal of Physics A: Mathematical and Theoretical, 47(34):342001, 2014.
  * [4] A. Fino and M. Kirane.  The Cauchy problem for heat equation with fractional Laplacian and exponential nonlinearity.  Communications on Pure and Applied Analysis, 19(7):3625–3650, 2020\.
  * [5] F. Flandoli, L. Galeati, and D. Luo.  Quantitative convergence rates for scaling limit of SPDEs with transport noise.  arXiv preprint 2104.01740, 2021.
  * [6] F. Flandoli, L. Galeati, and D. Luo.  Eddy heat exchange at the boundary under white noise turbulence.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 380(2219):20210096, 2022.
  * [7] F. Flandoli, R. Huang, and A. Papini.  Turbulence enhancement of coagulation: The role of eddy diffusion in velocity.  Physica D: Nonlinear Phenomena, 448:133726, 2023.
  * [8] F. Flandoli and E. Luongo.  The dissipation properties of transport noise.  In Stochastic Transport in Upper Ocean Dynamics, pages 69–85, Cham, 2023. Springer International Publishing.
  * [9] F. Flandoli and E. Luongo.  Stochastic partial differential equations in fluid mechanics, volume 2330 of Lecture Notes in Mathematics.  Springer, Singapore, 2023.
  * [10] F. Flandoli, S. Morlacchi, and A. Papini.  Effect of transport noise on Kelvin–Helmholtz instability.  In Stochastic Transport in Upper Ocean Dynamics II, pages 29–52, Cham, 2024. Springer Nature Switzerland.
  * [11] F. Flandoli and F. Russo.  Reduced dissipation effect in stochastic transport by Gaussian noise with regularity greater than 1/2.  arXiv preprint 2305.19293, 2023.
  * [12] L.-S. Hartmann and I. Pavlyukevich.  First order linear Marcus SPDEs.  arXiv preprint 2303.00674, 2023.
  * [13] S. Marcus.  Modeling and analysis of stochastic differential equations driven by point processes.  IEEE Transactions on Information Theory, 24(2):164–172, 1978.
  * [14] M. Mächler, D. Würtz, and Rmetrics core team members.  stabledist: Stable distribution functions.  2016-09-12.


# The deep learning and statistical physics applications to the problems of
combinatorial optimization.111National Supercomputing Forum (NSCF-2019),
Pereslavl-Zalessky, Program Systems Institute of the RAS, Russia.

Semyon Sinchenko  Moscow Aviation Institute, Volokolamskoe shosse, 4, Moscow,
125993, Russian Federation  Dmitry Bazhanov  Moscow Aviation Institute,
Volokolamskoe shosse, 4, Moscow, 125993, Russian Federation  Faculty of
Physics, Moscow State University, GSP-1, Leninskie Gory, Moscow, 119991,
Russian Federation  Institution of Russian Academy of Sciences, Dorodnicyn
Computing Centre of RAS, Vavilov St. 40, Moscow, 119333, Russian Federation

###### Abstract

We present herein a new approach based on the simultaneous application of the
deep learning and statistical physics methods to solve the combinatorial
optimization problems. The recent modern advanced techniques, such as an
artificial neural network, demonstrate their efficiency for solving various
physical tasks for the quantum many-body systems, which may be related
directly to the problems of combinatorial optimization. One of them is a
classical Maximum Cut (MaxCut) problem, which we ascribe here to the search of
the ground state of a quantum many-body system using an artificial neural
network and deep learning. We found that the exact solution received for a
quantum system corresponds to its counterpart in a classical MaxCut problem.
As a proof, we have realized our approach for two random graphs of different
size, containing 60 vertices (885 edges) and 100 vertices (2475 edges), and
achieved for them the total performances $p_{60}\sim 0.99$ and $p_{100}\sim
0.97$ of known maximal cuts, respectively.

##  1 Introduction

###  1.1 The MaxCut problem

The MaxCut problem is nowadays a well-studied problem of combinatorial
optimization and integer programming. The problem can be defined as follows:
for a given Graph $G=(V,E)$222Herein, we consider for simplicity the
unweighted and undirected graphs only. However, we will discuss finally about
the generalization of our ansatz., where $V$ is the set of vertices and $E$ is
the set of edges, the MaxCut problem seeks to find the subsets of vertices
$V_{1}$ and $V_{2}$ in such a way that the number of edges, which incident to
vertices from different sets, are maximal:

$\sum_{(i,j)\in E}(1-x_{i}x_{j})\to Max,$

where $x_{i}=1$ if the vertex from $V_{1}$ set and $x_{i}=-1$ if the vertex
from $V_{2}$ one.

In order to describe MaxCut problem one can define also the performance
parameter $p_{A}$ for a given algorithm $A$ that exhibits the cut of size $C$
close to the optimal value of cut $MCS$:

$p_{A}=\frac{C}{MCS}$

###  1.2 The well known approaches

The MaxCut problem is known to be NP-hard problem and its exact solution, that
can be found by the branch-and-bound algorithms, is unfeasible in the case of
medium or large graph. Some of these algorithms are successful in solving the
MaxCut problem, but they achieve the worth solution for performance parameter
about $p_{A}\sim 0.90$ [1]. Besides that, there was a theoretically predicted
upper bound that guarantees the performance of algorithms to be not more than
$p_{A}^{b}\sim 0.941$[2]. Nevertheless, during the last years the novel
approaches based on the deep Neural Network (NN) and Machine Learning (ML)
have been developed and intensively applied for solving the combinatorial
optimization problems[3]. Some of them [4][5] were employing, so-called, Graph
Convolutional Networks (GCN), which represent the new class of NN introduced
specifically to solve the graph semi-supervised learning problems[6]. Another
ones, based on Reinforcement-Learning (RL)[7], were involving the action space
and reward for the applied algorithm. The latter helps the algorithm to learn
and seek the best actions that maximize the reward. Besides, there is a wide
class of algorithms, which are built up on Monte Carlo Markov Chains (MCMC)
[8]. This class of algorithms does not guarantee anyone the nice solution in a
reasonable amount of time, but demonstrates a very good performance in some
combinatorial optimization problems, when the built-in hyperparameters are
properly tuned. One of the best algorithm of this class is a well known and
named as Simulated Annealing (SA)[9]. In the framework of SA one can employ
the random permutation of the current state at each iteration step. In our
case this state can be ascribed to $V_{1}$ or $V_{2}$ set. Then a new state is
drawn from the probability distribution $p$:

$p=P_{Gibbs}(E_{cur},E_{new},T),$

where $E_{cur}(E_{new})$ presents the current (new) state endowed by some
energy variable and $P_{Gibbs}$ is a conventional Gibbs distribution, written
as:

$P_{Gibbs}(E_{cur},E_{new},T)=\begin{cases}1,&\mbox{if }E_{new}\geq E_{cur}\\\
\exp{\frac{E_{new}-E_{cur}}{T}},&\mbox{otherwise}\end{cases}$

where $T$ is a temperature parameter, which should be tending slowly down, for
example as a function ($T\sim\frac{1}{N}$) of learning steps $N$. For the
MaxCut problem the energy variable in SA algorithm can be defined as a cut
size of the considered graph. A flowchart in the Fig.1 illustrates a workflow
of SA process.

New permuted sample Calculate the energy Accept with probability $p\sim
P_{Gibbs}(E_{cur},E_{new},T)$ Decrease the temperature Figure 1: The flowchart
of Simulated Annealing process.

However, there were some cases reported before where the treatment based on SA
algorithm failed to find a good solution in a reasonable amount of time [10].
Specifically, some difficulties can arise within MCMC-based algorithms due to
their poor convergence at meaningful temperatures for real conditions.
Therefore the main goal of this work is an attempt to eliminate such drawbacks
of MCMC-based algorithms and to improve them by implementation of the latest
achievements in statistical physics and quantum theory regarded to the
problems of combinatorial optimization.

##  2 Brief description of our ansatz

In this work we present a novel physical-inspired approach to solve the MaxCut
problem for graphs by means of deep NNs. In contrast to GCN-based methods,
this approach can be easily generalized to a wide class of optimization
problems and not only for graphs, but for any other instances. In comparison
to RL-based methods we suggest to use more simple NN with the classical
gradient descent and minimization procedures that leads to more stable
learning process than RL-learning process. In our approach we treat the
problem of combinatorial optimization as a physical task that should be solved
for quantum many-body system with a special, so-called, cost Hamiltonian by
means of Deep Learning based techniques.

###  2.1 Quantum many-body problem

The searching of the wave function $\Psi$ of the many-body system is one of
the most important, but one of the hardest problems in quantum physics. The
wave function contains the whole information about the state of any quantum
many-body system (particles, spins, etc.) and the dynamic of the wave function
allows to predict the state of such system at the each moment. Formally the
wave function $\Psi$ is the complex-valued function and its square modulus
relates to the probability to find the system at state $\ket{\psi}$ (here
$\ket{x}$ is a ket notation for the vector in a complex Hilbert space) as the
result of the measurements. The wave function of the system can be found
analytically as a solution of the Schrödinger equation:

$i\hbar\frac{d}{dt}\ket{\Psi(t)}=\hat{H}\ket{\Psi(t)},$

where $\hat{H}$ is the Hamiltonian operator that involves all interactions
between particles of the whole system. Therefore the amount of information
that is necessary to be stored for full encoding of the system increases
exponentially with the number of particles. Thereby, the necessity in
computation the integrals or sums (in discrete cases) over all Hilbert space
fails to use the direct approaches for solving the problems of quantum many-
body systems. One of the most popular solution is a Quantum Monte Carlo method
that uses the variations of the Monte-Carlo method to estimate integrals over
configuration space in polynomial time[11]. Another popular approach named
Matrix Product State is the generalization of numeric renormalization-group
procedure to quantum lattice, and also the more general Tensor Network
approach. However, there are some instances of quantum systems where the
mentioned approaches fail[13].

###  2.2 The link between MaxCut and quantum Hamiltonian

Several attempts to solve the quantum many-body problem by means of
combinatorial optimization methods are known from the state-of-the-art
study[14].

Let us assume that any arbitrary quantum system, that consists of $|V|$
particles with quantum spins $s=\\{+1,-1\\}$, is ascribed for a given graph
$G$ and defines the cost Hamiltonian for this graph in the following manner:

$\hat{H}=\sum_{(i,j)\in E}A_{ij}\sigma^{z}_{i}\sigma^{z}_{j},$

where $A$ is the adjacency matrix of the graph and $\sigma^{z}$ is the Pauli
matrix:

$\sigma^{z}=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}$

The Hamiltonian is diagonalizable in $z$-basis that allows to solve the MaxCut
problem more efficiently.

It is well known, that the ground state of a quantum system is the state with
minimum energy (or minimized Hamiltonian). Thereby, if from the given graph
the subsets of $V_{1}$ and $V_{2}$ vertices are determined as
$V_{1}=\\{x_{i}:s_{i}=1\\}$ and $V_{2}=\\{x_{i}:s_{i}=-1\\}$, where $s$ is the
particle spin, then the search for the ground state energy
$\Large\varepsilon_{0}$ of the quantum system with Hamiltonian $\hat{H}$ leads
to the solution of the MaxCut problem:

$C=-\frac{\Large\varepsilon-\sum_{(i,j)\in E}1}{2},$

where $\Large\varepsilon$ is an eigenvalue of $\hat{H}$. The reverse statement
is also true: if the MaxCut problem is solved for the given graph, the ground
state of the corresponding quantum system is also found. Such approach is
being actively developed nowadays in the area of Quantum Information Theory
and Quantum Algorithms and shows an excellent results in experimental quantum
computations for small graphs[15][16].However, no really scalable industrial
quantum computer still exists.

###  2.3 Neural Quantum States

The new approach to solve the quantum many-body problem, which was reported in
[17], is based on the Neural Quantum States (NQS). According to this work, one
can define a NN with one hidden layer as function $F(\bf{S})$ where $S$ is the
spin configuration vector: $S=\\{s_{i}\\}$. That function is fully described
by the matrix of weights $\bf{W}$, vector of biases $\bf{b}$ and the
activation functions $\sigma$333For example sigmoid-function:
$\sigma(x)=\frac{1}{1+e^{-x}}$:

$F(\bf{S})=\sigma(\bf{b}+\bf{W}\dotproduct{\bf{S}})$

Owing to Cybenko universal approximation theorem any function
$[0,1]^{n}\to{[0,1]}$ can be approximated by NN with one hidden layer with any
precision[18]. There is the generalization of that theorem to complex
plane[19]. Thus, hypothetically, for every real quantum system, there is a NN
with weights $\bf{W}$ and biases $\bf{b}$, so that the output of this network
approaches the $\Psi$ function of this system.

We can define the deep NN (in our case it was a simple Feed-Forward Neural
Network (FFNN)), which brings us a predictor444A black-box model that receives
a state as input and returns a $\Psi$ function as output. This predictor can
give the probability of any state of our system. Thereby, we can get the
samples of system configurations $S=\\{s_{i}\\}$ drawn from the probability
distribution defined by our variational $\Psi(S)$ function by means of this
predictor and employing MCMC-methods[8]. From physical principles, we know
that every physical system strives to have a gain in energy and tends to its
ground state with minimum energy. It means that we will get the states close
to energy minimum if we draw the samples of system configurations from the
true $\Psi$ function. Consequently, this allows us to determine the energy-
minimization procedure:

  * •

To draw the samples of system configurations $S=\\{s_{i}\\}$ from the
probability distribution by means of predictor and employing FFNN;

  * •

To estimate an energy of the samples $E_{loc}(S)$ and the stochastic
expectation value of this energy
$\langle\hat{H}\rangle\simeq\langle\langle{E_{loc}}\rangle\rangle$;

  * •

To estimate the gradient of energy $\partial_{p}\langle{\hat{H}}\rangle$ (see
in Appendix $A$ for more details);

  * •

To update the weights $\bf{W}$ and biases $\bf{b}$ of FFNN by using the
Gradient-Descent methods;

###  2.4 In comparison with Simulated Annealing

In principle, the closest approach to our NQS-based ansatz is SA algorithm,
where the heuristic Gibbs distribution is essentially replaced by the FFNN,
which receives a state as input and returns a complex-valued $\Psi$ function
as output. We use the squared modulus of $\Psi$ function to calculate the
acceptance probability $p$:

$p=\Bigl{|}\frac{\Psi(S^{\prime})}{\Psi(S)}\Bigr{|}^{2}$

Using this probability, we can redefine the sampling procedure. To do so far,
we apply an analogue of the annealing process (here we used Metropolis-
Hastings algorithm) at each learning epoch under implementation of deep NN and
calculate the acceptance probability instead of using the Gibbs one. Besides,
at the end of learning epoch we update the weights of deep NN based on
stochastic reconfiguration procedure [20]:

$\displaystyle G=\frac{\partial_{W}{\Psi(S)}}{\Psi(S)}$ $\displaystyle
E_{loc}=\sum_{S`}H_{SS`}\frac{\Psi(S`)}{\Psi(S)}\text{}$ $\displaystyle
SR=\langle{G^{T}G}\rangle-\langle{G}\rangle\langle{G}\rangle^{T}$
$\displaystyle
F=\langle{GE_{loc}}\rangle-\langle{G}\rangle\langle{E_{loc}}\rangle$
$\displaystyle W_{n+1}=W_{n}-\alpha SR^{-1}F\text{}$

We also discard some first samples of system configurations at each learning
epoch to make the learning process more stable. A detailed flowchart of the
approach presented is shown in the Fig.2.

New permuted state: $S^{\prime}$ Calculate $\Psi(S^{\prime})$ function that is
the output of deep NN Accept with probability
$p\sim\Bigl{|}\frac{\Psi(S^{\prime})}{\Psi(S)}\Bigr{|}^{2}$ Add a new sample
to collection Estimate the gradient of energy Update deep NN weights and
biases Figure 2: The flowchart of the developed approach.

##  3 Workflow details

We have employed two graphs from the Biq-Mac collection: the first graph named
g05_60.0 (60 vertices, 884 edges) and the second graph named g05_100.0 (100
vertices, 2474 edges)[21]. Both graphs are illustrated in the Fig.3 using
NetworkX library[22]. For each graph the corresponding quantum Hamiltonian has
been determined and the NQS approach was applied to find its ground state
energy $E_{0}$. The detailed description of the hyperparameters for NQS is
performed in Appendix $A$. At the end of the learning process, when deep NN
yields a nearly true $\Psi$ function, we have drawn 2000 samples from the
probability distribution $p$. Then the cut size ($C$) was computed for each
sample and compared with benchmark solution ($MCS$) of the MaxCut for a given
graph.

(a) 60-vertex graph

(b) 100-vertex graph

Figure 3: The sketches of employed graphs

We have used also the NetKet toolkit in order to estimate the ground state
energy $E_{0}$ and $\Psi$ function of a given quantum Hamiltonian[23]. This
toolkit presents a high-performance library written in $C++$ and containing
the Python-bindings. The fully developed source code and Dockerfile with
NetKet are available in the GitHub repository[24].

##  4 Results

We present in the Fig.4 and Fig.5 the progression curves corresponding to the
learning processes that 60-vertex and 100-vertex graphs possessed,
respectively. The left plots of both figures demonstrate the progression of
cut sizes ($C$) of both graph instances, that corresponds to the expectation
values of energy $\langle\hat{H}\rangle$ of the cost Hamiltonians. Moreover,
the dashed lines on these plots represent the estimated optimal solutions of
MaxCut ($MCS$) for considering graphs (or exact ground state energies of the
cost Hamiltonians). It is well seen in both plots that, when the learning
process progresses with subsequent NN updates, the probability of sampling the
high quality solutions increases rapidly for both graph instances (learning
epochs $\gtrsim$ 100). The detailed energy variations with learning iterations
are demonstrated in the middle plots of the Fig.4 and Fig.5. It is worth to
note here, that for the physical systems described by ”good” Ising or
Heisenberg models of Hamiltonians the energy variations should be vanished,
since the samples produced during the NN learning are close to the optimal
solution for the ground state with a high probability[17]. However, for a
given graph the sample that differ from the ground state by only one vertex
can exhibit the significant energy variations, since the degree of this vertex
of a graph is not the same. Nevertheless, as we can see, there is no
difficulty with a large energy variations, because the expectation value of
the energy reaches an optimal solution anyway. The right plots of the Fig.4
and Fig.5 depict the acceptance ratio of Metropolis-Hastings algorithm. It can
be seen that the diversity of states decreases rapidly with the learning
process of our NN. This result indicates also that we are going towards to a
nearly true $\Psi$ function.

Figure 4: The learning curve (left), energy variation (middle) and acceptance
rate (right) for 60-vertex graph. Figure 5: The learning curve (left), energy
variation (middle) and acceptance rate (right) for 100-vertex graph.

Note that in contrast to a classical deep learning, where the learning process
can be controlled by an error parameter only (for instance, the energy
difference between the well known value and the value estimated for the
samples), in present approach there are three significantly different key
parameters: error, energy variation, and acceptance rate. If the energy
expectation value is within a plateau of the learning curve and the acceptance
rate is close to zero, one can conclude that we have achieved a very good fit
to the true $\Psi$ function and can start the sampling from it to generate the
optimal or nearly optimal solution to a given MaxCut problem. This is a very
important advantage of the approach realized in current work, since the
control of the fitting process in unsupervised learning (the optimum value of
cut ($MCS$) is unknown in reality) is a hard problem. Therefore, having looked
at the learning curves in the Fig.4 and Fig.5, we can conclude that the
learning progression up to 800 epochs is our overhead cost and can be
interrupted significantly earlier after $\sim 100-200$ learning epochs.

As previously claimed, we should return for solving a classical MaxCut problem
in order to receive the optimal solution for a given graph based on the
learning of NQS $\Psi(S)$ function. To do so, for each graph up to 2000
samples were drawn from the probability distribution corresponding to
$\Psi(S)$ function using the Metropolis-Hasting algorithm. Then, these samples
were used to find an optimal ordering of the high-quality solutions through
the calculation of the cut size (see Fig.6). As shown in the Fig.6, the
sampling up to 2000 samples reveals that most of the obtained solutions are
the high-quality solutions, especially for 60-vertex graph. These findings
result on the remarkable values of performance parameter obtained for both
given graphs across the considered sets of their samplings:

$p_{60}\sim 0.99\quad\mbox{for 60-vertex graph}$ $p_{100}\sim
0.97\quad\mbox{for 100-vertex graph}$

We want to emphasize that the similar results one can achieve by using SA
algorithm (or an other MCMC-based algorithms), but such approach can fail for
more sophisticated problems than those we have studied here[10]. We expect
that the employing of the NNs instead of Gibbs-sampling or any other prior
distribution can improve the convergence significantly.

(a) 60graph

(b) 100graph

Figure 6: The computed cut size ($C$) for 2000 samples drawn from the
probability distribution of 60-vertex graph(a) and 100-vertex graph(b),
respectively. The dotted line corresponds to the optimal value of cut ($MCS$)
for a given graph.

##  5 Outlook

To conclude, our findings unveil a new insight into the application of a novel
physical-inspired NQS-based approach to solve the various combinatorial
optimization problems by implementation of an artificial NN and deep learning
means. The presented approach reveals a nice convergence towards a high-
quality solutions of the MaxCut problem with a remarkable performance $p_{A}$
on two typical in size graph instances (60-vertex and 100-vertex graphs). It
demonstrates a noticeable supremacy over a classical heuristic algorithms of
optimization in the performance, acceptance rate, solution time, and
computational efficiency of simulation procedure. Moreover, the method of NQS
has а reliable theoretical background based on quantum and perturbation
theory[25] and has applied successfully to study the problems of the quantum
many-body systems [26, 27, 28, 29]. In spite of that, the rigorous theoretical
proof is still needed to ensure that our approach guarantees the performance
of the obtained results for more complex instances or optimization problems.

Finally, it is easy to show how to generalize our approach to the directed
weighted graphs:

$\hat{H}=\sum_{(i,j)\in E}A_{ij}\sigma^{z}_{i}\sigma^{z}_{j}$

Herein, we can simply define $A$ as an adjacency matrix with weights taken as
values and $E$ as a set of directed edges. Thereby, in a more general sense,
we can ascribe any unconstrained binary optimization problem to the problem of
minimizing the cost Hamiltonian of the corresponding Ising-like spin glass
model. This kind of approach was demonstrated in some works on the application
of quantum annealing and quantum optimization methods to a wide range of
combinatorial optimization problems[30].

It is worth to note, that there are also many advanced tools in the fast-
growing field of deep learning that can be successfully applied in NQS-based
approaches of combinatorial optimization. For example, one can use
convolutional layers[31] instead of fully connected layers, which can
significantly reduce the number of learning parameters and, therefore, reduce
the Vapnik-Chervonenkis dimension of NN[32] and significantly simplify the
computation complexity of learning. We assume that one can move from complex
weights of NN to real weights, which allow us to use modern deep learning
libraries, such as Tensorflow[33] and PyTorch[34]. These libraries make it
easy to determine the structure of NNs and to use fast GPU computing[35].

##  6 Acknowledgments

The authors are grateful for the financial support from Russian Foundation for
Basic Research (Grant No.19-29-03051-mk).

## References

  * [1] Bian, Y., Gronskiy, A., & Buhmann, J. M. (2016). Information-theoretic analysis of MaxCut algorithms. In 2016 Information Theory and Applications Workshop (ITA) (pp. 1–5).
  * [2] Håstad, J. (2001). Some optimal inapproximability results. Journal of the ACM, 48(4), 798–859.
  * [3] Bengio, Y., Lodi, A., & Prouvost, A. (2018). Machine Learning for Combinatorial Optimization: a Methodological Tour d’Horizon. ArXiv Preprint ArXiv:1811.06128.
  * [4] Gasse, M., Chételat, D., Ferroni, N., Charlin, L., & Lodi, A. (2019). Exact Combinatorial Optimization with Graph Convolutional Neural Networks. ArXiv Preprint ArXiv:1906.01629.
  * [5] Amizadeh, S., Matusevych, S., & Weimer, M. (2019). PDP: A General Neural Framework for Learning Constraint Satisfaction Solvers. ArXiv Preprint ArXiv:1903.01969.
  * [6] Kipf, T. N., & Welling, M. (2017). Semi-Supervised Classification with Graph Convolutional Networks. International Conference on Learning Representations.
  * [7] Amizadeh, S., Matusevych, S., & Weimer, M. (2019). PDP: A General Neural Framework for Learning Constraint Satisfaction Solvers. ArXiv Preprint ArXiv:1903.01969.
  * [8] Gilks, W. R., Richardson, S., & Spiegelhalter, D. (1997). Markov Chain Monte Carlo in Practice. Technometrics, 39(3), 338–338.
  * [9] Kirkpatrick, S., Gelatt, C. D., & Vecchi, M. P. (1983). Optimization by simulated annealing. Science, 220(4598), 671–680.
  * [10] Greig, D. M., Porteous, B. T., & Seheult, A. H. (1989). Exact Maximum A Posteriori Estimation for Binary Images. Journal of the Royal Statistical Society Series b-Methodological, 51(2), 271–279.
  * [11] Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21(6), 1087–1092.
  * [12] White, S. R. (1992). Density matrix formulation for quantum renormalization groups. Physical Review Letters, 69(19), 2863–2866.
  * [13] Troyer, M., & Wiese, U.-J. (2005). Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations. Physical Review Letters, 94(17), 170201.
  * [14] Barahona, F., Grötschel, M., Jünger, M., & Reinelt, G. (1988). An application of combinatorial optimization to statistical physics and circuit layout design. Operations Research, 36(3), 493–513.
  * [15] Crooks, G. E. (2018). Performance of the Quantum Approximate Optimization Algorithm on the Maximum Cut Problem. ArXiv Preprint ArXiv:1811.08419.
  * [16] Glen Bigan Mbeng, Rosario Fazio, Giuseppe Santoro (2019). Quantum Annealing: a journey through Digitalization, Control, and hybrid Quantum Variational schemes. ArXiv Preprint ArXiv:1906.08948
  * [17] Carleo, Giuseppe, and Matthias Troyer. ”Solving the quantum many-body problem with artificial neural networks.” Science 355.6325 (2017): 602-606.
  * [18] Cybenko, George. ”Approximations by superpositions of a sigmoidal function.” Mathematics of Control, Signals and Systems 2 (1989): 183-192.
  * [19] Kim, Taehwan, and Tülay Adalı. ”Approximation by fully complex multilayer perceptrons.” Neural computation 15.7 (2003): 1641-1666.
  * [20] Sorella, Sandro. ”Green function Monte Carlo with stochastic reconfiguration.” Physical review letters 80.20 (1998): 4558.
  * [21] Wiegele, A. (2007). Biq Mac Library—A collection of Max-Cut and quadratic 0-1 programming instances of medium size. Preprint.
  * [22] Hagberg, A. A., Schult, D. A., & Swart, P. J. (2008). Exploring Network Structure, Dynamics, and Function using NetworkX.
  * [23] Carleo, Giuseppe, et al. (2019). NetKet: A Machine Learning Toolkit for Many-Body Quantum Systems. ArXiv Preprint ArXiv:1904.00031
  * [24] https://github.com/SemyonSinchenko/QOpt
  * [25] Borin, A., & Abanin, D. A. (2019). Approximating power of machine-learning ansatz for quantum many-body states. ArXiv Preprint ArXiv:1901.08615.
  * [26] Choo, K., Neupert, T., & Carleo, G. (2019). Study of the Two-Dimensional Frustrated J1-J2 Model with Neural Network Quantum States. ArXiv Preprint ArXiv:1903.06713.
  * [27] McBrian, K., Carleo, G., & Khatami, E. (2019). Ground state phase diagram of the one-dimensional Bose-Hubbard model from restricted Boltzmann machines. ArXiv Preprint ArXiv:1903.03076.
  * [28] Yoshioka, N., & Hamazaki, R. (2019). Constructing neural stationary states for open quantum many-body systems. ArXiv Preprint ArXiv:1902.07006.
  * [29] Torlai, G., Mazzola, G., Carrasquilla, J., Troyer, M., Melko, R., & Carleo, G. (2018). Neural-network quantum state tomography. Nature Physics, 14(5), 1.
  * [30] S. Hadfield (2018). Quantum Algorithms for Scientific Computing and Approximate Optimization. ArXiv Preprint ArXiv:1805.03265.
  * [31] Cun, Y. L., & Bengio, Y. (1994). Word-level training of a handwritten word recognizer based on convolutional neural networks. In Proceedings of the 12th IAPR International Conference on Pattern Recognition, Vol. 3 - Conference C: Signal Processing (Cat. No.94CH3440-5) (Vol. 2, pp. 88–92).
  * [32] Vapnik, V., & Chervonenkis, A. J. (1974). Theory of Pattern Recognition.
  * [33] Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., … Isard, M. (2016). TensorFlow: a system for large-scale machine learning. Operating Systems Design and Implementation, 265–283.
  * [34] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., … Lerer, A. (2017). Automatic differentiation in PyTorch.
  * [35] Right now we are doing tests regarding the implementation of the NQS-solver for MaxCut problems in Tensorflow 2.0, in which the abilities of the cost Hamiltonian to be diagonalized and GPU-computed are used.

## Appendix A. Hyperparameters of NQS

In the original paper of G.Carleo and M. Troyer the Neural Quantum States
(NQS) approach was proposed to represent the wave function $\Psi$ of the many-
body quantum system using the Restricted Boltzmann Machine (RBM)[17]. Using
RBM there are some difficulties when the number of particles in a quantum
system increases essentially. For example, if the density of the hidden units
is only four, for a graph of 100 vertices there will be 400 hidden units. Such
amount of units produces 40000 weights in the hidden layer, which become an
operating problem for stochastic reconfiguration method: for 10000 samples at
each learning epoch the dense matrices of 40000$\crossproduct$10000 dimension
must be operated. To get around this problem we used simple Multi-Layer
perceptrons instead of RBM.

### Network Architecture for 60-vertex graph

For small graph we used the following structure:

  * •

Fully-Connected layer: input shape 60, output shape 30

  * •

Fully-Connected layer: input shape 30, output shape 20

  * •

Fully-Connected layer: input shape 20, output shape 10

  * •

LnCosh layer

  * •

Sum layer

Learning Hyperparameters:

Optimizer | Nesterov momentum ($lr=0.008,\beta=0.9$)
---|---
Samples at each step | 6000
Discarded samples at each step (first N) | 3000
Learning epochs | 400

### Network Architecture for 100-vertex graph

For small graph we used the following structure:

  * •

Fully-Connected layer: input shape 60, output shape 40

  * •

Fully-Connected layer: input shape 40, output shape 30

  * •

Fully-Connected layer: input shape 30, output shape 10

  * •

LnCosh layer

  * •

Sum layer

Learning Hyperparameters:

Optimizer | Nesterov momentum ($lr=0.008,\beta=0.9$)
---|---
Samples at each step | 10000
Discarded samples at each step (first N) | 5000
Learning epochs | 800