Complexity of Bose-Hubbard Model :
Quantum Phase Transition

Wung-Hong Huang*

Department of Physics, National Cheng Kung University,

No.1, University Road, Tainan 701, Taiwan

Abstract

The operator approach is applied to investigate the complexity of Bose-Hubbard
model. We present a systematic method to expand the quantum complexity in
series of coupling constant. We first study 2-sites system. For the ground
state we can find the exact value of complexity which is the summation of all
order. The complexity is divergent at critical value of coupling constant that
indicates the quantum phase transition at this point. We then generalize the
method to the N-sites closed chain and any dimensional system. The found
properties are similar to those in 2-sites system. We also study the excited
state and present the general formulas of Bose-Hubbard model complexity, which
shows a similar form as that in $\lambda\phi^{4}$ theory studied in our
previous paper.

* Retired Professor of NCKU, Taiwan.

* E-mail<EMAIL_ADDRESS>
###### Contents

  1. 1 Introduction
  2. 2 Bose Hubbard Model
  3. 3 Complexity of Two-Sites Bose Hubbard Model
    1. 3.1 Complexity of Free Theory
      1. 3.1.1 Ground State
      2. 3.1.2 Excited States
    2. 3.2 Complexity of Interacting Theory
      1. 3.2.1 First Order
      2. 3.2.2 n’th-Order
    3. 3.3 Complexity and Quantum Phase Transition
  4. 4 Complexity of N-Sites Bose-Hubbard Model
    1. 4.1 Complexity of Free Theory
    2. 4.2 Complexity of Interacting Theory
      1. 4.2.1 First Order
      2. 4.2.2 n’th-Order
    3. 4.3 Complexity and Quantum Phase Transition
  5. 5 Concluding Remarks

##  1 Introduction

Complexity is a quantity expected to probe the non-local property of a quantum
system besides the entangle entropy [1, 2]. In the context of the holography
the complexity=volume (CV) conjecture [3] and complexity=action (CA)
conjecture [4, 5] were proposed to calculate the quantity. The complexity in
there is the volume of an extremal codimension-one bulk surface and
gravitational action evaluated on the Wheeler-DeWitt (WDW) patch anchored on a
time slice of the boundary respectively.

Complexity have been studied in the quantum field theory [6, 7, 8, 9, 10, 11].
The complexity therein is defined as the number of operations $\\{\cal
O^{I}\\}$ needed to transform a reference state $|\psi_{R}\rangle$ to a target
state $|\psi_{T}\rangle$. These operators are also called as quantum gates :
the more gates one needs, the more complexity the target state has. To
calculate the complexity one defines an affine parameter “s” associated to an
unitary operator $U(s)$ and use a set of function, $Y^{I}(s)$, to character
the quantum circuit. The unitary operation connecting the reference state and
target state is

$\displaystyle U(s)=\vec{\cal P}\,e^{\int_{0}^{s}\,Y_{I}(s)\,{\cal
O}_{I}},~{}~{}~{}~{}|\psi_{R}\rangle=U(0)|\psi_{R}\rangle,~{}~{}~{}~{}~{}|\psi_{T}\rangle=U(1)|\psi_{R}\rangle~{}~{}~{}\
,$ (1.1)

where $\vec{\cal P}$ is a time ordering along $s$. The complexity ${\cal C}$
and circuit depth $D[U]$ (cost function) are [6]

$\displaystyle{\cal C}$ $\displaystyle=$
$\displaystyle\underset{\rm\\{Y^{I}\\}}{\rm
Min}\,D[U],~{}~{}~{}~{}~{}~{}D[U]=\int_{0}^{1}ds\,{\sum_{I}|Y^{I}(s)|^{2}}\
.~{}~{}~{}$ (1.2)

Above definitions were shown to be consistent with a gravitational computation
[6].

The initial studies in field theory considered the Gaussian ground state
wavefunctions in reference state and target state [6, 8, 9]. The theories
studied so far are the free field theory or exponential type wavefunction in
interacting model [10]. In this (wavefunction) approach, since that the
excited-state wavefunction of harmonic oscillation is not pure exponential
form the wavefunction approach is hard to work.

In our first paper [11] we adopt the operator approach, in which the
transformation matrix between the second quantization operators of reference
state and target state is regarded as the quantum gate, to evaluate the
complexity in free scalar field theory111The operator approach had also been
used in [7, 9] to study the complexity of fermion theory.. We examined the
system in which the reference state is two oscillators with same frequency
$\omega_{f}$ while the target state is two oscillators with frequency
$\omega_{1}$ and $\omega_{2}$. The complexity in excited states is calculated
and find that that the square of geodesic length in the general state $|{\rm
N_{1},N_{2}}\rangle$ is

$\displaystyle
D_{\rm(N_{1},N_{2})}^{2}={\rm(N_{1}+1)}\left(\ln{\sqrt{\omega_{1}\over\omega_{f}}}\,\right)^{2}+{\rm(N_{2}+1)}\left(\ln{\sqrt{\omega_{2}\over\omega_{f}}}\,\right)^{2}$
(1.3)

The results was furthermore extended to the N couple harmonic oscillators
which correspond to the lattice version of free scalar field.

In our next paper [12] we included interaction to further study the complexity
using the operator approach. We present a systematic method to evaluate the
complexity of the $\lambda\phi^{4}$ field theory by the perturbation of small
coupling constant. We describes the lattice scalar field as coupled
oscillators. In two coupled oscillators, to the $\lambda^{n}$ order, the
square distance of excited state between target and reference state is

$\displaystyle
D^{(n)2}_{(N_{1},N_{2})}=(N_{1}+1)\left(\ln\Big{(}\sqrt{R^{(n)}_{1}}\Big{)}\right)^{2}+(N_{2}+1)\left(\ln\Big{(}\sqrt{R^{(n)}_{2}}\Big{)}\right)^{2}$
(1.4)

in which $R^{(n)}_{1}$ and $R^{(n)}_{2}$ are described by simple recurrent
relations

$\displaystyle R^{(n)}_{1}={\omega_{1}+{3\lambda\over
2}\,\Big{(}{1+N_{1}R^{(n-1)}_{1}\over
2\omega^{2}_{1}}+{1+N_{2}R^{(n-1)}_{2}\over\omega_{1}\omega_{2}}\Big{)}\over\omega_{f}+{3\lambda\over
2}\,{1+N_{1}\over
2\omega^{2}_{f}}},~{}~{}R^{(n)}_{2}={\omega_{2}+{3\lambda\over
2}\Big{(}{1+N_{1}R^{(n-1)}_{1}\over\omega_{2}\omega_{1}}+{1+N_{2}R^{(n-1)}_{2}\over
2\omega^{2}_{2}}\Big{)}\over\omega_{f}+{3\lambda\over 2}\,{1+N_{2}\over
2\omega^{2}_{f}}}\ ~{}~{}~{}$ (1.5)

with initial values $R^{(0)}_{i}={\omega_{i}\over\omega_{f}}$. We had also
generalize it to the case of N coupled oscillators which correspond to the
lattice version of $\lambda\phi^{4}$ theory.

Our method in [12] is very general and can be applied to many-body model in
condense matter. In this paper, along the same method, we will to study the
complexity of Bose Hubbard model [13, 14].

In section II we briefly review the Bose Hubbard model. In section III we
describe the method and use it to study the complexity of two-sites Bose-
Hubbard model. A general perturbative formula of the associated complexity is
presented. We see that the formula of complexity of Bose-Hubbard model is very
similar to that of $\lambda\phi^{4}$ field theory. Use it we evaluate the
complexity of Bose-Hubbard model and find that the complexity is divergent at
a critical value of parameter $g$, which can be identified as the quantum
phase transition point. In section IV we extend the method to the N-sites
closed chain Bose-Hubbard model and any dimensional system. We conclude in
last section.

##  2 Bose Hubbard Model

Quantum Phase transitions is known to occur when a coupling constant $g$ in
the Hamiltonian is varied across some critical value $g=g_{c}$ at which the
ground state dependence on $g$ becomes non-analytic [15, 16]. In this paper we
are interesting in the quantum phase transition happens in a system of Bosons
in the background of a periodic potential and a two-particle repulsive
interaction. The Hamiltonian is given by

$\displaystyle H=\int
d^{d}x\left(-\frac{\hbar^{2}}{2m}\psi^{\dagger}(x)\nabla^{2}\psi(x)\right)+\frac{{\cal
U}_{0}}{2}\int d^{d}x\psi^{\dagger}(x)\psi^{\dagger}(x)\psi(x)\psi(x)$ (2.1)

The constant ${\cal U}_{0}$ is the short-ranged repulsion between the Bosons.
At low energies, we can keep only the lowest vibrational state at each minima
of ${\cal V}_{0}(x)$ and the dynamics is given by the Bose Hubbard model [13]:

$\displaystyle{\cal H}$ $\displaystyle=$ $\displaystyle{J\over N}\sum_{\langle
i,j\rangle}\tilde{a}^{\dagger}_{i}\tilde{a}_{j}+{1\over
N}\sum_{i}\epsilon_{i}\tilde{n}_{i}+{U\over
N}\sum_{i}\tilde{n}_{i}(\tilde{n}_{i}-1),~{}~{}~{}~{}~{}~{}~{}\tilde{n}_{i}=\tilde{a}^{\dagger}_{i}\tilde{a}_{i}~{}~{}$
(2.2)

Here $\langle i,j\rangle$ refers to all nearest neighbor pairs, $J$ is a
hopping term while $U$ is an on-site repulsion proportional to ${\cal U}_{0}$.
The factor ${1\over N}$ is used to render ${\cal H}$ to be energy per site.

We consider the homogeneous case and set $\epsilon_{i}=0$ at each site. This
system has a superfluid-Mott insulator transition first studied by [17]. The
insulating phase is characterized by zero compressibility and has a gap in the
excitation spectrum. The Bose-Hubbard model can describe optical lattices and
its properties, such as the excitation spectrum, had been experimentally
probed [18, 19].

##  3 Complexity of Two-Sites Bose Hubbard Model

Hamiltonian of 2-sites Bose Hubbard Model can be written as

$\displaystyle{\cal H}_{1,2}={J\over
2}(\tilde{a}^{\dagger}_{1}\tilde{a}_{2}+\tilde{a}^{\dagger}_{2}\tilde{a}_{1})+{U\over
2}\left((\tilde{a}_{1}^{\dagger}\tilde{a}_{1})^{2}+(\tilde{a}_{2}^{\dagger}\tilde{a}_{2})^{2}-(\tilde{a}_{1}^{\dagger}\tilde{a}_{1})-(\tilde{a}_{2}^{\dagger}\tilde{a}_{2})\right)$
(3.1)

in which we consider the periodic boundary as shown in figure 1.

Figure 1: Two sites Bose-Hubbard Model model. We choose the model with
periodic boundary and there has two links. The associated complexity is
plotted in figure 2.

Now, we define new operators $a$ and $b$

$\displaystyle
a\equiv{1\over\sqrt{2}}(\tilde{a}_{1}+\tilde{a}_{2}),~{}~{}b\equiv{1\over\sqrt{2}}(\tilde{a}_{1}-\tilde{a}_{2})$
(3.2)

then, by $[\tilde{a}_{i},\tilde{a}_{j}^{\dagger}]=\delta_{i,j}$, we have
commutation relations

$\displaystyle[a,a^{\dagger}]={1\over
2}[\tilde{a}_{1}+\tilde{a}_{2},\tilde{a}^{\dagger}_{1}+\tilde{a}^{\dagger}_{2}]={1\over
2}+{1\over 2}=1$ (3.3) $\displaystyle[b,b^{\dagger}]={1\over
2}[\tilde{a}_{1}-\tilde{a}_{2},\tilde{a}^{\dagger}_{1}-\tilde{a}^{\dagger}_{2}]={1\over
2}+{1\over 2}=1$ (3.4) $\displaystyle[a,b^{\dagger}]={1\over
2}[\tilde{a}_{1}+\tilde{a}_{2},\tilde{a}^{\dagger}_{1}-\tilde{a}^{\dagger}_{2}]={1\over
2}-{1\over 2}=0=[a^{\dagger},b]$ (3.5)
$\displaystyle[a,b]=[a^{\dagger},b^{\dagger}]=0$ (3.6)

Thus $\\{a,a^{\dagger}\\}$ commute with $\\{b,b^{\dagger}\\}$ while each set
has standard commutation relation. We can then define the number operators
$n_{a}$ and $n_{b}$

$\displaystyle n_{a}$ $\displaystyle=$ $\displaystyle
a^{\dagger}a,~{}~{}~{}~{}n_{b}=b^{\dagger}b,~{}~{}~{}~{}~{}[n_{a},n_{b}]=0$
(3.7)

which define the Fock space

$\displaystyle n_{a}|N_{a}\rangle$ $\displaystyle=$ $\displaystyle
N_{a}|N_{a}\rangle,~{}~{}~{}~{}~{}n_{b}|N_{b}\rangle=N_{b}|N_{b}\rangle$ (3.8)

and kinematic term ${\cal H}^{(0)}$ becomes

$\displaystyle{\cal H}^{(0)}_{1,2}$ $\displaystyle=$ $\displaystyle{J\over
2}(\tilde{a}^{\dagger}_{1}\tilde{a}_{2}+\tilde{a}^{\dagger}_{2}\tilde{a}_{1})={J\over
2}(a^{\dagger}a-b^{\dagger}b)$ (3.9)

To proceed we begin to calculated the interaction terms :

$\displaystyle\tilde{a}_{1}^{\dagger}\tilde{a}_{2}+\tilde{a}_{2}^{\dagger}\tilde{a}_{1}$
$\displaystyle=$ $\displaystyle{1\over
2}(a^{\dagger}+b^{\dagger})(a-b)+{1\over
2}(a^{\dagger}-b^{\dagger})(a+b)=a^{\dagger}a-b^{\dagger}b=n_{a}-n_{b}$ (3.10)
$\displaystyle\tilde{a}_{1}^{\dagger}\tilde{a}_{1}+\tilde{a}_{2}^{\dagger}\tilde{a}_{2}$
$\displaystyle=$ $\displaystyle{1\over
2}(a^{\dagger}+b^{\dagger})(a+b)+{1\over
2}(a^{\dagger}-b^{\dagger})(a-b)=a^{\dagger}a+b^{\dagger}b=n_{a}+n_{b}$ (3.11)

and

$\displaystyle(\tilde{a}_{1}^{\dagger}\tilde{a}_{1})^{2}$ $\displaystyle=$
$\displaystyle\Big{(}{1\over
2}(a^{\dagger}+b^{\dagger})(a+b)\Big{)}^{2}={1\over
4}\Big{(}(a^{\dagger}a+b^{\dagger}b)+(a^{\dagger}b+ab^{\dagger})\Big{)}^{2}$
(3.12) $\displaystyle=$ $\displaystyle{1\over
4}\Big{(}n_{a}^{2}+n_{b}^{2}+2n_{a}n_{b}+(a^{\dagger}bab^{\dagger}+ab^{\dagger}a^{\dagger}b)+\text{irrelevant
terms}\Big{)}$ $\displaystyle=$ $\displaystyle{1\over
4}\Big{(}n_{a}^{2}+n_{b}^{2}+2n_{a}n_{b}+(n_{a}(1+n_{b})+(1+n_{a})n_{b})+\text{irrelevant
terms}\Big{)}$ $\displaystyle=$ $\displaystyle{1\over
4}\Big{(}n_{a}^{2}+n_{b}^{2}+4n_{a}n_{b}+n_{a}+n_{b}+\text{irrelevant
terms}\Big{)}$ $\displaystyle(\tilde{a}_{2}^{\dagger}\tilde{a}_{2})^{2}$
$\displaystyle=$ $\displaystyle\Big{(}{1\over
2}(a^{\dagger}-b^{\dagger})(a-b)\Big{)}^{2}={1\over
4}\Big{(}(a^{\dagger}a+b^{\dagger}b)-(a^{\dagger}b+ab^{\dagger})\Big{)}^{2}$
(3.13) $\displaystyle=$ $\displaystyle{1\over
4}\Big{(}n_{a}^{2}+n_{b}^{2}+2n_{a}n_{b}-(a^{\dagger}bab^{\dagger}+ab^{\dagger}a^{\dagger}b)+\text{irrelevant
terms}\Big{)}$ $\displaystyle=$ $\displaystyle{1\over
4}\Big{(}n_{a}^{2}+n_{b}^{2}+2n_{a}n_{b}-(n_{a}(1+n_{b})+(1+n_{a})n_{b})+\text{\it
irrelevant terms}\Big{)}$ $\displaystyle=$ $\displaystyle{1\over
4}\Big{(}n_{a}^{2}+n_{b}^{2}-n_{a}-n_{b}+\text{\it irrelevant terms}\Big{)}$

in which the ”irrelevant terms” are those with different power of $a$ and
$a^{\dagger}$ and/or with different power of $b$ and $b^{\dagger}$.

In this paper we investigate the quantity $\langle N_{a},N_{b}|{\cal
H}_{1.2}|N_{a},N_{b}\rangle$ for the states $|N_{a},N_{b}\rangle$ with fixed
$N_{a}$ and $N_{b}$, and calculate the complexity associated to the states. In
this way, only the terms that have the same power of $a$ and $a^{\dagger}$ and
that have the same power of $b$ and $b^{\dagger}$ are relevant since that
$\langle N_{a},N_{b}|{\it irrelevant\,terms}|N_{a},N_{b}\rangle=0$. Therefore
we neglect the irrelevant terms hereafter.

Collect all then Hamiltonian of 2-sites Bose-Hubbard model becomes

$\displaystyle{\cal H}_{1,2}$ $\displaystyle=$ $\displaystyle{J\over
2}(\tilde{a}^{\dagger}_{1}\tilde{a}_{2}+\tilde{a}^{\dagger}_{2}\tilde{a}_{1})+{U\over
2}\left((\tilde{a}_{1}^{\dagger}\tilde{a}_{1})^{2}-(\tilde{a}_{1}^{\dagger}\tilde{a}_{1})+(\tilde{a}_{2}^{\dagger}\tilde{a}_{2})^{2}-(\tilde{a}_{2}^{\dagger}\tilde{a}_{2})\right)$
(3.14) $\displaystyle=$ $\displaystyle{J\over 2}(n_{a}-n_{b})+{U\over
2}\Big{(}n_{a}^{2}+n_{b}^{2}+2n_{a}n_{b}-n_{a}-n_{b}\Big{)}+\text{\it
irrelevant terms}~{}~{}~{}$

We will use above Hamiltonian to find the complexity of the Bose-Hubbard
model. We first consider the free theory without the interaction and then the
interacting theory.

###  3.1 Complexity of Free Theory

Hamiltonian (3.14) is regarded as the target system which has the kinetic term
and potential term. The target kinetic term is

$\displaystyle{\cal K}^{\rm(tar)}_{a,b}={J\over
2}(\tilde{a}^{\dagger}_{1}\tilde{a}_{2}+\tilde{a}^{\dagger}_{2}\tilde{a}_{1})={J_{a}\over
2}\,a^{\dagger}a-{J_{b}\over 2}\,b^{\dagger}b$ (3.15)

in which, for later convenience, we define

$\displaystyle J_{a}=J_{b}=J$ (3.16)

To calculate the complexity we choose the following kinetic term as that of
the reference state

$\displaystyle K^{\rm(ref)}_{a^{\rm(ref)},\,b^{\rm(ref)}}$ $\displaystyle=$
$\displaystyle{J_{f}\over 2}(a^{\rm(ref)})^{\dagger}\,a^{\rm(ref)}-{J_{f}\over
2}\,(b^{\rm(ref)})^{\dagger}\,b^{\rm(ref)}$ (3.17)

where ${J_{f}}$ is the coupling strength of the reference state.

Now we see that with the replacement

$\displaystyle a^{\rm(ref)}\rightarrow\sqrt{J_{a}\over
J_{f}}\,a,~{}~{}b^{\rm(ref)}\rightarrow\sqrt{J_{b}\over J_{f}}\,b~{}~{}~{}$
(3.18)

one can obtain $K^{\rm(tar)}$ from $K^{\rm(ref)}$, i.e.

$\displaystyle K^{\rm(\rm ref)}\rightarrow K^{\rm(tar)}$ (3.19)

In the operator approach the gate matrix defined in (1.1) is constructed by
the transformation from target operator to reference operator in above
relation.

####  3.1.1 Ground State

Consider first the ground state. Target ground state is annihilated by $a,b$,
i.e. $a|0,0\rangle=b|0,0\rangle=0$, and reference ground state is annihilated
by $a^{\rm(ref)},b^{\rm(ref)}$, i.e. $a^{\rm(ref)}|0,0\rangle_{\rm
ref}=b^{\rm(ref)}|0,0\rangle_{\rm ref}=0$. In the operator approach to
complexity the gate matrix $U(s)$, defined in (1.1), connecting the target
operator to reference operator can be read from (3.18)

$\displaystyle\left(\begin{array}[]{cc}a&b\\\
\end{array}\right)=U(1)\left(\begin{array}[]{c}a^{\rm(ref)}\\\ b^{\rm(ref)}\\\
\end{array}\right)~{}~{}~{},~{}~{}~{}~{}U(1)=\left(\begin{array}[]{cc}\sqrt{J_{a}\over
J_{f}}&0\\\ 0&\sqrt{J_{a}\over J_{f}}\\\ \end{array}\right)$ (3.25)

in which the initial condition is $U(0)=diag(1,1)$. Since the transformation
matrix $U(1)$ is diagonal we can choose ${\cal O_{I}}=1$ in (1.1) and have a
simple relation

$\displaystyle
U(1)=\left(\begin{array}[]{cc}e^{\int_{0}^{1}\,ds\,Y_{a}(s)}&0\\\
0&e^{\int_{0}^{1}\,ds\,Y_{b}(s)}\\\
\end{array}\right)=\left(\begin{array}[]{cc}\sqrt{J_{a}\over J_{f}}&0\\\
0&\sqrt{J_{a}\over J_{f}}\\\ \end{array}\right)$ (3.30)

The associated solutions of $Y_{(1,2)}(s)$ are

$\displaystyle Y_{a}(s)$ $\displaystyle=$
$\displaystyle\ln\left(\sqrt{J_{a}\over
J_{f}}\right),~{}~{}~{}Y_{b}(s)=\ln\left(\sqrt{J_{b}\over
J_{f}}\right)~{}~{}~{}~{}$ (3.31)

which satisfied the initial condition. For the convention used in later we
define new functions

$\displaystyle R^{(0)}_{a}={J_{a}\over J_{f}},~{}~{}~{}R^{(0)}_{b}={J_{b}\over
J_{f}},~{}~{}~{}~{}R^{(0)}={J\over J_{f}}~{}~{}~{}~{}~{}$ (3.32)

in which $(0)$ means a quantity of the zero order of interaction. The squared
distance for ground state, denoted as $(D^{(0)}_{(0,0)})^{2}$, between target
and reference state can be calculated by formula (1.2)

$\displaystyle(D^{(0)}_{(0,0)})^{2}=Y_{a}(1)^{2}+Y_{b}(1)^{2}=\left(\ln\left(R^{(0)}_{a}\right)\right)^{2}+\left(\ln\left(R^{(0)}_{b}\right)\right)^{2}~{}~{}~{}~{}$
(3.33)

in which $(0,0)$ means a quantity of the ground state of two sites system.
This formula has a similar form as that in scalar field theory [6, 11].

####  3.1.2 Excited States

Consider next the $\\{N_{a}^{th},\,N_{b}^{th}\\}$ excited state which is
defined by

$\displaystyle|N_{a},N_{b}\rangle={(a^{\dagger})^{N_{a}+1}(b^{\dagger})^{N_{b}+1}\over\sqrt{(N_{a}+1)!\,(N_{b}+1)!}}|0,0\rangle,~{}~{}~{}a^{N_{a}+1}\,b^{N_{b}+1}|N_{a},N_{b}\rangle=0$
(3.34)

In the operator approach to complexity the gate matrices can be read from the
transformations

$\displaystyle\overbrace{a^{\rm(ref)}\,....a^{\rm(ref)}}^{N_{q}+1}\,\overbrace{b^{\rm(ref)}\,....b^{\rm(ref)}}^{N_{b}+1}\,\,\Rightarrow\,\,\overbrace{\sqrt{J_{a}\over
J_{f}}\,a\,...\sqrt{J_{a}\over
J_{f}}\,a}^{N_{a}+1}\,\overbrace{\sqrt{J_{b}\over
J_{f}}\,b\,....\sqrt{J_{b}\over J_{f}}\,b}^{N_{b}+1}\,$ (3.35)

Then, the gate matrix connecting the target operator with referenct operator
in (1.1) becomes a $(N_{a}+1)\times(N_{b}+1)$ diagonal matrix $U(s)$

$\displaystyle U(1)=\text{Diag}\left(\overbrace{\sqrt{J_{a}\over
J_{f}},...,\sqrt{J_{a}\over J_{f}}}^{N_{a}+1}\,,\overbrace{\sqrt{J_{b}\over
J_{f}},...,\sqrt{J_{b}\over J_{f}}}^{N_{b}+1}\right)$ (3.36)

which redues to (3.25) in the case of ground state $N_{a}=N_{b}=0$.

Follow the discussions in before the gate matrix $U_{i}$ defined in (1.1) now
becomes

$\displaystyle U_{i}(1)$ $\displaystyle=$ $\displaystyle\sqrt{J\over
J_{f}}\,,~{}~{}{\rm with}~{}~{}U_{i}(0)=1,~{}~{}1\leq i\leq N_{2}+N_{1}+2$
(3.37)

after use the relation $J_{a}=J_{b}=J$. The associated functions of $Y_{i}(s)$
in (1.1) are

$\displaystyle Y_{i}(s)$ $\displaystyle=$ $\displaystyle\ln\Big{(}\sqrt{J\over
J_{f}}\Big{)},~{}~{}1\leq i\leq N_{a}+N_{b}+2$ (3.38)

The squared distance for excited state, denoted as $D^{2}_{(N_{a},N_{b})}$,
between target and reference state calculated from (1.2) is

$\displaystyle(D^{(0)}_{(N_{a},N_{b})})^{2}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N_{a}+N_{b}+2}Y_{i}(1)^{2}=(N_{a}+N_{b}+2)\left(\ln\Big{(}\sqrt{R^{(0)}}\Big{)}\right)^{2},~{}~{}~{}R^{(0)}={J\over
J_{f}}~{}~{}$ (3.39)

This has a similar form as the result of scalar theory obtained earlier in
[11].

Since that state wavefunction is described by
$\Psi_{n}(x)={1\over\sqrt{n!}}\langle x|(a^{\dagger})^{n}|0\rangle$ the gate
matrix of wavefunction, $\Psi_{n}(x)$, is thus related to the gate matrix of
field operators, $(a^{\dagger})^{n}$.

###  3.2 Complexity of Interacting Theory

We next study the correction to the complexity due to the interaction term. We
write the Hamiltonian as the summation of kinetic (free) term and potential
(interaction) term

$\displaystyle{\cal H}_{1,2}$ $\displaystyle=$ $\displaystyle{J\over
2}(\tilde{a}^{\dagger}_{1}\tilde{a}_{2}+\tilde{a}^{\dagger}_{2}\tilde{a}_{1})+{U\over
2}\left((\tilde{a}_{1}^{\dagger}\tilde{a}_{1})^{2}-(\tilde{a}_{1}^{\dagger}\tilde{a}_{1})+(\tilde{a}_{2}^{\dagger}\tilde{a}_{2})^{2}-(\tilde{a}_{2}^{\dagger}\tilde{a}_{2})\right)$
(3.40) $\displaystyle=$ $\displaystyle{J\over 2}(n_{a}-n_{b})+{U\over
2}\Big{(}n_{a}^{2}+n_{b}^{2}+2n_{a}n_{b}-n_{a}-n_{b}+\text{irrelevant
terms}\Big{)}$ $\displaystyle=$ $\displaystyle{\cal K}_{a,b}^{\rm(tar)}+{\cal
V}_{a,b}$

The interacting term is

$\displaystyle{\cal V}_{a,b}^{\rm(tar)}$ $\displaystyle=$
$\displaystyle{U\over
2}\Big{(}n_{a}^{2}+n_{b}^{2}+2n_{a}n_{b}-n_{a}-n_{b}\Big{)}~{}~{}~{}$ (3.41)

To consider $\langle N_{a},N_{b}|{\cal V}_{a,b}|N_{a},N_{b}\rangle$ for the
excited state $|N_{a},N_{b}\rangle$ we write above equation as

$\displaystyle{\cal V}_{a,b}^{\rm(tar)}$ $\displaystyle=$
$\displaystyle{U\over
2}\Big{(}N_{a}\,a^{\dagger}a+N_{b}\,b^{\dagger}b+N_{a}\,b^{\dagger}b+N_{b}\,a^{\dagger}a-a^{\dagger}a-b^{\dagger}b\Big{)}$
(3.42)

Therefore

$\displaystyle{\cal H}_{a,b}^{\rm(tar)}={J\over 2}\left(1-{U\over
J}\,(N_{a}+N_{b}-1\Big{)}\right)a^{\dagger}a+{J\over 2}\left(-1-{U\over
J}\,(N_{a}+N_{b}-1\Big{)}\right)b^{\dagger}b~{}~{}$ (3.43)

The associated Hamiltonian of the reference state is

$\displaystyle H^{\rm(ref)}_{a^{\rm(ref)},\,b^{\rm(ref)}}$ $\displaystyle=$
$\displaystyle{J_{f}\over 2}\left(1-{U\over
J_{f}}\,(N_{a}-1\Big{)}\right)(a^{\rm(ref)})^{\dagger}a^{\rm(ref)}+{J_{f}\over
2}\left(-1-{U\over
J_{f}}\,(N_{b}-1\Big{)}\right)(b^{\rm(ref)})^{\dagger}b^{\rm(ref)}~{}~{}$

Above choice satisfies a desirable property of the reference state that it
does not contain any entanglement between operators
$\\{N_{a},\,a^{\rm(ref)},\,(a^{\rm(ref)})^{\dagger}\\}$ and
$\\{N_{b},\,b^{\rm(ref)},\,(b^{\rm(ref)})^{\dagger}\\}$. The property in
operator approach is like as that in coordinate approach in which a desirable
property of the reference state is that it should not contain any entanglement
between the original coordinates $x_{1}$ and $x_{2}$, as discussed in [6, 12].

####  3.2.1 First Order

In the case of zero-order of $U$ the Hamiltonian remain only the kinetic term.
This is the free case discussed in previous section. Now consider the
perturbation to the complexity for the two-sites Bose-Hubbard Model. At the
first order of $U$ we can use (3.43) and (3.2) to find transformations

$\displaystyle\left\\{\begin{array}[]{ccc}J\left(1-{U\over
J}\,(N_{a}+N_{b}-1)\right)a^{\dagger}a&\rightarrow&J_{f}\left(1-{U\over
J_{f}}\,(N_{a}-1)\right)(a^{\rm(ref)})^{\dagger}a^{\rm(ref)}\\\ \\\
J\left(1+{U\over
J}\,(N_{a}+N_{b}-1)\right)b^{\dagger}b&\rightarrow&J_{f}\left(1+{U\over
J_{f}}\,(N_{b}-1)\right)(b^{\rm(ref)})^{\dagger}b^{\rm(ref)}\end{array}\right.~{}~{}~{}$
(3.48)

To proceed we notice that the factors $N_{(a,b)}$ are within the coupling
term, i.e. ${U\over 2J}$ and we only need to consider their zero-order
transform. Thus, by (3.32), we can put following replacement

$\displaystyle N_{a}\rightarrow R^{(0)}_{a}\,N_{a}$ (3.49) $\displaystyle
N_{b}\rightarrow R^{(0)}_{b}\,N_{b}$ (3.50)

into (3.48) and final formulas in first-order transformations are

$\displaystyle\left\\{\begin{array}[]{ccc}R^{(1)}_{a}{(N_{a},\,N_{b})}&=&{J\left(1-{U\over
J}\,\left(R^{(0)}_{a}N_{a}+R^{(0)}_{b}N_{b}-1\right)\right)\over
J_{f}\left(1-{U\over J_{f}}\,(N_{a}-1)\right)}\\\
R^{(1)}_{b}{(N_{a},\,N_{b})}&=&{J\left(1+{U\over
J}\,\left(R^{(0)}_{a}N_{a}+R^{(0)}_{b}N_{b}-1\right)\right)\over
J_{f}\left(1+{U\over
J_{f}}\,(N_{b}-1)\right){}{}{}}\end{array}\right.~{}~{}~{}$ (3.53)

The first-order square distance is

$\displaystyle\Big{(}D^{(1)}_{(N_{a},N_{b})}\Big{)}^{2}=(N_{a}+1)\left(\ln\Big{(}\sqrt{R^{(1)}_{a}(N_{a},N_{b})}\Big{)}\right)^{2}+(N_{b}+1)\left(\ln\Big{(}\sqrt{R^{(1)}_{b}(N_{a},N_{b})}\Big{)}\right)^{2}~{}~{}$
(3.54)

which reduce to the square distance formula of free theory in (3.33) for
ground (i.e.$N_{a}=N_{b}=0$)

####  3.2.2 n’th-Order

Extending to higher-order interactions is straightforward. In terms of new
variables

$\displaystyle g={U\over J},~{}~{}~{}~{}\gamma={J\over J_{f}}~{}~{}~{}$ (3.55)

the recursion relations are

$\displaystyle\left\\{\begin{array}[]{ccc}R^{(n)}_{a}{(N_{a},\,N_{b})}&=&{\gamma\left(1-{g}\,\left(R^{(n-1)}_{a}N_{a}+R^{(n-1)}_{b}N_{b}\,-1\right)\right)\over\left(1-{g\,\gamma}\,(N_{a}-1)\right)},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}R^{(0)}_{a}={\gamma}\\\
R^{(n)}_{b}{(N_{a},\,N_{b})}&=&{\gamma\left(1+{g}\,\left(R^{(n-1)}_{a}N_{a}+R^{(n-1)}_{b}N_{b}\,-1\right)\right)\over
J_{f}\left(1+{g\,\gamma}\,(N_{b}-1)\right)},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}R^{(0)}_{b}={\gamma}\end{array}\right.~{}~{}$
(3.58)

with initial values $R^{(0)}_{(a,b)}$ defined in (3.32) with $J_{a}=J_{b}=J$.
For excited states, the $n$-order square distance is

$\displaystyle\Big{(}D^{(n)}_{(N_{a},N_{b})}\Big{)}^{2}=(N_{a}+1)\left(\ln\Big{(}\sqrt{R^{(n)}_{a}(N_{a},N_{b})}\Big{)}\right)^{2}+(N_{b}+1)\left(\ln\Big{(}\sqrt{R^{(n)}_{b}(N_{a},N_{b})}\Big{)}\right)^{2}~{}~{}$
(3.59)

which is the general formula of $n^{\prime}th$-order complexity of two-sites
Bose-Hubbard model.

###  3.3 Complexity and Quantum Phase Transition

We now use above formula to investigate the quantum phase transition in two-
sites Bose-Hubbard model. Consider first the ground in which $N_{a}=N_{b}=0$.
Then (3.58) become

$\displaystyle R^{(n)}_{a}{(0,\,0)}$ $\displaystyle=$
$\displaystyle{\gamma\left(1+g\right)\over
1+g\,\gamma},~{}~{}~{}~{}~{}~{}~{}R^{(n)}_{b}{(0,\,0)}={\gamma\left(1-g\right)\over
1-g\,\gamma}~{}~{}~{}~{}$ (3.60)

in which $g={U\over J},\,\gamma={J\over J_{f}}$, defined in (3.55), describes
the ratio of one-site repulsion and hopping constant. We expect that
complexity will be divergent at critical value $g_{c}$, which indicates the
quantum phase transition occurs at this point.

To proceed we make following comments about above relation :

1\. $R^{(n)}_{a,b}{(0,\,0)}$ are the qualities at n’th order of interaction
strength g. Thus (3.60) gives

$\displaystyle R^{(0)}_{a,b}{(0,\,0)}={J\over J_{f}}=\gamma$ (3.61)

which are those in (3.32). To first order (3.60) gives

$\displaystyle R^{(1)}_{a}{(0,\,0)}$ $\displaystyle=$
$\displaystyle\gamma+\gamma(1-\gamma)g,~{}~{}~{}~{}~{}~{}~{}R^{(1)}_{b}{(0,\,0)}=\gamma-\gamma(1-\gamma)g$
(3.62)

and using (3.60) we can find any order of $R^{(n)}_{a,b}{(0,\,0)}$ by series
expansion of $g$.

2\. In fact the formulas (3.60) are the exact values and we can let
$n\to\infty$. Note that being able to obtain the exact formula for ground
state complexity is a special property in our method. Now we see that
$R^{(n\to\infty)}_{b}{(0,\,0)}$ in (3.60) becomes zero at critical value
$g=g_{c}=1$. This implies that the associated complexity,
$(D^{(n\to\infty)}_{(0,0)})^{2}$, is divergent and thus it is the quantum
phase transition point.

3\. We plot the figure 2 to explicitly see this property.

Figure 2: Dependence of complexity $(D^{(n\to\infty)}_{(0,0)})^{2}$ on $``g"$
in 2-sites Bose-Hubbard model. The divergence at $g=g_{c}=1$ indicates the
quantum phase transition at this point. The figure is plotted with
$\gamma_{L}=1/2$ while $\gamma_{R}=3/2$.

4\. The value of $g$ we considered is positive. Thus the factor
$(1-g\,\gamma)$ in formulas (3.60) could become zero at $g_{0}=\gamma^{-1}$
which depends on the reference state parameter $\gamma$. To avoid $"g_{0}"$
become a critical point the following scheme is proposed. For range of
$g>g_{c}$ we choose $\gamma=\gamma_{R}>1$ while for range of $g<g_{c}$ we
choose $\gamma=\gamma_{L}<1$. Under this prescription the value of
$(1-g\,\gamma)$ is always positive and it remains only one critical point at
$g=g_{c}$.

5\. Notice that, in this way, the complexity calculated for $g>g_{c}$ is that
with $\gamma=\gamma_{R}>1$ while the complexity calculated with $g<g_{c}$ is
that with $\gamma=\gamma_{L}<1$. Therefore, one could compare the magnitude of
complexity between the systems of $g>g_{c}$. One could also compare the
magnitude of complexity between the systems of $g<g_{c}$. However, it makes no
sense to compare the magnitude of complexity between the systems of $g<g_{c}$
and of systems $g>g_{c}$ since the reference state parameter $\gamma$ is
different between them. The figure 2 is plotted with $\gamma_{L}=1/2$ while
$\gamma_{R}=3/2$.

6\. Consider next the excited ground in which $N_{a}=N_{b}=1$. For $n\geq 1$
(3.58) becomes

$\displaystyle R^{(n)}_{a}{(1,\,1)}$ $\displaystyle=$
$\displaystyle\gamma\left(1+g(1-2\gamma)\right),~{}~{}~{}~{}~{}~{}~{}R^{(n)}_{b}{(1,\,1)}=\gamma\left(1-g(1-2\gamma)\right),~{}~{}~{}n\geq
1~{}~{}~{}~{}$ (3.63)

As that in ground case, the formulas (3.63) are the exact values and we can
let $n\to\infty$. We also see that either $R^{(n\to\infty)}_{a}{(1,\,1)}$ or
$R^{(n\to\infty)}_{a}{(1,\,1)}$ becomes zero at critical value
$g=g^{(N=1)}_{c}={\pm 1\over 1-2\gamma}$. This implies that the associated
complexity, $(D^{(n\to\infty)}_{(1,1)})^{2}$, is divergent at the critical
value and thus it is the quantum phase transition point. Note that the
critical point is a physical property and it shall not depend on the reference
state parameter $\gamma$. However, the calculated value $g^{(N=1)}_{c}$
depends on $\gamma$. This feels short of author expectation. The reason may be
that the model described in (2.2) is that keeps only the lowest vibrational
state and could only describe the ground state property. The analysis in
above, however, have shown a possible way to calculate excited-state
complexity in other many-body models.

##  4 Complexity of N-Sites Bose-Hubbard Model

Hamiltonian of N-sites closed chain of Bose-Hubbard model defined in (2.2) can
be written as

$\displaystyle{\cal H}_{N}={J\over
2N}\sum_{j=0}^{N}(\tilde{a}^{\dagger}_{j+1}\tilde{a}_{j}+\tilde{a}^{\dagger}_{j-1}\tilde{a}_{j})+\frac{U}{N}\sum_{j=0}^{N}\tilde{n}_{j}(\tilde{n}_{j}-1),~{}~{}~{}~{}~{}~{}~{}\tilde{n}_{j}=\tilde{a}^{\dagger}_{j}\tilde{a}_{j}~{}~{}$
(4.1)

in which we consider the periodic boundary as shown in figure 3.

Figure 3: Six-sites Bose-Hubbard Model model. We choose the model with
periodic boundary and there has six links. The associated complexity is
plotted in figure 4.

We can generalize method in previous section to define the new operators

$\displaystyle a_{k}$ $\displaystyle=$
$\displaystyle{1\over\sqrt{N}}\sum_{j=0}^{N-1}\,exp\Big{(}{2\pi ik\over
N}\,j\Big{)}\,\tilde{a}_{j},~{}~{}a^{\dagger}_{k}={1\over\sqrt{N}}\sum_{j=0}^{N-1}\,exp\Big{(}{-2\pi
ik\over N}\,j\Big{)}\,\tilde{a}^{\dagger}_{j}$ (4.2)
$\displaystyle\tilde{a}_{k}$ $\displaystyle=$
$\displaystyle{1\over\sqrt{N}}\sum_{j=0}^{N-1}\,exp\Big{(}{-2\pi ik\over
N}\,j\Big{)}\,a_{j},~{}~{}\tilde{a}^{\dagger}_{k}={1\over\sqrt{N}}\sum_{j=0}^{N-1}\,exp\Big{(}{2\pi
ik\over N}\,j\Big{)}\,a^{\dagger}_{j}$ (4.3)

in which we impose a periodic boundary condition $\tilde{a}_{\rm
k+N+1}=\tilde{a}_{k}$ and $\tilde{a}^{\dagger}_{\rm
k+N+1}=\tilde{a}^{\dagger}_{k}$.

Note that the relative sign between the Fourier series of $a_{k}$ and
$a^{\dagger}_{k}$ is important to have standard commuation relation [6]. Then,
above definitions leads to simple relations

$\displaystyle~{}[\tilde{a}_{i},\tilde{a}_{j}]=[\tilde{a}^{\dagger}_{i},\tilde{a}^{\dagger}_{j}]=0,~{}~{}[\tilde{a}_{i},\tilde{a}^{\dagger}_{j}]=\delta_{i,j}~{}~{}~{}\Rightarrow~{}~{}~{}[a_{i},a_{j}]=[a^{\dagger}_{i},a^{\dagger}_{j}]=0,~{}~{}[a_{i},a^{\dagger}_{j}]=\delta_{i,j}$
(4.4)

and

$\displaystyle[a_{j},a^{\dagger}_{k}]$ $\displaystyle=$ $\displaystyle{1\over
N}\sum_{m,\ell}\,e^{{2\pi\,i(m\,j-k\,\ell)\over
N}}\,[a_{m},a^{\dagger}_{\ell}]={1\over
N}\sum_{m,\ell}\,e^{{2\pi\,i(m\,j-k\,\ell)\over N}}\,\delta_{m,\ell}={1\over
N}\sum_{\ell}\,e^{{2\pi\,i\ell(m-k)\over N}}=\delta_{j,k}~{}~{}~{}~{}~{}$
(4.5)

The general states can then be constructed by creation operators
$a_{k}^{\dagger}$

$\displaystyle|{N}_{0},\,{N}_{1}\cdot\,\cdot\,\cdot\,{N}_{N-1}\rangle=a_{0}^{{N}_{0}}\,a_{1}^{{N}_{1}}\cdot\,\cdot\,\cdot\,a_{N-1}^{{N}_{N-1}}\,|0\rangle$
(4.6)

due to the commutation relation between the new operators.

We now express N-sites Hamiltonian ${\cal H}_{N}$ in terms of new operators
$a_{j}$ and $a^{\dagger}_{j}$. Kinematic term is

$\displaystyle\sum_{j}\tilde{a}^{\dagger}_{j+1}\tilde{a}_{j}+\sum_{j}\tilde{a}^{\dagger}_{j-1}\tilde{a}_{j}$
$\displaystyle=$ $\displaystyle{1\over
N}\sum_{j}\sum_{m,k}e^{{2\pi\,i(-m(j+1)+kj)\over
N}}\,a^{\dagger}_{m}a_{k}+{1\over
N}\sum_{j}\sum_{m,k}e^{{2\pi\,i(-m(j-1)+kj))\over N}}\,a^{\dagger}_{m}a_{k}$
(4.7) $\displaystyle=$ $\displaystyle\sum_{k}e^{{-2\pi\,ik\over
N}}\,a^{\dagger}_{k}a_{k}+\sum_{k}e^{{2\pi\,im\over N}}\,a^{\dagger}_{k}a_{k}$
$\displaystyle=$ $\displaystyle 2\sum_{k=0}^{N-1}\,\cos\left({2\pi k\over
N}\right)\,a^{\dagger}_{k}a_{k}$ (4.8)

The kinematic term in ${\cal H}_{N}$ can thus be expressed as

$\displaystyle{\cal K}_{N}={J\over
2N}\left(\sum_{j}\tilde{a}^{\dagger}_{j+1}\tilde{a}_{j}+\sum_{j}\tilde{a}^{\dagger}_{j-1}\tilde{a}_{j}\right)$
$\displaystyle=$ $\displaystyle{J\over N}\sum_{k=0}^{N-1}\cos\left({2\pi
k\over N}\right)\,a^{\dagger}_{k}a_{k}$ (4.9)

which reduces to the formula of N=2 case in previous section. The interaction
terms has two parts. One is

$\displaystyle\sum_{j}\tilde{n}_{j}=\sum_{j}\tilde{a}^{\dagger}_{j}\tilde{a}_{j}={1\over
N}\sum_{j}\sum_{m,k}e^{{2\pi\,i(-m+k)j\over
N}}\,a^{\dagger}_{m}a_{k}=\sum_{m,k}\delta_{m,k}\,a^{\dagger}_{m}a_{k}=\sum_{k}n_{k}$
(4.10)

Another one is

$\displaystyle\sum_{j}\tilde{n}_{j}^{2}$ $\displaystyle=$
$\displaystyle\sum_{j}\tilde{a}^{\dagger}_{j}\tilde{a}_{j}\tilde{a}^{\dagger}_{j}\tilde{a}_{j}={1\over
N^{2}}\sum_{j}\sum_{\alpha,\beta,\gamma,\delta}e^{2\pi
j(-\alpha+\beta-\gamma+\eta)\over
N}a^{\dagger}_{\alpha}a_{\beta}a^{\dagger}_{\gamma}a_{\eta}$ (4.11)
$\displaystyle=$ $\displaystyle{1\over
N}\sum_{\alpha,\beta,\gamma,\eta}a^{\dagger}_{\alpha}a_{\beta}a^{\dagger}_{\gamma}a_{\eta}\,\delta_{\alpha-\beta+\gamma-\eta,0}$
(4.12)

The delta function in here means that once $\alpha=\beta$ then $\gamma=\eta$
and, if $\alpha\neq\beta$ then $\gamma\neq\eta$. In the first case we have
term

$\displaystyle{1\over
N}\sum_{\alpha,\gamma}a^{\dagger}_{\alpha}a_{\alpha}a^{\dagger}_{\gamma}a_{\gamma}={1\over
N}\sum_{\alpha,\gamma}n_{\alpha}\,a^{\dagger}_{\gamma}a_{\gamma}$ (4.13)

In the second case, following the discuss in previous section we have to
choose the relevant term which is

$\displaystyle{1\over
N}\sum_{\alpha\neq\gamma}a^{\dagger}_{\alpha}a_{\gamma}a^{\dagger}_{\gamma}a_{\alpha}$
$\displaystyle=$ $\displaystyle{1\over
N}\sum_{\alpha\neq\gamma}n_{\alpha}\,a_{\gamma}a^{\dagger}_{\gamma}={1\over
N}\sum_{\alpha\neq\gamma}n_{\alpha}\,(1+a^{\dagger}_{\gamma}a_{\gamma})$
(4.14) $\displaystyle=$ $\displaystyle\left({N-1\over
N}\sum_{\alpha}n_{\alpha}\right)+{1\over
N}\sum_{\alpha\neq\gamma}n_{\alpha}a^{\dagger}_{\gamma}a_{\gamma}$

Collect all and the interaction term in ${\cal H}_{N}$ becomes

$\displaystyle{\cal V}_{N}$ $\displaystyle=$ $\displaystyle{U\over
N}\Big{(}\sum_{\alpha,\gamma}n_{\alpha}a^{\dagger}_{\gamma}a_{\gamma}+\sum_{\alpha\neq\gamma}n_{\alpha}a^{\dagger}_{\gamma}a_{\gamma}-\sum_{\alpha}n_{\alpha}\Big{)}+\text{\it
irrelevant terms}$ (4.15) $\displaystyle=$ $\displaystyle{U\over
N}\Big{(}\sum_{k=0}^{N-1}n_{k}a^{\dagger}_{k}a_{k}+2\sum_{(k,j)=0}^{N-1}n_{j}a^{\dagger}_{k}a_{k}-\sum_{k=0}^{N-1}n_{k}\Big{)}+\text{\it
irrelevant terms}$

where $(k,j)$ denotes any pair at site $k$ and site $j$. When N=2 above
formula reduces to that of 2-sites formula in (3.41).

###  4.1 Complexity of Free Theory

Let us first discuss the case of free theory. Following the scheme detailed in
previous section the target system kinematic term and reference system
kinematic term are

$\displaystyle{\cal K}^{\rm(tar)}_{N}$ $\displaystyle=$ $\displaystyle{J\over
N}\sum_{k=0}^{N-1}\cos\left({2\pi k\over
N}\right)\,a^{\dagger}_{k}\,a_{k}~{}~{}$ (4.16) $\displaystyle{\cal
K}^{\rm(ref)}_{N}$ $\displaystyle=$ $\displaystyle{J_{f}\over
N}\sum_{k=0}^{N-1}\cos\left({2\pi k\over
N}\right)a_{k}^{\rm(ref){\dagger}}\,a_{k}^{\rm(ref)}~{}~{}$ (4.17)

where $J_{f}$ is the coupling strength of the reference state.

Now we see that with the replacement

$\displaystyle a^{\rm(ref)}_{k}\rightarrow\sqrt{J\over J_{f}}\,a_{k}$ (4.18)

one can obtain $K^{\rm(tar)}$ from $K^{\rm(ref)}$, i.e.

$\displaystyle K^{\rm(\rm ref)}\rightarrow K^{\rm(tar)}$ (4.19)

As mentioned before, in the operator approach the gate matrix defined in (1.1)
is constructed by the transformation from target operator to reference
operator in above relation. Thus the complexity, square distance, for the
general excited state $\\{{N}_{0},{N}_{1}\cdot\cdot\cdot,{N}_{N-1}\\}$ is easy
to be found by the method discussed in the 2-site case. The result is

$\displaystyle D^{2}_{\\{{N}_{0},{N}_{1}\cdot\cdot\cdot,{N}_{N-1}\\}}$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{N-1}\,({N}_{k}+1)\left(\ln\Big{(}\sqrt{R^{(0)}_{k}}\Big{)}\right)^{2},~{}~{}~{}~{}~{}~{}R^{(0)}_{k}={J\over
J_{f}}~{}~{}~{}$ (4.20)

which is a simple extension of 2-sites formula in (3.39).

###  4.2 Complexity of Interacting Theory

To include the interaction we first collect the total Hamiltonian

$\displaystyle{\cal H}_{N}$ $\displaystyle=$ $\displaystyle{1\over
N}\sum_{k=0}^{N-1}J\cos\left({2\pi k\over
N}\right)\,a^{\dagger}_{k}a_{k}+{U\over
N}\Big{(}\sum_{k=0}^{N-1}n_{k}a^{\dagger}_{k}a_{k}+2\sum_{(k,j)=0}^{N-1}n_{j}a^{\dagger}_{k}a_{k}-\sum_{k=0}^{N-1}a^{\dagger}_{k}a_{k}\Big{)}$
(4.21) $\displaystyle=$ $\displaystyle{1\over
N}\sum_{k=0}^{N-1}\left(J\cos\left({2\pi k\over
N}\right)\,+{U}\Big{(}-1+n_{k}+\sum_{j\neq
k}^{N-1}n_{j}\Big{)}\right)a^{\dagger}_{k}a_{k}$ (4.22)

where $n_{k}=a^{\dagger}_{k}a_{k}$ is the number operator. Then, following the
scheme described in previous section we furthermore write the total
Hamiltonian of target state and reference state as

$\displaystyle{\cal H}^{\rm(tar)}_{N}$ $\displaystyle=$ $\displaystyle{1\over
N}\sum_{k=0}^{N-1}\left(J\cos\left({2\pi k\over
N}\right)\,+{U}\Big{(}-1+{N}_{k}+\sum_{j\neq
k}^{N-1}{N}_{j}\Big{)}\right)a^{\dagger}_{k}a_{k}$ (4.23) $\displaystyle{\cal
H}^{\rm(ref)}_{N}$ $\displaystyle=$ $\displaystyle{1\over
N}\sum_{k=0}^{N-1}\left(J_{f}\cos\left({2\pi k\over
N}\right)\,+{U}\left(-1+{N}_{k}\right)\right)a_{k}^{\rm(ref){\dagger}}\,a_{k}^{\rm(ref)}$
(4.24)

where ${N}_{k}$ is an integral number defined by
$n_{k}|0\rangle=|{N}_{k}\rangle$. The reference operator in above satisfies a
desirable property of the reference state that it does not contain any
entanglement between operators $a_{k}^{\rm(ref)}$ and $a_{j}^{\rm(ref)}$ if
$k\neq j$, as discussed in [6, 12].

####  4.2.1 First Order

In the case of zero-order of $U$ the Hamiltonian remain only the kinetic term.
This is the free case. Now consider the perturbation to the complexity for the
N sites Bose Hubbard Model. At the first order of $U$ the transformation from
target operator to reference operator is chosen as following

$\displaystyle\left(J\cos\left({2\pi k\over
N}\right)\,+{U}\Big{(}-1+{N}_{k}+\sum_{j\neq
k}^{N-1}{N}_{j}\Big{)}\right)a^{\dagger}_{k}a_{k}$ $\displaystyle\rightarrow$
$\displaystyle\left(J_{f}\cos\left({2\pi k\over
N}\right)\,+{U}\left(-1+{N}_{k}\right)\right)a_{k}^{\rm(ref){\dagger}}\,a_{k}^{\rm(ref)}~{}~{}$

Now, we notice that the factors $N_{k},\,N_{j}$ are within the coupling term,
i.e. ${U}$, and we only need to consider their zero-order transform. Thus, by
(4.20), we shall put following replacement

$\displaystyle N_{k}\rightarrow R^{(0)}_{k}\,N_{k}$ (4.26)

into (4.2.1) and final formulas in first-order transformation is

$\displaystyle R^{(1)}_{k}({\\{N_{i}\\}})$ $\displaystyle=$
$\displaystyle{J\cos\left({2\pi k\over
N}\right)\,+{U}\Big{(}-1+R^{(0)}_{k}\,{N}_{k}+\sum_{j\neq
k}^{N-1}{N}_{j}\Big{)}\over J_{f}\cos\left({2\pi k\over
N}\right)\,+{U}\left(-1+{N}_{k}\right)}$ (4.27)

Therefore the first-order complexity or square distance is

$\displaystyle\Big{(}D^{(1)}_{\\{N_{i}\\}}\Big{)}^{2}=\sum_{k=0}^{N_{1}}(N_{k}+1)\,\left(\ln\Big{(}\sqrt{R^{(1)}_{k}({\\{N_{i}\\}})}\Big{)}\right)^{2}$
(4.28)

The case of ${\\{N_{i}\\}}=\\{N_{a},N_{b}\\}$ above formula reduce to the
square distance formula of 2-site theory in (3.54).

####  4.2.2 n’th-Order

Extending to higher-order interactions is straightforward. The recursion
relation is

$\displaystyle R^{(n)}_{k}({\\{N_{i}\\}})$ $\displaystyle=$
$\displaystyle{J\cos\left({2\pi k\over
N}\right)\,+{U}\Big{(}-1+R^{(n-1)}_{k}\,{N}_{k}+\sum_{j\neq
k}^{N-1}{N}_{j}\Big{)}\over J_{f}\cos\left({2\pi k\over
N}\right)\,+{U}\left(-1+{N}_{k}\right)}$ (4.29) $\displaystyle=$
$\displaystyle{\gamma\left(\cos\left({2\pi k\over
N}\right)\,+{g}\left(-1+R^{(n-1)}_{k}\,{N}_{k}+\sum_{j\neq
k}^{N-1}{N}_{j}\right)\right)\over\cos\left({2\pi k\over
N}\right)\,+{\gamma\,g}\left(-1+{N}_{k}\right)}~{}~{}$ (4.30)

with initial values $R^{(0)}_{(a,b)}$ defined in (4.20) and the $n$-order
square distance is

$\displaystyle\Big{(}D^{(n)}_{\\{N_{i}\\}}\Big{)}^{2}=\sum_{k=0}^{N_{1}}(N_{k}+1)\,\left(\ln\Big{(}\sqrt{R^{(n)}_{k}({\\{N_{i}\\}})}\Big{)}\right)^{2}~{}~{}$
(4.31)

These complete our derivations.

###  4.3 Complexity and Quantum Phase Transition

Consider first the ground in which $N_{i}=0$. Then (4.30) become

$\displaystyle R^{(n)}_{k}({\\{0\\}})$ $\displaystyle=$
$\displaystyle{\gamma\,\left(cos\left({2\pi\,k\over N}\right)-{g}\right)\over
cos\left({2\pi\,k\over
N}\right)-{\gamma\,g}}\,,~{}~{}~{}~{}k=0,\cdot\cdot\cdot\cdot,\,N-1~{}~{}~{}~{}$
(4.32)

which reduces to formulas (3.60) for N=2 case. Likes as that in two-sites
model we expect that complexity will be divergent at critical value $g_{c}$,
which indicates the quantum phase transition occurs at this point. Let us
discuss the property of (4.32).

1\. As that in 2-sites system, above formula is the exact value and we can let
$n\to\infty$. Being able to obtain the exact formula for ground state
complexity is a special property in our method.

2\. We see that $R^{(n\to\infty)}_{k=0}({\\{0\\}})$ in (4.32) becomes zero at
critical value $g=g_{c}=1$. This implies that the associated complexity,
$(D^{(n\to\infty)}_{(0,0)})^{2}$, is divergent and thus it is the quantum
phase transition point. The critical value $g_{c}=1$ is same as that in two-
sites system.

3\. We plot the figure 4 to explicitly see this property.

Figure 4: Dependence of complexity $(D^{(n\to\infty)}_{(0,0)})^{2}$ on $``g"$
in 6-sites Bose-Hubbard model. The divergence at $g=g_{c}=1$ indicates the
quantum phase transition at this point. The figure is plotted with
$\gamma_{L}=3/4$ while $\gamma_{R}=2$.

4\. In two-sites system, for range of $g>g_{c}$ we choose
$\gamma=\gamma_{R}>1$ while for range of $g<g_{c}$ we choose
$\gamma=\gamma_{L}<1$. In this algorithm, we see that the complexity becomes
infinity only at $g=g_{c}$. In the N-sites system, however, (4.32) tells us
that complexity could become infinity at the point $g=g_{k}<1$ if
$cos\left({2\pi\,k\over N}\right)-{g_{k}}=0$. While we know that these are not
physical critical points we have not yet found a simple algorithm to rule it
out. Figure 4 is plotted for 6-sites Bose-Hubbard model with $\gamma=2$. The
complexity for $g<1$ has three divergence points which depend on $\gamma$ and
are un-physical.

5\. When consider the excited state in which $\\{N\\}\neq\\{0\\}$ we can see
that the associated complexity become divergent at the critical value which
depends on the reference state parameter $\gamma$. As that in 2-sites system
the un-physical property reason may be that the model described in (2.2) is
that keeps only the lowest vibrational state and could only describe the
ground state property.

6\. Generalize the method to N dimensional Bose-Hubbard model is
straightforward, as that in quantum field theory [6]. We only need to extend
the N-sites closed chain Hamiltonian (4.1) to

$\displaystyle{\cal H}_{N}={J\over
2}\sum_{\vec{j}}(\tilde{a}^{\dagger}_{\vec{j}+\vec{1}}\,\tilde{a}_{\vec{j}}+\tilde{a}^{\dagger}_{{\vec{j}}-{\vec{1}}}\tilde{a}_{\vec{j}})+{U}\sum_{{\vec{j}}}\tilde{n}_{\vec{j}}\,(\tilde{n}_{\vec{j}}-1),~{}~{}~{}~{}~{}~{}~{}\tilde{n}_{\vec{j}}=\tilde{a}^{\dagger}_{\vec{j}}\,\tilde{a}_{\vec{j}}$
(4.33)

in which ${\vec{j}}=(j_{1},\,j_{2},\,\cdot\cdot\cdot,j_{N})$ is the
N-dimensional vector with integral number. In this way, to find the formula of
complexity we only need to extend the single parameter Fourier transform to N
parameters transform and can get a general formula of complexity. For example
the relation (4.31) becomes

$\displaystyle\Big{(}D^{(n)}_{\\{N_{i}\\}}\Big{)}^{2}=\sum_{\vec{k}}(N_{\vec{k}}+1)\,\left(\ln\Big{(}\sqrt{R^{(n)}_{\vec{k}}({\\{N_{i}\\}})}\Big{)}\right)^{2},~{}~{}~{}~{}{\vec{k}}=(k_{1},\,k_{2},\,\cdot\cdot\cdot,k_{N})$
(4.34)

which is the complexity of N dimensional Bose-Hubbard model.

##  5 Concluding Remarks

We adopt the operator approach to compute the complexity of the Bose-Hubbard
Model. Our prescription is to consider the value of $\langle
N_{1}\cdot\cdot\cdot N_{n}|{\cal H}|N_{1}\cdot\cdot\cdot N_{n}\rangle$. In
this way, after proper mode expansion of ${\cal H}$, many irrelevant terms
become zero and the contributions of complexity from the remained terms could
be calculated exactly through a perturbation algorithm. We then obtain the
general formula of complexity therein. It is interesting to see that the
formula of complexity of Bose-Hubbard model is very similar to that of
$\lambda\phi^{4}$ field theory. Our analysis leads to a result that the
critical point $g_{c}=1$ is found in 2-sites model and N-sites system. The
special property that the critical value being independent of lattice number
tells that our method is mean-field like trick. The more general method to
compute complexity in many-body models remained to be found.

Finally, our algorithm can be applied to many quantum field theories and
several many-body models in condense matter. We will study the complexity of
Sachdev-Ye-Kitaev (SYK) model [20, 21, 22, 23] in the next paper [24] to see
how the complexity probes the topological phase transition [25].

References

## References

  * [1] T. Faulkner, M. Guica, T. Hartman, R. C. Myers, and M. Van Raamsdonk, “Gravitation from Entanglement in Holographic CFTs,” JHEP 03 (2014) 051, arXiv:1312.7856 [hep-th].
  * [2] T. Hartman and J. Maldacena, “Time evolution of entanglement entropy from black hole interiors,” JHEP 1305 (2013) 14, arXiv:1303.1080 [hep-th].
  * [3] L. Susskind, “Computational Complexity and Black Hole Horizons,” Fortsch. Phys. 64 (2016) 24, arXiv:1403.5695 [hep-th].
  * [4] A. R. Brown, D. A. Roberts, L. Susskind, B. Swingle, and Y. Zhao, “Holographic Complexity Equals Bulk Action? ” Phys. Rev. Lett. 116 (2016) 191301, arXiv:1509.07876 [hep-th].
  * [5] S. Chapman, H. Marrochio, and R. C. Myers, “Complexity of formation in holography, “ JHEP1701 (2017) 62, arXiv:1610.08063 [hep-th].
  * [6] R. A. Jefferson and R. C. Myers, “Circuit complexity in quantum field theory,” JHEP 10 (2017) 107 arXiv:1707.08570 [hep-th].
  * [7] S. Chapman, M. P. Heller, H. Marrochio, F. Pastawski, “Towards Complexity for Quantum Field Theory States,” Phys. Rev. Lett. 120 (2018) 121602 arXiv:1707.08582 [hep-th].
  * [8] R. Khan, C. Krishnan, and S. Sharma,“ Circuit Complexity in Fermionic Field Theory,” Phys. Rev. D 98 (2018) 126001 , arXiv:1801.07620 [hep-th].
  * [9] L. Hackl and R. C. Myers, “Circuit complexity for free fermions,” JHEP07(2018)139, arXiv:1803.10638 [hep-th].
  * [10] A. Bhattacharyya, A. Shekar, A. Sinha, “Circuit complexity in interacting QFTs and RG flows,” JHEP 1810 (2018) 140, arXiv:1808.03105 [hep-th].
  * [11] W.-H. Huang, “Operator Approach to Complexity : Excited States,” Phys. Rev. D 100 (2019) 066013, arXiv:1905.02041 [hep-th].
  * [12] W.-H. Huang, “Perturbative Complexity of Interacting Theory,” Phys. Rev. D 103 (2021) 065002, arXiv:2008.05944 [hep-th].
  * [13] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller,” Cold Boseic Atoms in Optical Lattices,” Phys.Rev. Lett. 81 (1998) 3108, arXiv:9805329 [cond-mat].
  * [14] U. Sood, M. Kruczenski,”Circuit complexity near critical points, arXiv:2106.12648 [quant-ph].
  * [15] S. Sachdev, “Quantum Phase Transitions,” Cambridge: Cambridge University Press (1999).
  * [16] V. Matthias, “Quantum Phase Transitions,” Reports on Progress in Physics 66 (2003) 2069, arXiv:0309604 [cond-mat].
  * [17] M. P. A. Fisher et al., Phys. Rev. B 40, 546 (1989).
  * [18] M. Greiner, O.Mandel, T.Esslinger, T.W.Hansch, I.Bloch,”Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms,” Nature 415 (2002) 39.
  * [19] T. Stoferle, H.Moritz, C.Schori, M.Kohl, T.Esslinger,”Transition from a Strongly Interacting 1D Superfluid to a Mott Insulator,” Phys.Rev.Lett. 92 (2004) , arXiv:0309604130403, 0312440 [cond-mat].
  * [20] S. Sachdev and J. Ye, Gapless spin fluid ground state in a random, quantum Heisenberg magnet, Phys. Rev. Lett. 70 (1993) 3339 arXiv:9212030 [cond-mat].
  * [21] A. Kitaev, “Talks given at the Fundamental Physics Prize Symposium and KITP seminars.” http://online.kitp.ucsb.edu/online/joint98/kitaev/
http://online.kitp.ucsb.edu/online/entangled15/kitae.

  * [22] J. Polchinski and V. Rosenhaus, The Spectrum in the Sachdev-Ye-Kitaev Model, JHEP 04 (2016) 001 arXiv:1601.06768 [hep-th].
  * [23] J. Maldacena and D. Stanford, Remarks on the Sachdev-Ye-Kitaev model, Phys. Rev. D 94 (2016) arXiv:106002 [hep-th].
  * [24] W.-H. Huang, “ Operator approach to complexity of SYK model”
  * [25] F. Liu, S. Whitsitt, J. B. Curtis, R. Lundgren, P. Titum, Z-C Yang, J. R. Garrison, A. V. Gorshkov, “Circuit Complexity across a Topological Phase Transition,” Phys. Rev. Research 2 (2020) 013323. arXiv:1902.10720 [quant-ph]


# Solving the Stock Option Forecast problem by a numerical method for the
Black-Scholes Equation with Machine Learning Classification Model

Benjamin Jiang, Matthieu Durieux and Kirill V. Golubnichiy  Department of
Mathematics  University of Washington, Seattle, USA<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We proposed classification models that utilize the result from the Quasi-
Reversibility Method[2], which solves the Black-Scholes equation to forecast
the option prices one day in advance. Combining the minimizer from QRM with
our machine learning classifications, we can classify the option as an
increase or decrease in value. Based on the different classifications of the
options, we can apply various trading strategies which we aim to figure out
ways to improve the results from QRM’s extrapolations. To further test the
viability of our model, we collected 23548 options data from the real-world
market for our model, and we will then feed in the data along with the
minimizer from QRM to form decision trees and random forests, which we will
later test for accuracy, precision, and recall.

##  1 Introduction and Overview

The stock market, also known as an equity market, contains buyer and seller
that trades the ownership of a business. In order to make a profit in such a
market, we need to make intuitive guesses about the market trend; in other
words, we want to buy low and sell it high. However, since the market is
volatile and we are solving ill-posed problems because the data are very noisy
and the model is very sensitive to the noise in the data, it is very hard to
correctly predict the exact price the stock will be even after just a day.
Thus, we can make a more generalized prediction, such as estimating a market
gain or loss. On the other hand, options trading is the contract that gives
you the right to sell or buy the stock in the future. If we can predict if the
option is going to trade higher, we can make a profit out of the market.
Hence, we came up with QRM(Quasi Reversibility Method) to capture the market’s
movement by solving an ill-posed problem(Black–Scholes model) to predict the
option price one day in advance.[2] In our research, we will utilize a
minimizer which is the prediction of future option data from the QRM along
with other features from the stock market to forecast the market movement. We
proposed using ensemble methods such as gradient boosting trees and random
forests to classify the future data. We use the classification model because
we believe that by combining it with QRM, we can capture some critical
features from the market that determine the movement of an option’s prices
with the assumption that the market will not be too volatile in a short period
which allows us to split the stock into different categories. We will first
design weak learners to classify the stock into different classes. Then, we
will improve the result by minimizing the loss function.

##  2 Quasi-reversibility model and Black–Scholes model

###  2.1 Black-Scholes equation

The Black-Scholes equation is a backward parabolic partial differential
equation that is used to solve questions in mathematical finances. The
proposing model is to use some initial boundary conditions for stock that we
are interested in and simulate a more accurate options trading strategy. Since
it is an ill-posed problem, we will introduce variables $\tau$ to solve the
equation backward, which at any given time t, we have

$\displaystyle\tau=T-t.$ (2.1)

Give the Black-Scholes equation to be:

$\displaystyle\frac{\partial{V}}{\partial{t}}+\frac{1}{2}\sigma^{2}S^{2}\frac{\partial^{2}{V}}{\partial{S^{2}}}+rS\frac{\partial{V}}{\partial{S}}-rV=0$
(2.2)

Assume that the function $u(s,\tau)$ satisfy the Black-Scholes equation with a
volatility coefficient $\sigma$. We can rewrite the equation as:

$\displaystyle\frac{\partial{u(x,\tau)}}{\partial{\tau}}+\frac{1}{2}\sigma^{2}x^{2}\frac{\partial^{2}{u(x,\tau)}}{\partial{x^{2}}}=0$
(2.3)

with the initial condition

$\displaystyle u(x,0)$ $\displaystyle=f(x)$ (2.4)

Given the payoff function $f(x)$, which x is the stock price, and at maturity
time, we have $\tau=0$. With the assumption that gives risk-free interest rate
to be zero, the function $u(x,\tau)$ is the price of an option at time t. [9]

###  2.2 Quasi-Reversibility method

We introduce the method to solve Black-Scholes equations with some
approximations. We want to solve the below equation

$\displaystyle Pu$ $\displaystyle=u_{t}+\frac{1}{2}\sigma^{2}(t)x^{2}u_{xx}=0$
(2.5)

with the Dirichlet boundary conditions

$\displaystyle u(x_{b},t)=o_{b}(t),u(x_{a},t)=o_{a}(t)$ (2.6)

and initial conditions

$\displaystyle u(x,0)=f(x)$ (2.7)

where

$\displaystyle f_{x}$ $\displaystyle=\frac{x_{a}}{x_{b}}-1$ (2.8)
$\displaystyle f_{o}$ $\displaystyle=\frac{o_{a}}{o_{b}}-1$ (2.9)

where $t\in[0,2\tau],x\in(x_{b}(0),x_{a}(0))$
P represents the Black-Scholes equation’s partial differential operator. $x$
represents the stock price. $\sigma$ represent the volatility of the option
from the data. $o_{b},o_{a}$ represent the bid and ask price of the option and
$x_{b},x_{a}$ represent the bid and ask price of the stock. Also, we have the
assumption that the bid price is always less than the ask price.

Given the above conditions, the method will take 5 steps to solve for the
Black-Scholes equations.

  * •

First, we will need to solve for the variable that is dimensionless, namely
the $x$’s.

  * •

Secondly, we want to perform some interpolation and extrapolation to forecast
the option price.

  * •

Thirdly, we want to apply the initial boundary condition to the PDE.

  * •

Fourthly, we would like to apply regularization to try to solve the method
numerically.

  * •

Lastly, solve the minimizing problem with finite difference and
regularization. [9]

###  2.3 Analysis

Assume we have $Y_{2\tau}=(x_{b}(0),x_{a}(0))\times(0,2\tau)$the function

$\displaystyle o_{b}(t),o_{a}(t)\in H^{2}[0,2\tau],\sigma(t)\in
C^{1}[0,2\tau]$ (2.10)

which we can find the solution $u\in H^{2}(Y_{2\tau})$ which satisfy the
following conditions

$\displaystyle Au=0\in Y_{2\tau}$ (2.11) $\displaystyle
u(0,t)=o_{b}(t),u(1,t)=o_{a}(t)$ (2.12) $\displaystyle u(x,0)=f(x),x\in(0,1)$
(2.13)

We are going to perform a more general convergence analysis base on a new
Carleman estimate, where we assume the time interval is within a short period.
Let’s formulate a rather general case with the given property.
Given an non-negative constant L and let $Z_{L}$ be
${(x,t)\in(0,1)\times(0,L)}$. Then, let’s have two non-negative number
$\gamma_{1},\gamma_{2}$ which $\gamma_{2}>\gamma_{1}$. Suppose there are
function $\Gamma(x,t)\in C^{1}(\bar{Z}_{L})$ that satisfy

$\displaystyle||b||_{C^{1}(\bar{Z}_{L})}\leq\gamma_{2},\Gamma(x,t)\geq\gamma_{1}\in
Z_{L}.$ (2.14)

Then, the rather general version of the original problem becomes

$\displaystyle Mk=k_{t}+\Gamma(x,t)k_{x}x=0\in Z_{L}$ (2.15)

With the boundary condition

$\displaystyle k(0,t)=o_{b}(t),k(1,t)=o_{a}(t)$ (2.16)

and the initial condition

$\displaystyle k(x,0)=k(x)=o_{b}(0)(1-x)+o_{a}(0)x,x\in(0,1)$ (2.17)

We are trying to find a convergence solution to the general solutions. [11]

##  3 Classification algorithms

###  3.1 Decision Tree:

A decision tree is a classification algorithm based on a tree model. Such a
tree consists of a root node with no incoming edge, a decision node with
precisely one incoming edge and one outgoing edge, and leaf nodes with only an
incoming edge with no outgoing edge. The incoming edge is the result of the
previous decision, and the outgoing edge is the result of current decisions.
Also, it is possible to have a branch that consists of multiple nodes. Note
that all the data will classify the data into a final class. The idea of a
decision tree is elementary. By having some features of your data, we will
split the data based on what features are in the node, where each edge coming
out from the node represent the possible answer from the decision we made. [1]
The decision tree does not have to be a binary tree where you can have more
than two decisions made out of one node. In our case, we will only use a
binary tree because we only wanted to capture the option movement as increase
or decrease. We will feed the tree with minimizer and market features and
split the data into training, validation, and test sets. Then, we will create
trees with different combinations of features of our data to generate outcomes
for trading strategy. Below is a graph showing the basic binary decision,
which we call a classification.1.

Figure 1: Binary Decision Tree

###  3.2 Data Split:

Data split is the procedure to determine the decisions on each node.
Generally, we will have the feature which splits the data roughly in half as
our root note which makes the tree balanced while making the leaf node have
bias toward one of the classifying classes so we can have a final decision.

####  3.2.1 Entropy:

Entropy is the measure of the disorder of the feature relative to our data. In
another word, it determines how well we can classify the data into a certain
class. The formula is given below:

$\displaystyle H(S)=\sum_{i=1}^{c}-p_{i}\log_{2}{p_{i}},$ (3.1)

where $H(S)$ is the Entropy, or disorder, of the feature, $p_{i}$ is the
probability of feature $i$, and c is the total feature that we have. Entropy
can range from $[0,1]$, with the extremes, 0 representing the lowest Entropy,
signaling no disorder, which is bad for training because everything got
classified into just one class. On the other hand, an entropy of 1 signals a
good learning set because it splits up the data evenly. Furthermore, we can
use Entropy to determine the information gain of each node which is how much
information we can learn from a feature.

####  3.2.2 Information Gain:

Information Gain is the measure of change in Entropy due to splitting the data
on a feature. The formula is as follows:

$\displaystyle\mbox{Information
Gain}(S,X)=H(S)-\sum_{i=1}\frac{x_{i}}{S}H(x_{i}),$ (3.2)

where $H(S)$ is the measure of the disorder before the split, $x_{i}$ is the
number of observations on the $ith$ split, and $H(x_{i})$ is the disorder for
the target feature on split $x_{i}$. By subtracting the combined Entropy of
all the child nodes from the Entropy of the parent node, we can observe how
much ”information” we have gained or lost. The greater the information gain,
the more Entropy, or disorder, decreases. Information Gain is utilized to tell
whether or not a split was more or less beneficial based on how much ”noise”
was removed due to the split.

####  3.2.3 Gini Impurity:

Gini Impurity is the measure of the probability that an instance of a random
variable is incorrectly classified if that new instance were randomly
classified according to the distribution of class labels from the data set.
The formula is presented below:

$\displaystyle Gini=1-\sum_{i=1}^{n}(p_{i})^{2},$ (3.3)

n where n is the number of classes and $p_{i}$ is the probability of an item
with label $i$. Thus, to calculate the Gini Impurity, one summates the
probability, $p_{i}$, times the probability of a mistake in classifying that
item, $\sum_{k\neq i}p_{k}$. This leads Gini Impurity to vary between values
of 0 and 1, where a 0 represents complete purity of classification, a 0.5
represents an equal distribution of elements over some classes, and a 1
represents the random distribution of elements over multiple classes.

###  3.3 Hyper Parameter Selection:

Hyper-parameters define a model much like how the rules of a sport define the
sport. One can differentiate between models through the difference in Hyper-
Parameters as Hyper-Parameters address questions regarding the model design.
These questions might include the maximum depth of a tree, the maximum amount
of nodes, and the learning rate of gradient descent. These Hyper-Parameters
are not features of the model, and instead make up the structure of the model.
These metrics determine how the algorithm behaves when it is ran, and so
tuning for these Hyper-parameters proves essential to minimize noise and
provide higher accuracy. The selection process can vary, as you can either
change thee Hyper-parameters yourself and hope to see better results, or one
can use a method such as Grid Search to tune for the most efficient
combinations of the Hyper-parameters. The Grid Search is essentially building
a model for each and every possible assortment of Hyper-parameters, grading
each model, and selecting the one with the most promising results. Another
selection method would be to generate random Hyper-Parameter values from a
data distribution for each parameter. This process is called Random Search,
due to the randomness created when not using a definite set of values, but
rather random sampling the values from a distribution. [6]

###  3.4 Ensemble methods:

The ensemble method uses groups of based learners or models to form a final
decision by taking the average among different learners or models. It is more
or less like voting, where the decision is based on the majority opinion.
Since individual models might create high bias or variance, using the model
combination will likely offer a lower bias and variance by taking the random
samples of the same model or proposing a different model to train the data. It
is the key to the bagging and boosting method, where bagging trains the weaker
learner in parallel, whereas boosting trains them sequentially. We would like
to use such methods because we would want to eliminate the possible high bias.

###  3.5 Gradient Descent:

Gradient descent is a first-order iterative optimization method that tries to
navigate the local minima of some algorithm. It is a boosting method that is
trying to reduce the inaccuracy of our tree where it turns a weak model into a
stronger model by fixing the weakness. The idea is to look at the gradient of
the current location and then take an opposite step to the gradient, which
results in the steepest step towards the minima. The method is analogous to
going down the hill with the steepest path. Also, we have something called the
learning rate, which is the step size we are taking each time we are trying to
go down the hill. A larger learning rate indicates a bigger step that might go
over the minima and get to uphills. While taking a small step will grant you
more accuracy but sufferers from inefficiency and saddle points which tricks
the machine into that we have reached a minima. In the method, we will use
something called a loss function to determine our movement on the hill. It is
essential the difference between the actual output of our tree compares to our
prediction, which we call the error. We will try to minimize the loss function
until it is about zero, which is hopefully the global minimum. There are also
challenges in using the algorithm. It is a greedy method that only considers
the current location, and it can lead to a local minima, which is not
necessarily the best solution. [4]

Figure 2: Binary Decision Tree

####  3.5.1 Gradient Classification:

Here, we would like to apply gradient descent to our decision tree to boost
our accuracy. First, we would like to introduce the mean squared error, which
is used to evaluate the quality of our splits. Since we are solving an
optimization problem, in particular, trying to have a classification to
determine the label of different classes, we would use the log-loss function
and minimize it for better accuracy. Needs to compute the log-loss derivative
here

The idea is we are going to have some input data and a differentiable loss
function. In the classification case, we are going to use a log-loss function,
which has the form:

$\displaystyle\mbox{total loss}=-\frac{1}{N}\sum_{i=1}^{N}y_{i}\cdot
log(p(y_{i}))+(1-y_{i})\cdot(1-p(yi))$ (3.4)

The reason is that the logarithmic loss is a metric that tests the prediction
probability corresponding to the true label. With the loss function, we can
now try to do the iterations.

Step 1

$\displaystyle F_{0}(x)=argmin\sum_{i=1}^{n}L(y_{i},\gamma)$ (3.5)

  1. 1.

$y_{i}$ \- The output(observed value)

  2. 2.

$L$ \- The loss function

  3. 3.

$\gamma$ \- value for log(odds)

Now we want to add up the loss function for each observed value. We want to
find the gamma value that minimizes the sum. We will simply take the
derivative of each loss function.

Step 2: for m = 1 to M

$\displaystyle r_{im}=-[\frac{\partial L(y_{i},F(x_{i})}{\partial
F(x_{i})}]_{F(x)=F(m-1)(x)}\forall i=1,2,...,n$ (3.6)

We want to calculate the pseudo residuals because it is the intermediate error
term between the actual value and the intermediate predicted value. We have

$\displaystyle Observed-Predicted=PseudoResidual$ (3.7)

Then, we want to fit a regression tree that is closed under scaling $h_{m}(x)$
to pseudo residuals. After that, we want to compute multiplier $\gamma_{m}$ by
solving the one-dimensional optimization problem

$\displaystyle\gamma_{m}=argmin\sum_{i=1}^{n}L(y_{i},F_{m-1}(x_{i})+\gamma
h_{m}(x_{i}))$ (3.8)

Finally, wee will update our model to be

$\displaystyle F_{m}(x)=F_{m-1}(x)+\gamma_{m}h_{m}(x)$ (3.9)

We will simply output the $F_{m}(x)$ to be our new decision tree model. [6]
By differentiating the log-loss function at very iterations, we are able to
improve our decision tree accuracy by a large margin. Also, the use of cross-
validation eliminates some of the over-fitting when we get deeper trees that
exhibit higher train accuracy.

###  3.6 Random Forest

We are using a Bagging method which selects a subset of the data and uses it
to create many different training sets. Then, we will train these sets
independently and average the majority of the prediction that has been made to
acquire a more accurate estimation while eliminating over-fitting during the
process. [5] In this project, we are intended to use a random forest
algorithm, which utilizes the bagging method to create many random
uncorrelated decision trees, which is a forest. When generating such decision
trees, we purposely generate trees with low correlation, known as ”the random
subspace method.” While the decision tree utilized all the provided features,
random forest only considered a fraction of the features. The data split is
done by either Gini, entropy, or log-loss criterion. We will decide how we
want to split the tree based on hyper-parameter turnings. On the other hand,
to extract an important feature, we will calculate how important a certain
feature is by calculating when the probability of reaching the node decreases
and how much-weighted node impurity will decrease.

$\displaystyle imp(i)=w_{i}N_{i}-w_{left(i)}\cdot
N_{left(i)}-w_{right(i)}\cdot N_{right(i)}$ (3.10)

  * •

$w_{i}$ = weighted samples that reaches node i

  * •

$N_{i}$ = impurity of node i

  * •

left = left split on node i

  * •

right = right split on node i [6]

Due to its correctness, the random forest can be an extremely time-consuming
process because it trains many trees simultaneously. On the other hand, it is
also hard to interpret the result of the decision tree due to the gigantic
forest that was being created. That is why we use feature importance to help
us interpret the result. Finally, We would like to apply random forest to our
data because we believe that the market follows certain patterns. If we can
make the right decision while the data exhibit certain feature, we are able to
extract from it and improve the prediction we made from QRM. Furthermore,
since our only goal is to classify future options as either increase or
decrease, having majority votes with different trees can reduce over-fitting
and bias within the short period of option price predictions.

##  4 Result and Data

###  4.1 Classification matrix

We define the option price of a certain day to be the average of the bid and
the asking price. All the result below was calculated based on the following
matrix

  * •

Predicted label of 1 represents that we predicted the option price would
increase the next day.

  * •

Predicted label of -1 represents that we predicted the option price would
decrease the next day.

  * •

True label of 1 represents that the real option price increase the next day.

  * •

True label of -1 represents that the real option price decrease the next day.

###  4.2 Decision Tree

Method | Accuracy | Precision | Recall
---|---|---|---
Decision Tree | 53.90% | 58.24% | 62.66%

Figure 3: Confusion Matrix Figure 4: Decision Tree

###  4.3 Decision Tree with Gradient Boosting

Method | Accuracy | Precision | Recall
---|---|---|---
Gradient Boosting | 65.74% | 74.26% | 56.50%

Figure 5: Confusion Matrix Figure 6: Gradient Boosting Tree

###  4.4 Random Forest

Method | Accuracy | Precision | Recall
---|---|---|---
Random Forest | 69.02% | 70.96% | 75.68%

Figure 7: Confusion Matrix Figure 8: Random Forest Tree

Since the decision tree is too deep for the purpose of demonstration. We will
only show a partial tree which represents the first few layers of the tree.

###  4.5 Stock market Data:

Figure 9: Empirical Data

The data that has been used for this project are from bloomberg.com[3].

###  4.6 Result

Method | Accuracy | Precision | Recall
---|---|---|---
Decision Tree | 53.90% | 58.24% | 62.66%
Gradient Boosting | 65.74% | 74.26% | 56.50%
Random Forest | 69.02% | 70.96% | 75.68%

##  5 Summary and Conclusions

Although we cannot correctly predict the exact option price due to the market
volatility, we can make a good prediction of whether the option price will
increase or decrease based on the market features and the minimizer we
acquired from QRM. Based on our assumptions, we will consider the percentage
of precision as our profitable buy for the option. Notice that we had 58.24%
recall when only the Decision Tree was being applied. Improving our result
with gradient boosting, we are able to get 74.26%. Last but not least, we can
achieve 70.96% precision with random forest. The data suggest that we are able
to make a profit out of the market by applying the above methods.

However, since the current metric only examines if a stock option will
increase or decrease in value one day in advance, it does not consider the
fees and the duration of the option. Furthermore, the data gathered is within
a short time frame, which does not reflect the market’s volatility. In the
future, we will implement methods that predict the percentage of increase or
decrease relative to the current option price, so there is a high margin and
chance for people to make money in options trading. We can improve the model
accuracy by applying different machine learning models than tree
classifications while supplying a finer minimizer from the QRM.

##  6 Appendix

Derivation of the log_loss function

First We have some Input for our classification:

$Data(x_{i},y_{i})_{i=1}^{n}$

Also, We have some lost function

$L(y_{i},F(x))$

  1. 1.

$x_{i}$ \- The input data that we feed into the model

  2. 2.

$y_{i}$ \- The output we are trying to predict

Likelihood Function: The likelihood function describes the joint probability
of the observed data as a function of parameters of our input.

$L(\theta)=L(\theta;y)=f_{y}(y;\theta)$

Likelihood function is going to give us the product of all such probability

$L(\theta)=\prod_{i=1}^{n}f_{i}(y_{i};\theta)$

Here we are using the log likelihood function which is

$log\\_likelihood=[y_{i}\cdot log(p)+(1-y_{i})\cdot log(1-p)]$

  1. 1.

we require $y_{i}$ to be binary values which in this case 0 or 1.

  2. 2.

p represent the predicted probability.

Since we want to have a better prediction which means that we would want to
maximize the log likelihood function. If we want to use this as our loss
function, we will simply invert the function which

$log\\_loss=-log\\_likelihood\cdot$

If we minimize the loss function, we essentially will have a better fitting.
Since the likelihood function is a function of probability, we want to make it
into a function of corresponding odds is the logit function which has the form

$logit(p)=log(\frac{p}{1-p})$

which $\frac{p}{1-p}$ represents the corresponding odds of the given
probability. Let’s do some basic arithmetic by plugging in variables

$\displaystyle log\\_loss$ $\displaystyle=-[y_{i}\cdot log(p)+(1-y_{i})\cdot
log(1-p)]$ (6.1) $\displaystyle=-[y_{i}\cdot log(p)]-log(1-p)+y_{i}\cdot
log(1-p)$ (6.2) $\displaystyle=-y_{i}\cdot(log(p)-log(1-p))-log(1-p)$ (6.3)
$\displaystyle=-y_{i}\cdot log(\frac{p}{1-p})-log(1-p)$ (6.4)

First, let’s have $log(\frac{p}{1-p})=y$ Then,

$\displaystyle y$ $\displaystyle=log(\frac{p}{1-p})$ (6.6) $\displaystyle
e^{y}$ $\displaystyle=\frac{p}{1-p}$ (6.7) $\displaystyle e^{y}$
$\displaystyle=-1+\frac{1}{1-p}$ (6.8) $\displaystyle e^{y}+1$
$\displaystyle=\frac{1}{1-p}$ (6.9) $\displaystyle 1-p$
$\displaystyle=\frac{1}{e^{y}+1}$ (6.10) $\displaystyle-p$
$\displaystyle=\frac{1-e^{y}-1}{e^{y}+1}$ (6.11) $\displaystyle p$
$\displaystyle=\frac{e^{y}}{1+e^{y}}$ (6.12) $\displaystyle p$
$\displaystyle=\frac{e^{log(\frac{p}{1-p})}}{1+e^{log(\frac{p}{1-p})}}$ (6.13)

Then, we have

$\displaystyle log(1-p)$
$\displaystyle=log(1-\frac{e^{log(\frac{p}{1-p})}}{1+e^{log(\frac{p}{1-p})}})$
(6.14)
$\displaystyle=log(\frac{1+e^{log(\frac{p}{1-p})-1}}{1+e^{log(\frac{p}{1-p})})}$
(6.15) $\displaystyle=log(\frac{1}{1+e^{log(\frac{p}{1-p})})}$ (6.16)
$\displaystyle=log(1)-log(1+e^{log(\frac{p}{1-p})})$ (6.17)
$\displaystyle=-log(1+e^{log(\frac{p}{1-p})})$ (6.18)

Let’s put the equation back into log_loss function.

$\displaystyle log\\_loss$ $\displaystyle=-y_{i}\cdot
log(\frac{p}{1-p})-log(1-p)$ (6.19) $\displaystyle=-y_{i}\cdot
log(\frac{p}{1-p})+log(1+e^{log(\frac{p}{1-p})})$ (6.20)

Let’s differential the loss function which

$\displaystyle\frac{d}{dlog(\frac{p}{1-p})}(-y_{i}\cdot
log(\frac{p}{1-p})+log(1+e^{log(\frac{p}{1-p})}))=-y_{i}+\frac{e^{log(\frac{p}{1-p})}}{1+e^{log(\frac{p}{1-p})}}$
(6.21)

Which we can rewrite it as, Predicted - Observed where

$\displaystyle Observed=y_{i}$ (6.22) $\displaystyle
predicted=\frac{e^{log(\frac{p}{1-p})}}{1+e^{log(\frac{p}{1-p})}}$ (6.23)

[10]

## References

  * [1] L. Breiman, J. Friedman, R. Olshen, and C. Stone, _Classification and Regression Trees_ Wadsworth, Belmont, CA, 1984.
  * [2] Michael V Klibanov and Andrey V Kuzhuget and Kirill V Golubnichiy, _An ill-posed problem for the Black–Scholes equation for a profitable forecast of prices of stock options on real market data_ , Inverse Problems, 32(1), 2016.
  * [3] https://bloomberg.com.
  * [4] Lemaréchal, C, _Cauchy and the Gradient Method_ Doc Math Extra: 251–254, (2012)
  * [5] Breiman, L, _Random Forests_ Machine Learning 45, 5–32, (2001).
  * [6] Pedregosa, F. and Varoquaux, G. and Gramfort, A. and Michel, V. and Thirion, B. and Grisel, O. and Blondel, M. and Prettenhofer, P. and Weiss, R. and Dubourg, V. and Vanderplas, J. and Passos, A. and Cournapeau, D. and Brucher, M. and Perrot, M. and Duchesnay, E., _Scikit-learn: Machine Learning in Python_ , Journal of Machine Learning Research, 12(2825–2930), 2011.
  * [7] Klibanov, Michael V. and Shananin, Aleksander A. and Golubnichiy, Kirill V. and Kravchenko, Sergey M., _Forecasting Stock Options Prices via the Solution of an Ill-Posed Problem for the Black-Scholes Equation_ , arXiv, 2022
  * [8] Pedro Santa-Clara, Alessio Saretto, _Option strategies: Good deals and margin calls_ , Journal of Financial Markets, Volume 12, Issue 3, 2009
  * [9] Mikhail V. Klibanov, Kirill V. Golubnichiy and Andrey V. Nikitin, _Quasi-Revesibility Method and Neural Network Machine Learning for Forecating of Stock Option Prices_ , Inverse Problem, 2022
  * [10] Jurafsky, Daniel and Martin, James H., _Speech and Language Processing: An Introduction to Natural Language Processing, Computational Linguistics, and Speech Recognition_ , Prentice Hall PTR (3rd ed.). pp 76-95, 2020
  * [11] M.V. Klibanov, A. A. Shananin, K. V. Golubnichiy and S. M. Kravchenko, _Forecasting Stock Options Prices via the Solution of an Ill-Posed Problem for the Black-Scholes Equation_ , arXiv preprint arXiv:2202.07174.