# Non-Bloch bands in two-dimensional non-Hermitian systems

Kazuki Yokomizo  Department of Physics, The University of Tokyo, 7-3-1 Hongo,
Bunkyo-ku, Tokyo, 113-0033, Japan  Shuichi Murakami  Department of Physics,
Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo, 152-8551,
Japan

###### Abstract

The non-Bloch band theory can describe energy bands in a one-dimensional (1D)
non-Hermitian system. On the other hand, whether the non-Bloch band theory can
be extended to higher-dimensional non-Hermitian systems is nontrivial. In this
work, we construct the non-Bloch band theory in two classes of two-dimensional
non-Hermitian systems, by reducing the problem to that for a 1D non-Hermitian
model. In these classes of systems, we get the generalized Brillouin zone for
complex wave vectors and investigate topological properties. In the model of
the non-Hermitian Chern insulator, as an example, we show the bulk-edge
correspondence between the Chern number defined from the generalized Brillouin
zone and the appearance of the edge states.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

A non-Hermitian Hamiltonian effectively describes a nonequilibrium system [1].
In recent studies, the non-Hermitian skin effect, i.e. the localization of
bulk eigenstates, plays a crucial role because it causes remarkable phenomena
[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. In particular,
with the non-Hermitian skin effect, energy eigenvalues under an open boundary
condition and those under a periodic boundary condition are different. Such
phenomena associated with the non-Hermitian skin effect are an obstacle to
investigate properties of a non-Hermitian system. For example, the bulk-edge
correspondence between a topological invariant defined from the conventional
Bloch wave vectors and existence of topological edge states seems to be
violated within the conventional definition of a topological invariant.

Recently, the non-Bloch band theory was proposed [2, 19, 20, 21, 22]. It was
shown that in a one-dimensional (1D) non-Hermitian system, the generalized
Brillouin zone $\beta=e^{ik}$ for the complex Bloch wave number $k$ becomes
different from that in a Hermitian system, and it reproduces energy
eigenvalues under an open boundary condition. Importantly, a topological
invariant defined from the generalized Brillouin zone restores the bulk-edge
correspondence. Additionally, some features of the generalized Brillouin zone
lead to unique phenomena, such as appearance of a topological semimetal phase
with exceptional points [23]. Thus, the non-Bloch band theory is a powerful
tool for studies on many properties in a non-Hermitian system. However, there
is a lack of studies on a two-dimensional (2D) non-Hermitian system in terms
of the non-Bloch band theory.

There have been some studies on topological phases and boundary states in a 2D
non-Hermitian system in terms of the conventional Bloch band theory. For
example, some previous works investigated how non-Hermiticity affects the
bulk-edge correspondence [24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35].
Furthermore, Refs. [36, 37, 38, 39] calculated the Hall conductance associated
with the topological edge states in the model of the non-Hermitian Chern
insulator. In recent years, non-Hermitian topology in a higher dimensional
system has been attracting much attention. In theory, Refs. [40, 41, 42, 43,
44, 45, 46, 47] proposed the higher-order topological non-Hermitian skin
effect in which eigenstates are localized at corners of the system, and in
experiment, it has been observed [48, 49, 50]. Thus, many topological aspects
of a 2D non-Hermitian system have been studied in terms of the real Bloch wave
vectors. However, as we learned from a 1D non-Hermitian system, we cannot
investigate physics of a 2D non-Hermitian system with an open boundary
condition in terms of the real Bloch wave vector. Therefore, it is necessary
to construct the non-Bloch band theory in a 2D non-Hermitian system.

So far, the non-Bloch bands in some 2D non-Hermitian models were investigated
[51, 52, 53, 54, 55, 56]. Nevertheless, it is unclear whether the non-Bloch
band theory in general 2D non-Hermitian systems can be constructed or not. In
this work, we construct the non-Bloch band theory in two classes of 2D non-
Hermitian systems, where we can reduce the problem to that of a 1D non-
Hermitian model. One class of systems has specific symmetries, such as a
mirror symmetry suppressing the non-Hermitian skin effect in one direction.
The other class is decoupled into two 1D non-Hermitian systems due to a
special form of the Hamiltonian. These properties indicate that the
eigenstates in the bulk are written as linear combination of a few plane
waves. Thereby, the condition for the generalized Brillouin zone can be
obtained. In this paper, we exemplify three systems, e.g., the model of the
non-Hermitian Chern insulator, the non-Hermitian Benalcazar–Bernevig–Hughes
(BBH) model, and the Okugawa-Takahashi-Yokomizo (OTY) model. We show that the
generalized Brillouin zone reproduces the energy eigenvalues in these models
with a full open boundary condition. Furthermore, we investigate topological
properties in terms of the generalized Brillouin zone.

This paper is organized as follows. In Sec. II, we introduce our 2D non-
Hermitian tight-binding model and propose two classes of systems where we can
apply the non-Bloch band theory. Then, we can get the condition for the
generalized Brillouin zone. In Sec. III, we analyse the three 2D non-Hermitian
models in terms of the non-Bloch band theory. In Sec. IV, we discuss the 2D
non-Bloch band theory proposed in this work from the viewpoint of the
formation of standing wave and comment on difficulty to construct the non-
Bloch band theory in general 2D non-Hermitian systems. Finally, in Sec. V, we
summarize our result.

##  II Model

Figure 1: Two-dimensional tight-binding system on a square lattice. The
geometry of the system is a rectangle, where the lengths in the $x$ direction
and in the $y$ direction are $L_{x}$ and $L_{y}$, respectively. The black
circles express the unit cells at the lattice point
${\bm{r}}=\left(n_{x},n_{y}\right)~{}\left(n_{x}=1,\dots,L_{x},n_{y}=1,\dots,L_{y}\right)$.

First of all, we introduce a 2D non-Hermitian tight-binding system. Throughout
this paper, the system lies on a square lattice, and the geometry of the
system is a rectangle, as shown in Fig. 1. The real-space Hamiltonian of our
system is given by

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{{\bm{r}}}\sum_{\mu,\nu=1}^{q}\left(\sum_{i=-N_{x}}^{N_{x}}t_{i,\mu\nu}^{\left(x\right)}c_{{\bm{r}}+i\hat{\bm{x}},\mu}^{\dagger}c_{{\bm{r}},\nu}\right.$
(1)
$\displaystyle\left.+\sum_{i=-N_{y}}^{N_{y}}t_{i,\mu\nu}^{\left(y\right)}c_{{\bm{r}}+i\hat{\bm{y}},\mu}^{\dagger}c_{{\bm{r}},\nu}\right),$

where $c_{{\bm{r}},\mu}^{\dagger}$ is a creation operator,
${\bm{r}}=\left(n_{x},n_{y}\right)~{}\left(n_{x}=1,\dots,L_{x},n_{y}=1,\dots,L_{y}\right)$
is a lattice point, $q$ is the number of internal degrees of freedom in a unit
cell, and $\hat{\bm{x}}$ and $\hat{\bm{y}}$ are unit vectors in the $x$
direction and in the $y$ direction, respectively. The particles hop up to the
$N_{x}$th and $N_{y}$th nearest-neighbor unit cells along the $x$ and $y$
directions, respectively. We note that when either
$t_{i,\nu\mu}^{\left(x\right)}\neq\left(t_{-i,\mu\nu}^{\left(x\right)}\right)^{\ast}$
or
$t_{i,\nu\mu}^{\left(y\right)}\neq\left(t_{-i,\mu\nu}^{\left(y\right)}\right)^{\ast}$
is satisfied, the system becomes non-Hermitian. Now, the real-space eigen-
equation is written as

$H\left|{\psi}\right\rangle=E\left|\psi\right\rangle,$ (2)

where the eigenstates are given by

$\left|\psi\right\rangle=\left(\dots,\psi_{\left(n_{x},n_{y}\right),1},\dots,\psi_{\left(n_{x},n_{y},1\right),q},\dots\right)^{\rm
T}.$ (3)

Next, we describe a way to construct the non-Bloch band theory in our system.
In the case of $q=1$, it is straightforward to construct the non-Bloch band
theory in the system. This is because the system can be decomposed into two 1D
non-Hermitian models. In Appendix A, we derive the condition for the
generalized Brillouin zone in this case. On the other hand, in the case of
$q\geq 2$, we find that when the system satisfies some conditions, one can
calculate the generalized Brillouin zone in terms of the conventional non-
Bloch band theory. In the following, we explain the two cases where the
description of the non-Bloch bands is possible.

###  II.1 Case A

In a 1D non-Hermitian system, when the energy of the eigenstate with the Bloch
wave number $+k$ and that of the eigenstate with $-k$ are degenerate, the non-
Hermitian skin effect is suppressed because $k$ becomes real [57, 58, 59]. The
degeneracy leading to the suppression of the non-Hermitian skin effect comes
from symmetries, e.g., a global symmetry, such as a PT symmetry, and a
crystalline symmetry, such as a mirror symmetry. Importantly, we here propose
that such suppression is realized also in a 2D non-Hermitian system. For
example, when the Bloch Hamiltonian $H\left(k_{x},k_{y}\right)$ satisfies

$U^{-1}H^{\rm T}\left(k_{x},-k_{y}\right)U=H\left(k_{x},k_{y}\right),$ (4)

where $U$ is a unitary matrix satisfying $U^{2}=1$, we can show that the non-
Hermitian skin effect in the $y$ direction is suppressed. In fact, when Eq.
(4) is satisfied, the wave with the Bloch wave vector
$\left(k_{x},k_{y}\right)$ is reflected to that with
$\left(k_{x},-k_{y}\right)$ at the boundary of the system parallel to the $x$
direction. Then, the condition for the formation of standing wave means that
the imaginary parts of wave vectors for the incident and reflected waves are
equal. Therefore, we can obtain ${\rm Im}\left(k_{y}\right)={\rm
Im}\left(-k_{y}\right)$, meaning that $k_{y}$ is real, and the non-Hermitian
skin effect in the $y$ direction is suppressed.

In this case, the bulk eigenstates extend over the $y$ direction because the
Bloch wave number $k_{y}$ becomes real. Hence, in the limit of a large system
size, the asymptotic behavior of the bulk eigenstates in the cylinder
geometry, where a periodic boundary condition in the $y$ direction is imposed
to the system, matches that in the full-open geometry. Therefore, the energy
bands in the cylinder geometry asymptotically reproduces the energy levels in
a finite open plane. Importantly, since the cylinder geometry can be regarded
as a pseudo-1D system, one can calculate the energy bands of this system in
terms of the conventional non-Bloch band theory. Thus, we can construct the
non-Bloch band theory in the system.

Based on the above concept, we establish the non-Bloch band theory in the 2D
non-Hermitian systems without the non-Hermitian skin effect in the $y$
direction. With Eq. (4), by assuming that the Bloch wave number $k_{y}$ is
real, we take

$\psi_{{\bm{r}},\mu}=\sum_{j}\left(\beta_{x,j}\right)^{n_{x}}e^{ik_{y}n_{y}}\phi_{\mu}^{\left(j\right)}$
(5)

as an ansatz for the eigen-equation in the cylinder geometry. Here, we define
the non-Bloch matrix as

$\left[{\cal
H}\left(\beta_{x},k_{y}\right)\right]_{\mu\nu}=\sum_{m=-N_{x}}^{N_{x}}t_{m,\mu\nu}^{\left(x\right)}\left(\beta_{x}\right)^{m}+\sum_{m=-N_{y}}^{N_{y}}t_{m,\mu\nu}^{\left(y\right)}e^{ik_{y}m}.$
(6)

Then, in Eq. (5), $\beta_{x}=\beta_{x,j}$ is the solution of the
characteristic equation

$\det\left[{\cal H}\left(\beta_{x},k_{y}\right)-E\right]=0.$ (7)

We note that Eq. (7) is an algebraic equation for $\beta_{x}$ of $2qN_{x}$th
degree. Therefore, the condition for the generalized Brillouin zone is given
by

$\left|\beta_{x,qN_{x}}\right|=\left|\beta_{x,qN_{x}+1}\right|,$ (8)

where the solutions of Eq. (7) are numbered as

$\left|\beta_{x,1}\right|\leq\dots\leq\left|\beta_{x,2qN_{x}}\right|.$ (9)

For a given real value of $k_{y}$, we can obtain trajectories of
$\beta_{x,qN_{x}}$ and $\beta_{x,qN_{x}+1}$ with Eq. (8) forming loops on the
complex plane. Thus, by changing the real value of $k_{y}$, the surface formed
by $\left(\beta_{x},\beta_{y}\right)~{}\left(\beta_{y}\equiv
e^{ik_{y}},k_{y}\in{\mathbb{R}}\right)$ is the generalized Brillouin zone.
Finally, from the generalized Brillouin zone, we can calculate the energy
bands by using Eq. (7). In this work, we will show that it matches the energy
levels in a finite open plane. This means that the non-Hermitian skin effect
is indeed suppressed in the $y$ direction in this case.

###  II.2 Case B

Next, we focus on the case that the characteristic equation of the Bloch
Hamiltonian $H\left(k_{x},k_{y}\right)$ is written in the form of separation
of variables. In this case, the motion of the bulk eigenstates in the $x$
direction is completely decoupled from that in the $y$ direction. Namely, the
2D non-Hermitian systems can be regarded as two 1D non-Hermitian systems. In
the following, we assume that the characteristic equation is given by

$\det\left[H\left(k_{x},k_{y}\right)-E\right]=R_{x}\left(e^{ik_{x}},E\right)+R_{y}\left(e^{ik_{y}}\right)=0.$
(10)

Then, by applying the conventional non-Bloch band theory to the 1D non-
Hermitian systems, we expect that one can get the energy bands in the 2D non-
Hermitian system, and indeed, we will show it in this work.

Now, we derive the condition for the generalized Brillouin zone to calculate
the energy bands. With Eq. (10), we have

$\psi_{{\bm{r}},\mu}=\sum_{j_{x},j_{y}}\left(\beta_{x,j_{x}}\right)^{n_{x}}\left(\beta_{y,j_{y}}\right)^{n_{y}}\phi_{\mu}^{\left(j_{x},j_{y}\right)}$
(11)

as an ansatz for Eq. (2). Here, $\beta_{x}=\beta_{x.j_{x}}$ and
$\beta_{y}=\beta_{y.j_{y}}$ in Eq. (11) are the solutions of the equation
$R_{x}\left(\beta_{x},E\right)=\lambda$ and the equation
$R_{y}\left(\beta_{y}\right)=-\lambda$, respectively, where $\lambda$ is an
arbitrary constant. We note that these equations are algebraic equations of
$2qN_{x}$th degree and $2qN_{y}$th degree, respectively. Thus, we can get the
generalized Brillouin zone spanned by $\left(\beta_{x},\beta_{y}\right)$
satisfying

$\left|\beta_{a,qN_{a}}\right|=\left|\beta_{a,qN_{a}+1}\right|,$ (12)

where the solutions of $R_{x}\left(\beta_{x},E\right)=\lambda$ and
$R_{y}\left(\beta_{y}\right)=-\lambda$ are numbered as

$\left|\beta_{a,1}\right|\leq\dots\leq\left|\beta_{a,2qN_{a}}\right|$ (13)

for $a=x,y$. Finally, the energy bands are obtained from the characteristic
equation of the non-Bloch matrix

$\left[{\cal
H}\left(\beta_{x},\beta_{y}\right)\right]_{\mu\nu}=\sum_{m=-N_{x}}^{N_{x}}t_{m,\mu\nu}^{\left(x\right)}\left(\beta_{x}\right)^{m}+\sum_{m=-N_{y}}^{N_{y}}t_{m,\mu\nu}^{\left(y\right)}\left(\beta_{y}\right)^{m}.$
(14)

We note that unlike Case A, the non-Hermitian skin effect appear in both
directions.

##  III Examples

In this section, we calculate the generalized Brillouin zone and the energy
bands in some examples, e.g., the model of the non-Hermitian Chern insulator,
the non-Hermitian BBH model, and the OTY model. The first model is in Case A,
and the second and third models are in Case B, as discussed in Sec. II.
Throughout this calculation, we show that the non-Bloch band theory in the
Case A and Case B indeed reproduces the energy levels in a finite open plane.
Furthermore, we study topological properties in theses systems in terms of the
generalized Brillouin zone.

###  III.1 Non-Hermitian Chern insulator

First of all, we focus on the model of a non-Hermitian Chern insulator
investigated in Ref. [51]. The real-space Hamiltonian of this system is given
by

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{\bm{r}}\left({\bm{c}}_{{\bm{r}}+\hat{\bm{x}}}^{\dagger}T_{x}^{\dagger}{\bm{c}}_{{\bm{r}}}+{\bm{c}}_{{\bm{r}}}^{\dagger}T_{x}{\bm{c}}_{{\bm{r}}+\hat{\bm{x}}}\right.$
(15)
$\displaystyle\left.+{\bm{c}}_{{\bm{r}}+\hat{\bm{y}}}^{\dagger}T_{y}^{\dagger}{\bm{c}}_{{\bm{r}}}+{\bm{c}}_{{\bm{r}}}^{\dagger}T_{y}{\bm{c}}_{{\bm{r}}+\hat{\bm{y}}}+{\bm{r}}_{{\bm{n}}}^{\dagger}M{\bm{c}}_{{\bm{r}}}\right),$

where

$\displaystyle T_{x}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\begin{array}[]{cc}-t_{x}&-iv_{x}\vspace{3pt}\\\
-iv_{x}&t_{x}\end{array}\right),~{}T_{y}=\frac{1}{2}\left(\begin{array}[]{cc}-t_{y}&-v_{y}\vspace{3pt}\\\
v_{y}&t_{y}\end{array}\right)$ (20) $\displaystyle M$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}m&i\gamma_{x}+\gamma_{y}\vspace{3pt}\\\
i\gamma_{x}-\gamma_{y}&-m\end{array}\right),$ (23)

and ${\bm{c}}_{\bm{r}}^{\dagger}=\left(c_{{\bm{r}},{\rm
A}}^{\dagger},c_{{\bm{r}},{\rm B}}^{\dagger}\right)$ are the creation
operators of particles on sublattices A and B, respectively. For simplicity,
all the parameters are set to be real. We note that this system is in the case
of $N_{x}=1,N_{y}=1$, and $q=2$ in Eq. (1). The real-space eigen-equation is
written as

$T_{x}{\bm{\psi}}_{{\bm{r}}+\hat{\bm{x}}}+T_{x}^{\dagger}{\bm{\psi}}_{{\bm{r}}-\hat{\bm{x}}}+T_{y}{\bm{\psi}}_{{\bm{r}}+\hat{\bm{y}}}+T_{y}^{\dagger}{\bm{\psi}}_{{\bm{r}}-\hat{\bm{y}}}+M{\bm{\psi}}_{\bm{r}}=E{\bm{\psi}}_{\bm{r}},$
(24)

where ${\bm{\psi}}_{\bm{r}}=\left(\psi_{{\bm{r}},{\rm A}},\psi_{{\bm{r}},{\rm
B}}\right)^{\rm T}$. In the following, we investigate the case of
$\gamma_{y}=0$. Then, the Bloch Hamiltonian $H\left(k_{x},k_{y}\right)$ of the
system has the symmetry written as Eq. (4) with $U=1$. Hence, this system is
in Case A as discussed in Sec. II.1, and we expect that the non-Hermitian skin
effect is suppressed along the $y$ direction. By following the discussion in
Sec. II.1, we can get the generalized Brillouin zone and the energy bands.

In the model of the non-Hermitian Chern insulator, the non-Bloch matrix is
given by

${\cal
H}\left(\beta_{x},k_{y}\right)=\beta_{x}T_{x}+\beta_{x}^{-1}T_{x}^{\dagger}+e^{ik_{y}}T_{y}+e^{-ik_{y}}T_{y}^{\dagger}+M.$
(25)

Since the characteristic equation $\det\left[{\cal
H}\left(\beta_{x},k_{y}\right)-E\right]=0$ is a quartic equation for
$\beta_{x}$, the condition for the generalized Brillouin zone is given by

$\left|\beta_{x,2}\right|=\left|\beta_{x,3}\right|,$ (26)

where the solutions of the characteristic equation satisfy Eq. (9). As an
example, we show the generalized Brillouin zone formed from the complex Bloch
wave number $k_{x}$ and the real Bloch wave number $k_{y}$ in Fig. 2(a).
Furthermore, the energy bands are calculated from the generalized Brillouin
zone as shown in Fig. 2(b). This model has two bands, one with ${\rm
Re}\left(E\right)>0$ and the other with ${\rm Re}\left(E\right)<0$. We also
calculate the energy levels of the system on a finite open plane as shown in
Fig. 2(c). Comparing Fig, 2(b) with Fig. 2(c), we confirm that our non-Bloch
bands reproduce the energy levels in a finite open plane.

Next, we investigate a topological aspect of the model of the non-Hermitian
Chern insulator in terms of the Chern number defined from the generalized
Brillouin zone. Let $\left|u_{R,n}\left({\bm{k}}\right)\right\rangle$ and
$\left\langle u_{L,n}\left({\bm{k}}\right)\right|$ denote the right and left
eigenvectors of Eq. (25), respectively, where ${\bm{k}}$ is the complex Bloch
wave vector with $k_{x}$ defined by $e^{ik_{x}}\equiv\beta_{x}$. Here, $n$
expresses a band index, and we set $n=+$ and $n=-$ to represent the energy
bands with ${\rm Re}\left(E\right)>0$ and ${\rm Re}\left(E\right)<0$,
respectively. The right and left eigenvectors satisfy $\left\langle
u_{L,m}\left({\bm{k}}\right)\middle|u_{R,n}\left({\bm{k}}\right)\right\rangle=\delta_{m,n}$.
Then, we define the Chern number as

$\displaystyle C_{n}$ $\displaystyle=$ $\displaystyle\frac{1}{2\pi
i}\int_{T_{\beta}}d{\bm{k}}{\cal B}_{n}\left({\bm{k}}\right),$ (27)
$\displaystyle{\cal B}_{n}\left({\bm{k}}\right)$ $\displaystyle=$
$\displaystyle\left\langle\frac{\partial
u_{L,n}\left({\bm{k}}\right)}{\partial k_{x}}\middle|\frac{\partial
u_{R,n}\left({\bm{k}}\right)}{\partial k_{y}}\right\rangle$ (28)
$\displaystyle-\left\langle\frac{\partial
u_{L,n}\left({\bm{k}}\right)}{\partial k_{y}}\middle|\frac{\partial
u_{R,n}\left({\bm{k}}\right)}{\partial k_{x}}\right\rangle,$

where ${\cal B}_{n}\left({\bm{k}}\right)$ is the Berry curvature. The
definition of Eq. (27) means that the Chern number is defined as an integral
of the Berry curvature over the generalized Brillouin zone denoted by
$T_{\beta}$. Now, we compute the Chern number with $n=-$ by using the Fukui-
Hatsugai-Suzuki method [60]. As shown in Fig. 2(d), the Chern number takes
$-1$ in some region of the parameter $m$. Importantly, the nonzero Chern
number corresponds to appearance of the topological edge states as we discuss
next. Indeed, the topological edge states shown in orange appear in Fig. 2(c).
Thus, the non-Bloch band theory shows the bulk-edge correspondence in the
model of the non-Hermitian Chern insulator.

In general, the Chern number (27) predicts existence of topological edge
states. In the following, we focus on a non-Hermitian system with a line gap
in which energy bands are separated into two sets, set 1 and set 2, by a line
on the complex energy plane. First, we note that the sum of Eq. (27) over all
the bands belonging to the set 1 is quantized, which is shown in Appendix B.
According to Ref. [57], a non-Hermitian system with a line gap can be
transformed to a Hermitian system without closing the gap. The transformation
corresponds to the deformation from the generalized Brillouin zone to the
conventional Brillouin zone. Then, the values of Eq. (27) remain quantized
under the deformation. Hence, after transforming the system, Eq. (27) becomes
the Chern number of the Hermitian system, defined by the conventional
Brillouin zone. In the Hermitian system, the bulk-edge correspondence for the
Chern number is established. Importantly, since the bulk gap remains open
before and after the transformation, the topological properties in the non-
Hermitian system and the Hermitian system are common. Therefore, the nonzero
value of Eq. (27) corresponds to appearance of topological edge states across
the gap between the set 1 and set 2. Thus, we have shown the bulk-edge
correspondence for Eq. (27).

Finally, we comment on the relationship between the present work and Ref.
[51]. This previous work investigated the model of the non-Hermitian Chen
insulator in terms of a low-energy continuum model in which the Bloch
Hamiltonian of Eq. (15) is expanded up to second order of the Bloch wave
vector around $\left(k_{x},k_{y}\right)=\left(0,0\right)$. In the low-energy
continuum model, one can obtain an approximated value of the complex-valued
Bloch wave vector only near $\left(k_{x},k_{y}\right)=\left(0,0\right)$. On
the other hand, in this work, we show a way to get the full generalized
Brillouin zone in the original lattice model. We note that near
$\left(k_{x},k_{y}\right)=\left(0,0\right)$, our generalized Brillouin zone
reproduces the approximated result in the previous work. In particular, the
previous work proposed that the band gap closes at

$m_{c}\simeq
t_{x}+t_{y}+\frac{t_{x}\gamma_{x}^{2}}{2v_{x}^{2}}+\frac{t_{y}\gamma_{y}^{2}}{2v_{y}^{2}},$
(29)

and the Chern number changes its value there. With the values of the system
parameters in our calculation, Eq. (29) gives $m_{c}\simeq 1.91$, and it
matches the value of $m$ at which the Chern number changes in Fig. 2(d).

Figure 2: Generalized Brillouin zone, energy eigenvalue, and Chern number in
the model of the non-Hermitian Chern insulator. (a) Generalized Brillouin zone
spanned by $\left(\beta_{x},\beta_{y}\right)~{}\left(\beta_{y}\equiv
e^{ik_{y}},k_{y}\in{\mathbb{R}}\right)$. (b) Energy band calculated from the
generalized Brillouin zone. (c) Energy level in the system with
$L_{x}=L_{y}=30$. We show the topological edge states in orange. (d) Chern
number for the band with ${\rm Re}\left(E\right)<0$ computed by using the
Fukui-Hatsugai-Suzuki method. The Chern number changes its value at $m=1.91$.
In (a)-(d), the values of the system parameters are set to be
$t_{x}=0.9,t_{y}=1,v_{x}=1.2,v_{y}=0.8,\gamma_{x}=0.2$, and $\gamma_{y}=0$,
and in (a)-(c), $m=1.6$.

###  III.2 Non-Hermitian Benalcazar–Bernevig–Hughes model

Next, we study the non-Hermitian BBH model with the asymmetric intercell
hopping amplitude and the next-nearest neighbor hopping amplitude by following
Ref. [61]. The real-space Hamiltonian of this system is given by

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{\bm{r}}\left({\bm{c}}_{{\bm{r}}+\hat{\bm{x}}}^{\dagger}T_{x}^{\dagger}{\bm{c}}_{{\bm{r}}}+{\bm{c}}_{{\bm{r}}}^{\dagger}T_{x}{\bm{c}}_{{\bm{r}}+\hat{\bm{x}}}\right.$
(30)
$\displaystyle\left.+{\bm{c}}_{{\bm{r}}+\hat{\bm{y}}}^{\dagger}T_{y}^{\dagger}{\bm{c}}_{{\bm{r}}}+{\bm{c}}_{{\bm{r}}}^{\dagger}T_{y}{\bm{c}}_{{\bm{r}}+\hat{\bm{y}}}+{\bm{c}}_{{\bm{r}}}^{\dagger}M{\bm{c}}_{{\bm{r}}}\right),$

where

$\displaystyle T_{x}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}0&0&t_{3}^{x}&0\vspace{3pt}\\\
0&0&0&-t_{3}^{x}\vspace{3pt}\\\ t_{2}^{x}&0&0&0\vspace{3pt}\\\
0&-t_{2}^{x}&0&0\end{array}\right),~{}T_{y}=\left(\begin{array}[]{cccc}0&t_{3}^{y}&0&0\vspace{3pt}\\\
t_{2}^{y}&0&0&0\vspace{3pt}\\\ 0&0&0&t_{3}^{y}\vspace{3pt}\\\
0&0&t_{2}^{y}&0\end{array}\right)$ (39) $\displaystyle M$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}0&t_{1}^{y}+\gamma^{y}/2&t_{1}^{x}+\gamma^{x}/2&0\vspace{3pt}\\\
t_{1}^{y}-\gamma^{y}/2&0&0&-t_{1}^{x}-\gamma^{x}/2\vspace{3pt}\\\
t_{1}^{x}-\gamma^{x}/2&0&0&t_{1}^{y}+\gamma^{y}/2\vspace{3pt}\\\
0&-t_{1}^{x}+\gamma^{x}/2&t_{1}^{y}-\gamma^{y}/2&0\end{array}\right),$ (44)

and ${\bm{c}}_{\bm{r}}^{\dagger}=\left(c_{{\bm{r}},{\rm
A}}^{\dagger},c_{{\bm{r}},{\rm B}}^{\dagger},c_{{\bm{r}},{\rm
C}}^{\dagger},c_{{\bm{r}},{\rm D}}^{\dagger}\right)$ are the creation
operators of particles on sublattices A, B, C, and D, respectively. For
simplicity, all the parameters are set to be real. We note that this system is
in the case of $N_{x}=1,N_{y}=1$, and $q=4$ in Eq. (1). The real-space eigen-
equation of this system is written in the form of Eq. (24) with
${\bm{\psi}}_{\bm{r}}=\left(\psi_{{\bm{r}},{\rm A}},\psi_{{\bm{r}},{\rm
B}},\psi_{{\bm{r}},{\rm C}},\psi_{{\bm{r}},{\rm D}}\right)^{\rm T}$.

In this system, we have the Bloch Hamiltonian as

$H\left({\bm{k}}\right)=H_{x}\left(e^{ik_{x}}\right)\otimes\sigma_{z}+\sigma_{0}\otimes
H_{y}\left(e^{ik_{y}}\right),$ (46)

where

$\displaystyle H_{a}\left(e^{ik_{a}}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}0&R_{+}^{a}\left(e^{ik_{a}}\right)\\\
R_{-}^{a}\left(e^{ik_{a}}\right)&0\end{array}\right),$ (49) $\displaystyle
R_{\pm}^{a}\left(e^{ik_{a}}\right)$ $\displaystyle=$ $\displaystyle
t_{3}^{a}e^{\pm
ik_{a}}+\left(t_{1}^{a}\pm\frac{\gamma^{a}}{2}\right)+t_{2}^{a}e^{\mp ik_{a}}$
(50)

for $a=x,y$, and $\sigma_{0}$ is an identity matrix, and $\sigma_{z}$ is the
$z$ component of the Pauli matrix. We note that $H_{a}$ satisfies
$\sigma_{z}H_{a}\sigma_{z}=-H_{a}$. Then, the characteristic equation of Eq.
(46) is obtained as

$\sum_{a=x,y}R_{+}^{a}\left(e^{ik_{a}}\right)R_{-}^{a}\left(e^{ik_{a}}\right)-E^{2}=0.$
(51)

Thus, we find that the system is in Case B as discussed in Sec. II.2.
Therefore, from the solutions of the equations

$R_{+}^{a}\left(\beta_{a}\right)R_{-}^{a}\left(\beta_{a}\right)=\lambda_{a},$
(52)

the condition for the generalized Brillouin zone is given by

$\left|\beta_{a,2}\right|=\left|\beta_{a,3}\right|,$ (53)

where we have

$\left|\beta_{a,1}\right|\leq\dots\left|\beta_{a,4}\right|$ (54)

for $a=x,y$. Then, we can calculate the energy bands from

$E=\pm\sqrt{\lambda_{x}+\lambda_{y}}.$ (55)

As an example, we show the energy bands of the system obtained from the non-
Bloch band theory in Fig. 3(a) and confirm that they match the energy levels
in a finite open plane shown in blue in Fig. 3(b).

We note that in Fig. 3(b), one can see that the energy levels shown in magenta
do not belong to the bulk bands. These originate from the hybrid skin-
topological states [40, 43, 45, 46, 49], which result from the hybridization
between topological edge states and non-Hermitian skin states. A general
method of the calculation of energy bands of such states is beyond the present
work. On the other hand, in some models, such as the OTY model studied in Sec.
III.3, we can obtain the energy bands of the hybrid skin-topological states.

Finally, we investigate a topological phase of this system. Based on the
result obtained from Ref. [61], we defined the topological invariant as

$\nu=w_{x}w_{y},$ (56)

where

$w_{a}=-\frac{1}{2\pi}\frac{\left[{\rm
arg}R_{+}\left(\beta_{a}\right)\right]_{C_{\beta_{a}}}-\left[{\rm
arg}R_{-}\left(\beta_{a}\right)\right]_{C_{\beta_{a}}}}{2}$ (57)

for $a=x,y$. Here, let $C_{\beta_{a}}$ denote a loop formed by Eq. (53) on the
complex $\beta_{a}$-plane. In Eq. (57), $\left[{\rm
arg}R_{\pm}\left(\beta_{a}\right)\right]_{C_{\beta_{a}}}$ means the change of
${\rm arg}R_{\pm}\left(\beta_{a}\right)$ as $\beta_{a}$ goes along
$C_{\beta_{a}}$ in a counterclockwise way. In the definition, we use the fact
that the non-Bloch matrix can be written in the form of Eq. (46). Importantly,
the nonzero value of $\nu$ corresponds to existence of the topological corner
states at $E=0$. Indeed, in the present case, the system has $w_{x}=1$ and
$w_{y}=1$, which agrees with appearance of the zero-energy modes in a finite
open plane as shown in Fig. 3(b).

Figure 3: Energy eigenvalue of the non-Hermitian Benalcazar-Bernevig-Hughes
model. (a) Energy band calculated from the generalized Brillouin zone. (b)
Energy level in the system with $L_{x}=L_{y}=40$. In (b), the energy levels of
the hybrid topological-skin states and the topological corner states are shown
in magenta and orange, respectively. The system parameters are set to be
$t_{1}^{x}=1.2,t_{2}^{x}=1,t_{3}^{x}=0.2,\gamma^{x}=4/3,t_{1}^{y}=0.8,t_{2}^{y}=0.9,t_{3}^{y}=0.2$,
and $\gamma^{y}=5/3$.

###  III.3 Okugawa-Takahashi-Yokomizo model

Figure 4: Eigenstate and energy eigenvalue in the Okugawa-Takahashi-Yokomizo
model. Spatial distribution of (a) the hybrid skin-topological states and (b)
the extended state in the system with $L_{x}=L_{y}=30$. The energy levels of
the hybrid skin-topological states and the extended state are $E_{1}=0.84$ and
$E_{2}=1.45-2.1i$, respectively. (c) Energy level in the system with
$L_{x}=L_{y}=30$. (d) Energy band calculated from the non-Bloch band theory
proposed in this work (blue) and that for the hybrid skin-topological states
(red). We set the system parameters to be $t_{x}=1.2,g_{x}=0.9,t_{y}=0.8$, and
$g_{y}=0.7$.

Finally, we study the OTY model proposed in Ref. [42]. The real-space
Hamiltonian of this system is given by

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{\bm{r}}\left({\bm{c}}_{{\bm{r}}+\hat{\bm{x}}}^{\dagger}T_{x}^{-}{\bm{c}}_{\bm{r}}+{\bm{c}}_{\bm{n}}T_{x}^{+}{\bm{c}}_{{\bm{r}}+\hat{\bm{x}}}\right.$
(58)
$\displaystyle\left.+{\bm{c}}_{{\bm{r}}+\hat{\bm{y}}}^{\dagger}T_{y}^{-}{\bm{c}}_{\bm{r}}+{\bm{c}}_{\bm{r}}T_{y}^{+}{\bm{c}}_{{\bm{r}}+\hat{\bm{y}}}\right),$

where

$\displaystyle T_{x}^{\pm}=\left(\begin{array}[]{cc}t_{x}\pm
g_{x}&0\vspace{3pt}\\\ 0&t_{x}\mp
g_{x}\end{array}\right),~{}T_{y}^{\pm}=\left(\begin{array}[]{cc}0&-t_{y}\pm
g_{y}\vspace{3pt}\\\ t_{y}\pm g_{y}&0\end{array}\right),$ (63) (64)

and ${\bm{c}}_{\bm{r}}^{\dagger}=\left(c_{{\bm{r}},{\rm
A}}^{\dagger},c_{{\bm{r}},{\rm B}}^{\dagger}\right)$ are the creation
operators of particles on sublattices A and B, respectively. All the
parameters are set to be positive real numbers. We note that this system is in
the case of $N_{x}=1,N_{y}=1$, and $q=2$ in Eq. (1). The real-space eigen-
equation is obtained as

$T_{x}^{-}{\bm{\psi}}_{{\bm{r}}+\hat{\bm{x}}}+T_{x}^{+}{\bm{\psi}}_{{\bm{r}}-\hat{\bm{x}}}+T_{y}^{-}{\bm{\psi}}_{{\bm{r}}+\hat{\bm{y}}}+T_{y}^{+}{\bm{\psi}}_{{\bm{r}}-\hat{\bm{y}}}=E{\bm{\psi}}_{\bm{r}},$
(65)

where ${\bm{\psi}}_{\bm{n}}=\left(\psi_{{\bm{n}},{\rm A}},\psi_{{\bm{n}},{\rm
B}}\right)^{\rm T}$. The previous work proposed that the OTY model exhibits
the second-order non-Hermitian skin effect. In this effect,
$\mathcal{O}\left(L_{x}\right)$ (or $\mathcal{O}\left(L_{y}\right)$)
eigenstates are localized at the corners while
$\mathcal{O}\left(L_{x}L_{y}\right)$ eigenstates extend over the system. For
example, among the energy levels of the OTY model on a finite open plane with
$L_{x}=L_{y}=30$ as shown in Fig. 4(c), the eigenstates with $E_{1}=0.84$ are
localized at the corners [Fig. 4(a)] while the eigenstate with
$E_{2}=1.45-2.1i$ extends over the system [Fig. 4(b)]. We note that the
second-order non-Hermitian skin state in this model is a kind of the hybrid
skin-topological state [43].

Now, we calculate the energy bands of the OTY model from the non-Bloch band
theory. In this system, we have the characteristic equation of the Bloch
Hamiltonian as

$\displaystyle R_{x}\left(e^{ik_{x}},E\right)+R_{y}\left(e^{ik_{y}}\right)=0,$
(66)

where

$\displaystyle R_{x}\left(e^{ik_{x}},E\right)$ $\displaystyle=$
$\displaystyle\left[t_{x}\left(e^{ik_{x}}+e^{-ik_{x}}\right)-E\right]^{2}$
$\displaystyle-g_{x}^{2}\left(e^{ik_{x}}-e^{-ik_{x}}\right)^{2},$
$\displaystyle R_{y}\left(e^{ik_{y}}\right)$ $\displaystyle=$ $\displaystyle
t_{y}^{2}\left(e^{ik_{y}}+e^{-ik_{y}}\right)^{2}-g_{y}^{2}\left(e^{ik_{y}}-e^{-ik_{y}}\right)^{2}.$

Hence, this system is in Case B as discussed in Sec. II.2. Then, the
conditions for the generalized Brillouin zone are given by the solutions of
the quadratic equation $R_{x}\left(\beta_{x},E\right)=\lambda$ for $\beta_{x}$
and those of the quadratic equation $R_{y}\left(\beta_{y}\right)=-\lambda$ for
$\beta_{y}$. In fact, the condition is obtained as

$\left|\beta_{a,2}\right|=\left|\beta_{a,3}\right|,~{}\left(a=x,y\right),$
(68)

where the solutions satisfy Eq. (3). At last, we can calculate the energy
bands by combining the generalized Brillouin zone and Eq. (66). Importantly,
the generalized Brillouin zone spanned by
$\left(\beta_{x},\beta_{y}\right)~{}\left(\beta_{x}\equiv
e^{ik_{x}},\beta_{y}\equiv e^{ik_{y}}\right)$ becomes a torus defined by
$\left|\beta_{x}\right|=1$ and $\left|\beta_{y}\right|=1$ in this case because
the Bloch wave vector $\left(k_{x},k_{y}\right)$ takes real values. Then, we
obtain the energy bands as

$\displaystyle E\left(k_{x},k_{y}\right)$ $\displaystyle=$ $\displaystyle
2t_{x}\cos k_{x}$ $\displaystyle\pm
2i\sqrt{g_{x}^{2}\sin^{2}k_{x}+t_{y}^{2}\cos^{2}k_{y}+g_{y}^{2}\sin^{2}k_{y}}.$

The results are shown in blue in Fig. 4(d). Comparing Fig. 4(c) with Fig.
4(d), we see that the blue energy bands reproduces the bulk energy levels in a
finite open plane. Meanwhile, the energy levels of the hybrid skin-topological
states are not obtained from the non-Bloch band theory proposed in this work.

Next, we explain a way to calculate the energy band of the hybrid skin-
topological state. In the following, we set $t_{y}>g_{y}$ for convenience. The
key ingredient is that $\beta_{y}$ can be determined as a constant so that the
non-Bloch matrix

${\cal
H}\left(\beta_{x},\beta_{y}\right)=\beta_{x}T_{x}^{-}+\beta_{x}^{-1}T_{x}^{+}+\beta_{y}T_{y}^{-}+\beta_{y}^{-1}T_{y}^{+}$
(70)

has eigenvectors independent of the value of $\beta_{x}$, such as
$\left(1,0\right)^{\rm T}$ and $\left(0,1\right)^{\rm T}$. In this case, a
non-Bloch wave for $\beta_{x}$ can be constructed. In fact, we can take an
ansatz

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle{\bm{\psi}}_{\bm{r}}=\Psi_{n_{x}}{\bm{\Phi}}_{n_{y}},\vspace{5pt}\\\
\displaystyle\Psi_{n_{x}}=\sum_{j}\left(\beta_{x,j}\right)^{n_{x}}\upsilon^{\left(j\right)},\vspace{5pt}\\\
\displaystyle{\bm{\Phi}}_{n_{y}}=\sum_{j}\left(\beta_{y,j}\right)^{n_{y}}\chi^{\left(j\right)}\left(\begin{array}[]{c}0\\\
1\end{array}\right)\end{array}\right.$ (76)

for Eq. (65), where $\upsilon^{\left(j\right)}$ and $\chi^{\left(j\right)}$
are constants. By imposing $\left(0,1\right)^{\rm T}$ to be an eigenvector of
${\cal H}\left(\beta_{x},\beta_{y}\right)$, the $\left(1,2\right)$ component
of ${\cal H}\left(\beta_{x},\beta_{y}\right)$ should vanish, and possible
values of $\beta_{y}$ are given by

$\beta_{y,1}=-i\sqrt{\frac{t_{y}-g_{y}}{t_{y}+g_{y}}},~{}\beta_{y,2}=i\sqrt{\frac{t_{y}-g_{y}}{t_{y}+g_{y}}}.$
(77)

In this case, by putting $\chi^{\left(1\right)}=-\chi^{\left(2\right)}$, the
wave function satisfies an open boundary condition along the $y$ direction

${\bm{\Phi}}_{0}={\bm{0}}.$ (78)

We note that the wave function is located near $n_{y}=1$ because
$t_{y}>g_{y}>0$. Then, the energy eigenvalue of ${\cal
H}\left(\beta_{x},\beta_{y}\right)$ is given by

$\left(t_{x}+g_{x}\right)\beta_{x}+\frac{t_{x}-g_{x}}{\beta_{x}}=E,$ (79)

Therefore, the condition for the generalized Brillouin zone for the hybrid
skin-topological state is obtained from the two solutions of Eq. (79) as

$\left|\beta_{x,1}\right|=\left|\beta_{x,2}\right|.$ (80)

Finally, we can get the energy band as

$E=2\sqrt{t_{x}^{2}-g_{x}^{2}}\cos\theta,$ (81)

where $\theta$ is real.

Importantly, we can take another ansatz of the hybrid skin-topological state

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle{\bm{\psi}}_{\bm{r}}=\Psi_{n_{x}}{\bm{\Phi}}_{n_{y}},\vspace{5pt}\\\
\displaystyle\Psi_{n_{x}}=\sum_{j}\left(\beta_{x,j}\right)^{n_{x}}\bar{\upsilon}^{\left(j\right)},\vspace{5pt}\\\
\displaystyle{\bm{\Phi}}_{n_{y}}=\sum_{j}\left(\beta_{y,j}\right)^{n_{y}-L_{y}-1}\bar{\chi}^{\left(j\right)}\left(\begin{array}[]{c}1\\\
0\end{array}\right).\end{array}\right.$ (87)

Here, $\bar{\upsilon}^{\left(j\right)}$ and $\bar{\chi}^{\left(j\right)}$ are
constants. By imposing $\left(1,0\right)^{\rm T}$ to be an eigenvector of
${\cal H}\left(\beta_{x},\beta_{y}\right)$, the allowed values of $\beta_{y}$
are given by

$\beta_{y,1}=-i\sqrt{\frac{t_{y}+g_{y}}{t_{y}-g_{y}}},~{}\beta_{y,2}=i\sqrt{\frac{t_{y}+g_{y}}{t_{y}-g_{y}}},$
(88)

In this case, by taking
$\bar{\chi}^{\left(1\right)}=-\bar{\chi}^{\left(2\right)}$, the wave function
satisfies an open boundary condition along the $y$ direction

${\bm{\Phi}}_{L_{y}+1}={\bm{0}},$ (89)

and it is located near $n_{y}=L_{y}$. Finally, the energy eigenvalue of ${\cal
H}\left(\beta_{x},\beta_{y}\right)$ is given by

$\left(t_{x}-g_{x}\right)\beta_{x}+\frac{t_{x}+g_{x}}{\beta_{x}}=E,$ (90)

and we can get the energy band of this hybrid skin-topological state in the
form of Eq. (81).

The energy band for the hybrid skin-topological state obtained above is shown
in red in Fig. 4(d). Comparing Fig. 4(c) with Fig. 4(d), we confirm that our
calculation reproduces the energy levels of the hybrid skin-topological states
in a finite open plain.

##  IV Discussion

Figure 5: Reflection of the plane waves at the boundaries in the two-
dimensional non-Hermitian system. Through the reflection at the boundary
parallel to the $y$ direction, the wave number $k_{x}$ changes to
$k_{x}^{\prime}$, but the wave number $k_{y}$ is unchanged. The situation is
similar at the boundary parallel to the $x$ direction.

In this section, we explain the reason why our non-Bloch band theory is
established in the 2D non-Hermitian system. The key ingredient of our work is
that the ansatz of the real-space eigen-equation can be written as Eqs. (5)
and (11) in Case A and Case B, respectively. This is related to the condition
for the standing-wave formation. In the bulk, the standing wave is formed by
interference between plane waves which are generated by reflection at the
boundaries of the system. In general, the process of the standing-wave
formation is complex. In some cases, the condition for the standing-wave
formation becomes straightforward. In Fig. 5, we show such a case in the
rectangle geometry, where the standing wave consists of four plane waves. The
plane wave $e^{ik_{x}x+ik_{y}y}$ preserves the wave number in the direction
parallel to the boundary where the plane wave is reflected. In the present
work, we propose the two cases in which the above situation can be realized.
In Case A, the reflection of the plane wave becomes simple because of the
suppression of the non-Hermitian skin effect in one direction. In Case B,
since the motions of the plane wave in the $x$ direction and the $y$ direction
are decoupled, we can regard the 2D system as two 1D systems. Finally, the
leading terms of the standing wave have the form of the linear combination of
a few plane waves in both cases.

Based on the above concept, it naturally follows that the geometry of the
system largely affects the condition for the standing-wave formation. In fact,
the non-Hermitian skin effect in a 2D system depends on geometry [44] Namely,
the non-Bloch band theory proposed in this work should depend crucially on
geometries. For example, in the model of the non-Hermitian Chern insulator as
discussed in Sec. III.1, the energy bands obtained from the non-Bloch band
theory for the rectangle geometry [Fig. 6(a), which is the same as Fig. 2(b)]
do not match the energy levels in a finite open plane [Fig. 6(b)] on the
diamond geometry [Fig. 6(c)]. This indicates that non-Bloch bands on the
rectangle geometry are different from those on the diamond geometry. Thus, we
conclude that the energy bands in a 2D non-Hermitian system depend on
geometries. Construction of a general method to calculate the energy bands in
a 2D non-Hermitian system with various geometries is left for a future work.

Finally, we focus on the model of the non-Hermitian Chern insulator with
$\gamma_{x}\neq 0$ and $\gamma_{y}\neq 0$. We note that the non-Hermitian skin
effect occurs in both the $x$ and the $y$ directions, and the system is not in
Case A. In this case, the plane wave is scattered at the boundary into
numerous waves. Since the reflection process of the plane wave becomes complex
compared to the case in Fig. 5, we cannot find the condition for the standing-
wave formation, and the standing wave cannot be written in the form of Eq.
(5). Hence, it is unclear whether we can construct the non-Bloch band theory
in the model of the non-Hermitian Chern insulator with arbitrary parameters,
and also in general 2D non-Hermitian systems.

Figure 6: Geometry dependence of the energy bands in the model of the non-
Hermitian Chern insulator (5). (a) Energy band calculated from the non-Bloch
band theory for the rectangle geometry proposed in this work. (b) Energy level
in a finite open plane on the diamond geometry with $L=59$. (c) Diamond
geometry. In (a) and (b), we set the system parameters same as Figs. 2(a)-(c).

##  V Summary

In this paper, we propose two classes of the 2D non-Hermitian tight-binding
systems in which the non-Bloch band theory can be constructed. In terms of the
present non-Bloch band theory, we calculate the generalized Brillouin zone and
the energy bands in the model of the non-Hermitian Chern insulator, the non-
Hermitian BBH model, and the OTY model. We confirm that the energy bands
reproduce the bulk energy levels in a finite open plane on the rectangle
geometry. Furthermore, we investigate topological properties of theses models
from the non-Bloch band theory. In the model of the non-Hermitian Chern
insulator and the non-Hermitian BBH model, it is shown that the topological
invariant defined from the generalized Brillouin zone can predict existence of
the topological edge states. In the OTY model, we propose the method to
calculate the energy bands for the hybrid skin-topological state. We note that
the construction of the non-Bloch band theory in general 2D non-Hermitian
system and a general method to calculate the energy bands for the hybrid skin-
topological state are left for a future work.

###### Acknowledgements.

We are grateful to Ryo Okugawa for variable discussion. This work was
supported by JSPS KAKENHI Grant No. 22H00108 and No. JP21J01409.

##  Appendix A Non-Bloch band theory in the case with $q=1$

In this appendix, we construct the non-Bloch band theory in the system
described by the Hamiltonian (1) with $q=1$. For example, the non-Hermitian
skin effect in this case was investigated in a synthetic photonic lattice
[62]. When $q=1$, the real-space eigen-equation (2) is explicitly written as

$\sum_{i=-N_{x}}^{N_{x}}t_{i}^{\left(x\right)}\psi_{{\bm{r}}+i\hat{\bm{x}}}+\sum_{i=-N_{y}}^{N_{y}}t_{i}^{\left(y\right)}\psi_{{\bm{r}}+i\hat{\bm{y}}}=E\psi_{\bm{r}},$
(91)

where ${\bm{r}}=\left(n_{x},n_{y}\right)$. When we take

$\psi_{\bm{r}}=\varphi_{n_{x}}^{\left(x\right)}\varphi_{n_{y}}^{\left(y\right)}$
(92)

as an ansatz for Eq. (91), we get

$\displaystyle\sum_{i=-N_{a}}^{N_{a}}t_{i}^{\left(a\right)}\varphi_{n_{a}+i}^{\left(a\right)}$
$\displaystyle=$
$\displaystyle\lambda_{a}\varphi_{n_{a}}^{\left(a\right)},~{}\left(a=x,y\right),$
(93) $\displaystyle\lambda_{x}+\lambda_{y}$ $\displaystyle=$ $\displaystyle
E.$ (94)

Here, the solutions of Eq. (93) can be given in the form of a linear
combination. Namely, we obtain

$\varphi_{n_{a}}^{\left(a\right)}=\sum_{j}\left(\beta_{a,j}\right)^{n_{a}}\phi_{j},$
(95)

where $\beta_{a}=\beta_{a,j}$ is the solution of the characteristic equation

$\sum_{m=-N_{a}}^{N_{a}}t_{i}^{\left(a\right)}\left(\beta_{a}\right)^{m}-\lambda_{a}=0$
(96)

for $a=x,y$. By following the conventional non-Bloch band theory, we can
obtain the condition for the generalized Brillouin zone in terms of the
solutions of Eq. (96), given by

$\left|\beta_{a,N_{a}}\right|=\left|\beta_{a,N_{a}+1}\right|,$ (97)

where

$\left|\beta_{a,i}\right|\leq\dots\left|\beta_{a,2N_{a}}\right|$ (98)

for $a=x,y$. Thus, a set of the trajectories of $\beta_{a,N_{a}}$ and
$\beta_{a,N_{a}+1}$ satisfying Eq. (97) forms the generalized Brillouin zone.
Finally, from $\beta_{x}$ and $\beta_{y}$ on the generalized Brillouin zone,
we can calculate the energy band by using Eq. (94).

##  Appendix B Quantization of the Chern number

Figure 7: Generalized Brillouin zone on the $k_{x}\mathchar 45\relax k_{y}$
plane for the complex Bloch wave vector ${\bm{k}}$. The imaginary part of
${\bm{k}}$ is omitted. $l$ represents the boundary between region I and region
II.

In this appendix, we show that the values of the Chern number defined in Eq.
(27) are quantized when the band $n$ is separated from other bands by a line
gap. We assume that the system is in Case A or Case B, and the generalized
Brillouin zone can be constructed. For the non-Bloch matrix obtained from Eq.
(1), let $\left|u_{R,n}\left({\bm{k}}\right)\right\rangle$ and $\left\langle
u_{L,n}\left({\bm{k}}\right)\right|$ denote the right and left eigenvectors,
respectively. Here, ${\bm{k}}$ is the complex Bloch wave vector, and $n$
expresses the band index. We impose that the right and left eigenvectors
satisfy

$\left\langle
u_{L,m}\left({\bm{k}}\right)\middle|u_{R,n}\left({\bm{k}}\right)\right\rangle=\delta_{m,n}.$
(99)

Furthermore, we define the Berry connection as

$\mathcal{\bm{A}}_{n}\left({\bm{k}}\right)=\left\langle
u_{L,n}\left({\bm{k}}\right)\middle|{\bm{\nabla}}_{\bm{k}}\middle|u_{R,n}\left({\bm{k}}\right)\right\rangle.$
(100)

First of all, we discuss the generalized Brillouin zone in this work. In our
2D non-Hermitian system, let us regard $k_{y}$ as a system parameter. Under a
given value of $k_{y}$, the system is regarded as a 1D non-Hermitian system
described by the complex Bloch wave number $k_{x}$. Then, a complex parameter
$e^{ik_{x}}$ forms a loop on the complex plane [63]. This leads to numerous
loops of $e^{ik_{x}}$ for all values of $k_{y}$. Similarly,
$e^{ik_{y}},~{}k_{y}\in{\mathbb{C}}$ forms a loop for all values of $k_{x}$.
Therefore, the whole generalized Brillouin zone in the 2D non-Hermitian system
forms a closed surface. For example, the generalized Brillouin zone in the
Chern number discussed in Sec. II.1 becomes a closed surface [Fig. 2(a)].

Topologically nontrivial phases can be associated with topological obstruction
for a single gauge choice over the whole ${\bm{k}}$ space. Therefore, in the
present case, we separate the generalized Brillouin zone into two regions so
that we take gauge I on one region and gauge II on the other region [Fig. 7].
Importantly, the two gauges are related by the gauge transformation

$\displaystyle\left|\psi_{R,n}^{\rm II}\left({\bm{k}}\right)\right\rangle$
$\displaystyle=$ $\displaystyle
r_{n}\left({\bm{k}}\right)e^{i\chi_{n}\left({\bm{k}}\right)}\left|\psi_{R,n}^{\rm
I}\left({\bm{k}}\right)\right\rangle,$ (101)
$\displaystyle\left\langle\psi_{L,n}^{\rm II}\left({\bm{k}}\right)\right|$
$\displaystyle=$
$\displaystyle\frac{1}{r_{n}\left({\bm{k}}\right)}e^{-i\chi_{n}\left({\bm{k}}\right)}\left\langle\psi_{L,n}^{\rm
I}\left({\bm{k}}\right)\right|,$ (102)

where $r_{n}\left({\bm{k}}\right)$ and $\chi_{n}\left({\bm{k}}\right)$ are
real functions [25]. Here, from the Stokes’ theorem, Eq. (27) can be rewritten
in terms of the corresponding Berry connection as

$C_{n}=\frac{1}{2\pi i}\oint_{l}d{\bm{k}}\cdot\left(\mathcal{\bm{A}}^{\rm
I}_{n}\left({\bm{k}}\right)-\mathcal{\bm{A}}^{\rm
II}_{n}\left({\bm{k}}\right)\right),$ (103)

where $l$ denote the boundary between region I and region II. Hence, we obtain

$C_{n}=\frac{1}{2\pi}\oint_{l}d{\bm{k}}\cdot{\bm{\nabla}}_{\bm{k}}\chi_{n}\left({\bm{k}}\right),$
(104)

where we use

$\mathcal{\bm{A}}^{\rm II}_{n}\left({\bm{k}}\right)=\mathcal{\bm{A}}^{\rm
I}_{n}\left({\bm{k}}\right)+i{\bm{\nabla}}_{\bm{k}}\chi_{n}\left({\bm{k}}\right)+\frac{{\bm{\nabla}}_{\bm{k}}r_{n}\left({\bm{k}}\right)}{r_{n}\left({\bm{k}}\right)}$
(105)

and

$\oint_{l}d{\bm{k}}\cdot{\bm{\nabla}}_{\bm{k}}\log
r_{n}\left({\bm{k}}\right)=0.$ (106)

The integral in Eq. (104) is always equal to an integer multiple of $2\pi$
because the wave functions are single-valued. As a result, the Chern number
defined in Eq. (27) takes quantized values when the band $n$ is separated from
other bands by a line gap.

## References

  * Ashida  _et al._ [2020] Y. Ashida, Z. Gong, and M. Ueda, Adv. Phys. 69, 249 (2020).
  * Yao and Wang [2018] S. Yao and Z. Wang, Phys. Rev. Lett. 121, 086803 (2018).
  * Song  _et al._ [2019] F. Song, S. Yao, and Z. Wang, Phys. Rev. Lett. 123, 170401 (2019).
  * Okuma  _et al._ [2020] N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, Phys. Rev. Lett. 124, 086801 (2020).
  * Zhang  _et al._ [2020] K. Zhang, Z. Yang, and C. Fang, Phys. Rev. Lett. 125, 126402 (2020).
  * Li  _et al._ [2020] L. Li, C. H. Lee, S. Mu, and J. Gong, Nat. Commun. 11, 5491 (2020).
  * Yokomizo and Murakami [2021a] K. Yokomizo and S. Murakami, Phys. Rev. B 104, 165117 (2021a).
  * Brandenbourger  _et al._ [2019] M. Brandenbourger, X. Locsin, E. Lerner, and C. Coulais, Nat. Commun. 10, 4608 (2019).
  * Xiao  _et al._ [2020] L. Xiao, T. Deng, K. Wang, G. Zhu, Z. Wang, W. Yi, and P. Xue, Nat. Phys. 16, 761 (2020).
  * Weidemann  _et al._ [2020] S. Weidemann, M. Kremer, T. Helbig, T. Hofmann, A. Stegmaier, M. Greiter, R. Thomale, and A. Szameit, Science 368, 311 (2020).
  * Helbig  _et al._ [2020] T. Helbig, T. Hofmann, S. Imhof, M. Abdelghany, T. Kiessling, L. Molenkamp, C. Lee, A. Szameit, M. Greiter, and R. Thomale, Nat. Phys. 16, 747 (2020).
  * Hofmann  _et al._ [2020] T. Hofmann, T. Helbig, F. Schindler, N. Salgo, M. Brzezińska, M. Greiter, T. Kiessling, D. Wolf, A. Vollhardt, A. Kabaši, C. H. Lee, A. Bilušić, R. Thomale, and T. Neupert, Phys. Rev. Research 2, 023265 (2020).
  * Ghatak  _et al._ [2020] A. Ghatak, M. Brandenbourger, J. van Wezel, and C. Coulais, Proc. Nat. Ac. Sc. USA 117, 29561 (2020).
  * Zhang  _et al._ [2021a] X. Zhang, Y. Tian, J.-H. Jiang, M.-H. Lu, and Y.-F. Chen, Nat. Commun. 12, 5377 (2021a).
  * Chen  _et al._ [2021] Y. Chen, X. Li, C. Scheibner, V. Vitelli, and G. Huang, Nat. Commun. 12, 5935 (2021).
  * Zhang  _et al._ [2021b] L. Zhang, Y. Yang, Y. Ge, Y.-J. Guan, Q. Chen, Q. Yan, F. Chen, R. Xi, Y. Li, D. Jia, _et al._ , Nat. Commun. 12, 6297 (2021b).
  * Wang  _et al._ [2022] W. Wang, X. Wang, and G. Ma, Nature 608, 50 (2022).
  * Liang  _et al._ [2022] Q. Liang, D. Xie, Z. Dong, H. Li, H. Li, B. Gadway, W. Yi, and B. Yan, Phys. Rev. Lett. 129, 070401 (2022).
  * Yokomizo and Murakami [2019] K. Yokomizo and S. Murakami, Phys. Rev. Lett. 123, 066404 (2019).
  * Yokomizo and Murakami [2020a] K. Yokomizo and S. Murakami, Prog. Theor. Exp. Phys. 2020, 12A102 (2020a).
  * Yokomizo and Murakami [2021b] K. Yokomizo and S. Murakami, Phys. Rev. B 103, 165123 (2021b).
  * Yokomizo  _et al._ [2022] K. Yokomizo, T. Yoda, and S. Murakami, Phys. Rev. Research 4, 023089 (2022).
  * Yokomizo and Murakami [2020b] K. Yokomizo and S. Murakami, Phys. Rev. Research 2, 043045 (2020b).
  * Leykam  _et al._ [2017] D. Leykam, K. Y. Bliokh, C. Huang, Y. D. Chong, and F. Nori, Phys. Rev. Lett. 118, 040401 (2017).
  * Shen  _et al._ [2018] H. Shen, B. Zhen, and L. Fu, Phys. Rev. Lett. 120, 146402 (2018).
  * Kawabata  _et al._ [2018] K. Kawabata, K. Shiozaki, and M. Ueda, Phys. Rev. B 98, 165148 (2018).
  * Ezawa [2019] M. Ezawa, Phys. Rev. B 100, 081401 (2019).
  * Luo and Zhang [2019] X.-W. Luo and C. Zhang, Phys. Rev. Lett. 123, 073601 (2019).
  * Wu  _et al._ [2019] H. C. Wu, L. Jin, and Z. Song, Phys. Rev. B 100, 155117 (2019).
  * Wu  _et al._ [2020] Y.-J. Wu, C.-C. Liu, and J. Hou, Phys. Rev. A 101, 043833 (2020).
  * Xue  _et al._ [2020] H. Xue, Q. Wang, B. Zhang, and Y. D. Chong, Phys. Rev. Lett. 124, 236403 (2020).
  * Ao  _et al._ [2020] Y. Ao, X. Hu, Y. You, C. Lu, Y. Fu, X. Wang, and Q. Gong, Phys. Rev. Lett. 125, 013902 (2020).
  * Song  _et al._ [2020a] A. Y. Song, X.-Q. Sun, A. Dutt, M. Minkov, C. Wojcik, H. Wang, I. A. D. Williamson, M. Orenstein, and S. Fan, Phys. Rev. Lett. 125, 033603 (2020a).
  * Xie  _et al._ [2021] L. C. Xie, H. C. Wu, X. Z. Zhang, L. Jin, and Z. Song, Phys. Rev. B 104, 125406 (2021).
  * Teo  _et al._ [2022] H. T. Teo, H. Xue, and B. Zhang, Phys. Rev. A 105, 053510 (2022).
  * Philip  _et al._ [2018] T. M. Philip, M. R. Hirsbrunner, and M. J. Gilbert, Phys. Rev. B 98, 155430 (2018).
  * Chen and Zhai [2018] Y. Chen and H. Zhai, Phys. Rev. B 98, 245130 (2018).
  * Hirsbrunner  _et al._ [2019] M. R. Hirsbrunner, T. M. Philip, and M. J. Gilbert, Phys. Rev. B 100, 081104 (2019).
  * Groenendijk  _et al._ [2021] S. Groenendijk, T. L. Schmidt, and T. Meng, Phys. Rev. Research 3, 023001 (2021).
  * Lee  _et al._ [2019] C. H. Lee, L. Li, and J. Gong, Phys. Rev. Lett. 123, 016805 (2019).
  * Kawabata  _et al._ [2020] K. Kawabata, M. Sato, and K. Shiozaki, Phys. Rev. B 102, 205118 (2020).
  * Okugawa  _et al._ [2020] R. Okugawa, R. Takahashi, and K. Yokomizo, Phys. Rev. B 102, 241202 (2020).
  * Fu  _et al._ [2021] Y. Fu, J. Hu, and S. Wan, Phys. Rev. B 103, 045420 (2021).
  * Zhang  _et al._ [2022] K. Zhang, Z. Yang, and C. Fang, Nat. Commun. 13, 2496 (2022).
  * Li  _et al._ [2022] Y. Li, C. Liang, C. Wang, C. Lu, and Y.-C. Liu, Phys. Rev. Lett. 128, 223903 (2022).
  * Zhu and Gong [2022] W. Zhu and J. Gong, Phys. Rev. B 106, 035425 (2022).
  * Fang  _et al._ [2022] Z. Fang, M. Hu, L. Zhou, and K. Ding, Nanophotonics 11, 3447 (2022).
  * Palacios  _et al._ [2021] L. S. Palacios, S. Tchoumakov, M. Guix, I. Pagonabarraga, S. Sánchez, and A. G Grushin, Nat. Commun. 12, 4691 (2021).
  * Zou  _et al._ [2021] D. Zou, T. Chen, W. He, J. Bao, C. H. Lee, H. Sun, and X. Zhang, Nat. Commun. 12, 7201 (2021).
  * [50] C. Shang, S. Liu, R. Shao, P. Han, X. Zang, X. Zhang, K. N. Salama, W. Gao, C. H. Lee, R. Thomale, A. Manchon, S. Zhang, T. J. Cui, and U. Schwingenschlogl, arXiv:2203.00484 .
  * Yao  _et al._ [2018] S. Yao, F. Song, and Z. Wang, Phys. Rev. Lett. 121, 136802 (2018).
  * Liu  _et al._ [2019] T. Liu, Y.-R. Zhang, Q. Ai, Z. Gong, K. Kawabata, M. Ueda, and F. Nori, Phys. Rev. Lett. 122, 076801 (2019).
  * Yu  _et al._ [2021] Y. Yu, M. Jung, and G. Shvets, Phys. Rev. B 103, L041102 (2021).
  * Lin  _et al._ [2021] Z. Lin, L. Ding, S. Chen, S. Li, S. Ke, X. Li, and B. Wang, Phys. Rev. A 103, 063507 (2021).
  * Xiao and Chan [2022] Y.-X. Xiao and C. T. Chan, Phys. Rev. B 105, 075128 (2022).
  * [56] J. Bartlett and E. Zhao, arXiv:2209.06774 .
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Phys. Rev. X 9, 041015 (2019).
  * Yi and Yang [2020] Y. Yi and Z. Yang, Phys. Rev. Lett. 125, 186802 (2020).
  * Okugawa  _et al._ [2021] R. Okugawa, R. Takahashi, and K. Yokomizo, Phys. Rev. B 103, 205205 (2021).
  * Fukui  _et al._ [2005] T. Fukui, Y. Hatsugai, and H. Suzuki, J. Phys. Soc. Jpn. 74, 1674 (2005).
  * Okugawa  _et al._ [2019] R. Okugawa, S. Hayashi, and T. Nakanishi, Phys. Rev. B 100, 235302 (2019).
  * Song  _et al._ [2020b] Y. Song, W. Liu, L. Zheng, Y. Zhang, B. Wang, and P. Lu, Phys. Rev. Applied 14, 064076 (2020b).
  * Yang  _et al._ [2020] Z. Yang, K. Zhang, C. Fang, and J. Hu, Phys. Rev. Lett. 125, 226402 (2020).


11institutetext: Department of Electrical Engineering and Electronics,
University of Liverpool, Liverpool L69 3GJ, UK; 11email: {E.Drousiotis,
A.M.Philips<EMAIL_ADDRESS>22institutetext: Department of
Computer Science, University of Liverpool, Liverpool L69 3BX, UK
22email<EMAIL_ADDRESS>

# Bayesian Decision Trees Inspired from Evolutionary Algorithms

Efthyvoulos Drousiotis 11 Alexander M  Phillips 11 Paul G  Spirakis 22 Simon
Maskell 11

###### Abstract

Bayesian Decision Trees (DTs) are generally considered a more advanced and
accurate model than a regular Decision Tree (DT) because they can handle
complex and uncertain data. Existing work on Bayesian DTs uses Markov Chain
Monte Carlo (MCMC) with an accept-reject mechanism and sample using naive
proposals to proceed to the next iteration, which can be slow because of the
burn-in time needed. We can reduce the burn-in period by proposing a more
sophisticated way of sampling or by designing a different numerical Bayesian
approach. In this paper, we propose a replacement of the MCMC with an
inherently parallel algorithm, the Sequential Monte Carlo (SMC), and a more
effective sampling strategy inspired by the Evolutionary Algorithms (EA).
Experiments show that SMC combined with the EA can produce more accurate
results compared to MCMC in 100 times fewer iterations.

###### Keywords:

Swarm Particle Optimisation Bayesian Decision Trees Machine Learning

##  1 Introduction and Relevant Work

In Bayesian statistics, obtaining and calculating random samples from a
probability distribution is challenging. Markov Chain Monte Carlo (MCMC) is a
widely used method to tackle this issue. MCMC can characterise a distribution
without knowing its analytic form by using a series of samples. MCMC has been
used to solve problems in various domains, including biology[24],
forensics[23], education[10, 9], and chemistry[14], among other areas. Monte
Carlo applications are generally considered embarrassingly parallel since each
chain can run independently on two or more independent machines or cores.
Nevertheless, the principal problem is that processing within each chain is
not embarrassingly parallel. When the feature space and the proposal are
computationally expensive, we can only do a little to improve the running
time. As mentioned before, running multiple independent MCMC chains in
parallel can be a strategy for parallelization, but it is not very effective.
Even though you increase the number of chains and decrease the number of
samples produced from each chain, the burn-in process for each chain would
remain unchanged.

Much research has been done to improve the efficiency of MCMC methods. The
improvements can be grouped into several categories[19]. One of these
categories is based on understanding the geometry of the target density
function. An example of this is Hamiltonian Monte Carlo[13] (HMC), which uses
an auxiliary variable called momentum and updates it using the gradient of the
density function. Different methods, such as symplectic integrators of various
precision levels, have been developed to approximate the Hamiltonian
equation[3]. HMC tends to generate less correlated samples than the
Metropolis-Hastings algorithm, but it requires the gradient of the density
function to be available and computationally feasible.

Another approach involves dividing complex problems into smaller and more
manageable components. For example, as discussed earlier, using multiple
parallel MCMC chains to explore the parameter space and then combining the
samples obtained from these chains[18]. However, this approach does not
achieve faster convergence of the chains to the stationary distribution. This
is because all chains have to converge independently of each other in MCMC. It
is also possible to partition the data space, which is already implemented in
the context of Bayesian DTs[12], or partition the parameter space [1]into
simpler pieces that can process independently which are proven not to be much
effective.

MCMC DTs simulations can take a long to converge on a good model when the
state space is large and complex. This is due to both the number of iterations
needed and the complexity of the calculations performed in each iteration,
such as finding the best features and structure of a DT. MCMC, in general, and
as will be explained detailed in section 2, generates samples from a specific
distribution by proposing new samples and deciding whether to accept them
based on the evaluation of the posterior distribution. Because the current
sample determines the next step of the MCMC, it can be not easy to process a
single MCMC chain simultaneously using multiple processing elements. A method
is [11] proposed to parallelise a single chain on MCMC Decision trees.
However, the speedup is not always guaranteed.

Another method of speeding up MCMC, focuses on enhancing the proposal
function, which is the approach we pursue in this paper. This can be achieved
through techniques such as simulated tempering[17], adaptive MCMC[8], or
multi-proposal MCMC[16]. Having a good proposal function in Markov Chain Monte
Carlo (MCMC) and, in general, Monte Carlo methods is crucial for the
efficiency and accuracy of the algorithm. Poor proposals can lead to slow
convergence, poor mixing, and biased estimates. Good proposal functions can
efficiently explore the target distribution and reduce the correlation between
successive samples.

Several papers describe novelties specifically on Bayesian DTs, focusing on
different improvements. For example, [12, 11, 22], contributed towards the
runtime enhancement, [27] explored a new novel proposal move called "radical
restructure," which changes the structure of the tree ($T$), without changing
the number of leaves or the partition of observations into leaves. Moreover,
[4] proposed different criteria for accepting or rejecting the $T$ on MCMC,
such as posterior probability, marginal likelihood, residual sum of squares,
and misclassification rates. The general approach to improving the proposal
function in Monte Carlo methods has been introduced previously. However, in
the context of Bayesian DTs, there is still enough space for exploration.

Our contribution is then:

  * •

We describe for the first time a novel algorithm inspired by the Evolutionary
Algorithms to improve the proposal function of the Bayesian DTs that uses an
inherently parallel algorithm, a Sequential Monte Carlo (SMC) sampler, to
generate samples.

##  2 Bayesian Decision Trees

A DT operates by descending a tree $T$. The process of outputting a
classification probability for a given datum starts at a root node (see Figure
1). At each non-leaf node, a decision as to which child node to progress to is
made based on the datum and the parameters of the node. This process continues
until a leaf node is reached. At the leaf node, a node-specific and datum-
independent classification output is generated.

Figure 1: Decision Tree

Our model describes the conditional distribution of a value for $Y$ given the
corresponding values for $x$, where $x$ is a vector of predictors and $Y$ the
corresponding values that we predict. We define the tree to be $T$ such that
the function of the non-leaf nodes is to (implicitly) identify a region, $A$,
of values for $x$ for which $p\left(Y|x\in A\right)$ can be approximated as
independent of the specific value of $x\in A$, ie as $p(Y|x)\approx
p(Y|\phi_{j},x\in A)$. This model is called a probabilistic classification
tree, according to the quantitative response $Y$.

For a given tree, $T$, we define its depth to be $d(T)$, the set of leaf nodes
to be $L(T)$ and the set of non-leaf nodes to be $\bar{L}(T)$. The $T$ is then
parameterised by: the set of features for all non-leaf nodes,
$k_{\bar{L}(T)}$; the vector of corresponding thresholds, $c_{\bar{L}(T)}$;
the parameters, $\phi_{L(T)}$, of the conditional probabilities associated
with the leaf nodes. This is such that the parameters of the $T$ are
$\theta(T)=[k_{\bar{L}(T)},c_{\bar{L}(T)},\phi_{L(T)}]$ and $\theta(T)_{j}$
are the parameters associated with the $j$th node of the $T$, where:

$\displaystyle\theta(T)_{j}=\left\\{\begin{array}[]{ll}\left[k_{j},c_{j}\right]&j\in\bar{L}(T)\\\
\phi_{j}&j\in L(T).\end{array}\right.$ (3)

Given a dataset comprising $N$ data, $Y_{1:N}$ and corresponding features,
$x_{1:N}$, and since DTs are specified by $T$ and $\theta(T)$, a Bayesian
analysis of the problem proceeds by specifying a prior probability
distribution, $p(\theta(T),T)$ and associated likelihood,
$p(Y_{1:N}|T,\theta(T),x_{1:N})$. Because $\theta(T)$ defines the parametric
model for $T$, it will usually be convenient to adopt the following structure
for the joint probability distribution of $N$ data, $Y_{1:N}$, and the $N$
corresponding vectors of predictors, $x_{1:N}$:

$\displaystyle p(Y_{1:N},T,\theta(T)|x_{1:N})=$ $\displaystyle
p(Y_{1:N}|T,\theta(T),x_{1:N})p(\theta(T),T)$ (4) $\displaystyle=$
$\displaystyle p(Y_{1:N}|T,\theta(T),x_{1:N})p(\theta(T)|T)p(T)$ (5)

which we note is proportional to the posterior,
$p(T,\theta(T)|Y_{1:N},x_{1:N})$, and where we assume

$\displaystyle p(Y_{1:N}|T,\theta(T),x_{1:N})$
$\displaystyle=\prod_{i=1}^{N}p(Y_{i}|x_{i},T,\theta(T))$ (6) $\displaystyle
p(\theta(T)|T)$ $\displaystyle=\prod_{j\in T}p(\theta(T)_{j}|T)$ (7)
$\displaystyle=\prod_{j\in T}p(k_{j}|T)p(c_{j}|k_{j},T)$ (8) $\displaystyle
p(T)$ $\displaystyle=\frac{a}{(1+d(T))^{\beta}}$ (9)

Equation 6 describes the product of the probabilities of every data point,
$Y_{i}$, being classified correctly given the datum’s features, $x_{i}$, the
$T$ structure, and the features/thresholds, $\theta(T)$, associated with each
node of the $T$. At the $j$th node, equation 8 describes the product of
possibilities of picking the $k_{j}$th feature and corresponding threshold,
$c_{j}$, given the $T$ structure. Equation 9 is used as the prior for the $T$.
This prior is recommended by [5] and three parameters specify this prior: the
depth of the $T$, $d(T)$; the parameter, $a$, which acts as a normalisation
constant; the parameter, ${\beta>0}$, which specifies how many leaf nodes are
probable, with larger values of $\beta$ reducing the expected number of leaf
nodes. $\beta$ is crucial as this is the penalizing feature of our
probabilistic $T$ which prevents an algorithm that uses this prior from over-
fitting and allows convergence to occur[20]. Changing $\beta$ allows us to
change the prior probability associated with “bushy” trees, those whose leaf
nodes do not vary too much in depth.

An exhaustive evaluation of equation 4 over all trees will not be feasible,
except in trivially small problems, because of the sheer number of possible
trees.

Despite these limitations, Bayesian algorithms can still be used to explore
the posterior. Such algorithms simulate a chain sequence of trees, such as:

$T_{0},T_{1},T_{2},....,T_{n}$ (10)

which converge in distribution to the posterior, which is itself proportional
to the joint distribution, $p(Y_{1:N}|T,\theta(T),x_{1:N})p(\theta(T)|T)p(T)$,
specified in equation 4. We choose to have a simulated sequence that
gravitates toward regions of the higher posterior probability. Such a
simulation can be used to search for high-posterior probability trees
stochastically. We now describe the details of algorithms that achieve this
and how they can be implemented.

###  2.1 Stochastic Processes on Trees

To design algorithms that can search the space of trees stochastically, we
first need to define a stochastic process for moving between trees. More
specifically, we consider the following four kinds of move from one $T$ to
another:

  * •

Grow(G) : we sample one of the leaf nodes, $j\in L(T)$, and replace it with a
new node with parameters, $k_{j}$ and a $c_{j}$, which we sample uniformly
from their parameter ranges.

  * •

Prune(P) : we sample the $j$th node (uniformly) and make it a leaf.

  * •

Change(C) : we sample the $j$th node (uniformly) from the non-leaf nodes,
$\bar{L}(T)$, and sample $k_{j}$ and a $c_{j}$ uniformly from their parameter
ranges.

  * •

Swap(S) : we sample the $j_{1}$th and $j_{2}$th nodes uniformly, where
$j_{1}\neq j_{2}$, and swap $k_{j_{1}}$ with $k_{j_{2}}$ and $c_{j_{1}}$ with
$c_{j_{2}}$.

We note that there will be situations (eg when moving from a $T$ with a single
node) when some moves cannot occur. We assume each ‘valid’ move is equally
likely and this then makes it possible to compute the probability of
transition from one $T$, to another, $T^{\prime}$, which we denote
$q\left(T^{\prime},\theta(T^{\prime})|T,\theta(T)\right)$.

##  3 Our Approach on Evolutionary Algorithms

Evolutionary Algorithms (EA) mimic living organisms’ behavior, using natural
mechanisms to solve problems[2]. In our approach, the optimisation problem is
represented as a multi-dimensional space on which the population lives (in our
case, the population is the total number of trees). Each location on the
solution space corresponds to a feasible solution to the optimisation problem.
The optimal solution is found by identifying the optimal location in the
solution space.

Pheromones play a crucial role in evolutionary algorithms as they are used for
communication among the population [7, 15]. When a member of the population
moves, it receives pheromones from the other member of the population and uses
them to determine its next move. Once all members of the population reach new
locations, they release pheromones; the concentration and type of pheromones
released depend on the objective function or fitness value at that location.
The solution space is the medium for transmitting pheromones, allowing
individuals to receive and be affected by the pheromones released by other
individuals, creating a global information-sharing model.

The population will gradually gain a rough understanding of global information
through their movements, which can significantly benefit the optimisation
process. In our approach, the EA can use the solution space as a memory to
record the best and worst solutions produced in each iteration. Once the
positioning stage is finished, all pheromones on the solution space are
cleared. To guide the optimisation process, the most representative extreme
solutions are selected from the recorded solutions, and the corresponding
locations are updated with permanent pheromones. Unlike permanent pheromones,
temporary pheromones only affect the movements of trees in the next iteration.

##  4 Methods

###  4.1 Conventional MCMC

One approach is to use a conventional application of Markov Chain Monte Carlo
to DTs, as found in [11].

More specifically, we begin with a tree, $T_{0}$ and then at the $i$th
iteration, we propose a new $T^{\prime}$ by sampling $T^{\prime}\sim
q\left(T^{\prime},\theta(T^{\prime})|T_{i},\theta(T_{i})\right)$. We then
accept the proposed $T^{\prime}$ by drawing $u\sim U([0,1])$ such that:

$\displaystyle
T_{i+1}=\left\\{\begin{array}[]{ll}T^{\prime}&u\leq\alpha(T^{\prime}|T)\\\
T_{i}&u>\alpha(T^{\prime}|T)\end{array}\right.$ (13)

where we define the acceptance ratio, $\alpha(T^{\prime},T)$ as:

$\displaystyle\alpha(T^{\prime}|T)=\frac{p(Y_{1:N}|T,\theta(T),x_{1:N})}{p(Y_{1:N}|T^{\prime},\theta(T^{\prime}),x_{1:N})}\frac{q\left(T,\theta(T)|T^{\prime},\theta(T^{\prime})\right)}{q\left(T^{\prime},\theta(T^{\prime})|T,\theta(T)\right)}$
(14)

.

This process proceeds for $n$ iterations. We take the opportunity to highlight
that this process is inherently sequential in its nature.

###  4.2 Evolutionary Algorithm in Bayesian Decision Trees

####  4.2.1 Initializing population

The initial population plays a crucial role in the solutions’ quality. An
initial population with a good mix of diversity and a substantial number of
trees can help improve the optimisation performance of the algorithm. To
create a diverse initial population, a random method is often employed. This
approach helps to ensure that the algorithm can perform a global search
efficiently, is easy to use, and has a diverse set of individuals in the
initial population.

The population size of trees is invariable, denoted as $n$. The location of
every $T_{i}$ in the $D$ \- dimensional space can be described as
$:T=(T_{1}^{1},T_{2}^{2},…,T_{n}^{d},…,T_{N}^{D})$. According to the value of
our objective function $p(Y1:N,T,\theta(T)|x1:N)$ , the fitness value of
location $d$ and $T_{i}$ can be measured. By comparing the current location of
each $T_{i}$ in the initial population, the optimal location and the worst
location in the initial population were obtained, and the value of the
objective function of the optimal location in the initial population was
recorded.

####  4.2.2 Positioning stage

In the positioning stage, in our use case, trees release permanent and
temporary pheromones. The solution space records the locations of the terrible
solution and the excellent solution produced by each iteration. While all
trees move to the new location, the pheromones are updated differently, which
will be discussed in section 4.2.4. The process in the positioning stage is
shown in Figure 2 In our case, the possible movements of the $T_{i}$ are those
described in section 2.1. When the proposed move is the $Grow$, we search for
possible solutions in a higher dimensional space; when the proposed move is
the $Prune$, we explore for possible solutions in lower dimensional space; and
when the proposed moves are $Change$ and $Swap$, we search for solutions on
the same dimensional space.

Figure 2: Positioning Stage

####  4.2.3 Permanent Pheromones

$Permanent\;pheromones$ have persistent effects. On each iteration, we
evaluate each $T_{i}$ to compute the value $\alpha_{n}$ (see equation 14). If
$\alpha_{n}$ is greater than a uniform number between $[0,1]$, we store the
$T_{i}$ and the stochastic move associated with $positive\;exploration$ and
$effective\;moves$, respectively. If $\alpha_{n}$ is less than a uniform
number between $[0,1]$ we store the $T_{i}$ and the stochastic move
associated, on $negative\;explartion$ and $innefective\;moves$
respectively($negative\;explartion$ and $innefective\;moves$ will be used on
$temporary\;pheromones$). We repeat the procedure above for all trees. We then
calculate the $permanent\;pheromones$ given the $effective\;moves$.
$Permanent\;pheromones$ is a single list with $4$ real numbers, representing
the possibilities of each one of the $4$ stochastic moves to be selected on
the next iteration. We update the $permanent\;pheromones$ by adding to each
possibility on the $permanent\;pheromones$, the number of times each move is
in the $effective\;moves$ list, and we then divide each element of the list by
the sum of the $effective\;moves$ (we normalise to add up to 1). We have also
designed a mechanism to avoid having a biased $permanent\;pheromones$ list.
For example, in the first stages of the algorithm, it is common the $Grow$
move to be a more useful proposal compared to the $Prune$. In such cases, if a
specific move has more than $80\%$ possibilities to be chosen, we re-weight
the list by setting the dominant move having $40\%$ possibilities to be
selected and the rest having $20\%$ possibilities to be chosen on the next
iteration.

####  4.2.4 Temporary Pheromones

$Temporary\;pheromones$ can only affect the movements of trees in the next
iteration. We discussed above how we end up with a list called
$effective\;moves$ and a list called $ineffective\;moves$. When the iteration
ends,$temporary\;pheromones$ clear, and each $T_{i}$ will release new
$temporary\;pheromones$, which will be recorded. Algorithm 1 shows how we
produce and store pheromones.

###  4.3 Sequential Monte Carlo with EA

We are considering a Sequential Monte Carlo(SMC) sampler[6] to handle the
problem of sampling the trees from the posterior. After we have collected all
the useful information from the movement of the trees, we now need to sample
using the pheromones produced. As the algorithm 2 shows, there are three
possible sampling techniques where one has a subcategory. We choose the
sampling technique by drawing a uniform number between $[0,1]$. At this point,
we need to specify that on each iteration, we draw a new uniform number, so
each $T_{i}$ has the possibility to sample with a different strategy.

The first sampling technique uses the $temporary\;pheromones$, and it has a
possibility of $45\%$ to be chosen. As mentioned earlier,
$temporary\;pheromones$ are produced during the previous iteration. This
sampling technique has a subcategory, where the samples classified in
$positive\;exploration$ use a different sampling technique from those listed
on $negative\;exploration$. If the $T_{i}$ is listed in
$positive\;exploration$, we pick a stochastic move $m$ from the list with the
$ineffective\;moves$ uniformly. We then remove all the identical $m$ from the
$ineffective\;moves$ and sample uniformly with equal probabilities from the
remaining $ineffective\;moves$. If the $T_{i}$ is in $negative\;exploration$,
we pick uniformly a stochastic move $m$ from the $effective\;moves$ list to
sample the particular $T_{i}$.

The second sampling technique uses the $permanent\;pheromones$, and it has a
possibility of $45\%$ to be chosen. As mentioned earlier,
$permanent\;pheromones$ are updated dynamically on every iteration,
considering all the previous iterations. We sample the $T_{i}$ using the
possibilities in the list we describe in subsection A$permanent\;pheromones$
in section 4.2.

The last sampling technique is straightforward. We use a list with stochastic
moves, where each move has a uniform probability of being chosen. This
sampling technique is unaffected by the pheromones. The main reason for
including this technique is to ensure our algorithm is not biased, as this is
the most common way of sampling in Bayesian DTs.

Algorithm 1 Pheromones Production Stage

Initialise $n$ trees($T$)

Sample trees $[T_{0},T_{1},...,T_{n}]$

Initialise initial possibilities =
$[p(G),p(P),p(C),p(S)]=[0.25,0.25,0.25,0.25]$

Initialise permanent pheromones = $[p(G),p(P),p(C),p(S)]$ $\triangleright$
permanent pheromones

Initialise $positive\;exploration$ list

Initialise $effective\;moves$ list $\triangleright$ Temporary Pheromones

Initialise $negative\;exploration$ list

Initialise $ineffective\;moves$ list $\triangleright$ Permanent Pheromones

iterations = $10$

for $(i\leq iterations,i++)$ do

Evaluate trees $[T_{0},T_{1},...,T_{n}]$

Store their acceptance probability $[\alpha_{0},\alpha_{1},...,\alpha_{n}]$

for $(s\leq n,s++)$ do

Draw a uniform number $u_{1}\sim u[0,1]$

if $\alpha_{s}>u_{1}$ then

append $T_{s}$ to list $positive\;exploration$

append $T_{s}move$ to list $effective\;moves$

else

append $T_{s}$ to list $negative\;expolation$

append $T_{s}move$ to list $ineffective\;moves$

end if

end for

update $permanent\;pheromones$ list given the $effective\;moves$ list

end for

Algorithm 2 SMC with EA

for $i\leq n,i++$ do

Draw a uniform number $u_{2}\sim u[0,1]$

if $u_{2}\leq 0.45$ then

if $T_{i}\;in\;positive\;exploration$ then

pick uniformly a move $m$ from $ineffective\;moves$

Remove every identical $m$ from the $ineffective\;moves$

Sample $T_{i}$ uniformly given the updated $ineffective\;moves$

end if

if $T_{i}\;in\;negative\;exploration$ then

Pick uniformly a move $m$ from $effective\;moves$

Sample $T_{i}$ by applying the selected move $m$

end if

else if $u_{2}>0.45$ and $u_{2}\leq 0.9$ then

Sample $T_{i}$ using $permanent\;pheromones$

else

Sample $T_{i}$ using $initial\;possibilities$

end if

end for

Empty $positive\;exploration$

Empty $effective\;moves$ $\triangleright$ Temporary Pheromones

Empty $negative\;exploration$

Empty $ineffective\;moves$ $\triangleright$ Permanent Pheromones

##  5 Experimental Setup and Results

To demonstrate the accuracy improvement we achieved through our proposed
methods, we experiment on three publicly available
datasets111https://archive.ics.uci.edu/ml/index.php listed in Table 1. We
acknowledge that the size of the datasets we conduct the experiments is small.
The main reason is to have a fair comparison between the MCMC algorithm and
SMC, as experiments show[12] that the former struggles to converge on an
adequate time on big datasets, compared to the latter. We also aim to show
that SMC, an inherently parallel algorithm combined with EA, can be a great
fit within the context of big data. For each dataset, we have multiple testing
hypotheses. Firstly, we compare the SMC-EA with MCMC on $1000$ iterations and
$10$ chains for MCMC and $10$ $T$ for the SMC-EA, $100$ iterations and $100$
chains for MCMC and $100$ $T$ for the SMC-EA, and $10$ iterations with $1000$
chains for MCMC and $1000$ $T$ for the SMC-EA.

Table 1: Datasets description Dataset | Attributes | Instances
---|---|---
Heart Disease | 75 | 303
Lung Cancer | 56 | 32
SCADI | 206 | 70

This section presents the experimental results obtained using the proposed
methods with a focus on accuracy improvement and the ability of the SMC-EA to
evolve smoothly in a very short period of iterations. We obtained the
following results using a local HPC platform comprising twin SP 6138
processors (20 cores, each running at 2GHz) with 384GB memory RAM. We use the
same hyper-parameters $\alpha$ and $\beta$ for every contesting algorithm for
testing purposes and a fair comparison and evaluation.

We tested both MCMC and SMC-EA with a 5-Fold Cross-Validation. Results
indicate what is discussed in Section 4.2.1. When we initialise more trees, in
our case trees, we introduce a more diverse set, which helps improve the
algorithm’s optimization performance. More specifically, SMC-EA with $1000$
$T$ and $10$ iterations running SCADI dataset has an accuracy improvement of
$\sim 2\%$ and $\sim 6\%$ compared to having $100$ $T$ with $100$ iterations,
and $10$ $T$ with $1000$ iterations respectively. On the Heart Disease dataset
SMC-EA with $1000$ $T$ and $10$ iterations has an accuracy improvement of
$\sim 3\%$ and $\sim 4\%$ compared to having $100$ $T$ with 100 iterations and
$10$ $T$ with $1000$ iterations respectively. On the Lung Cancer dataset SMC-
EA with $1000$ $T$ and $10$ iterations, has an accuracy improvement of $\sim
8\%$ and $\sim 14\%$ compared to having $100$ $T$ with $100$ iterations, and
$10$ $T$ with $1000$ iterations respectively.

On the other hand, MCMC performs poorly when we have fewer iterations and more
chains compared to SMC-EA. This is expected as MCMC needs adequate time to
converge. When MCMC ran for more iterations, the algorithm made better
predictions, and the accuracy achieved cannot be accepted as an acceptable
threshold. SMC-EA has a $\sim 12\%$ better predictive accuracy on the SCADI
dataset compared to MCMC, $\sim 7\%$ on Heart Disease, and $\sim 17\%$ on Lung
Cancer(see Tables 2, 3 and 4).

On the SMC-EA algorithm, trees do not have a leading $T_{i}$, and their
movements are guided by interactions between them rather than one individual
$T$ dominating the group of trees. This helps to prevent individualism and
stagnation in the population’s evolution. Furthermore, a diverse population of
trees is more effectively handled by the approach we suggest, as we avoid
falling into local optimisation. The stochastic nature of the SMC-EA algorithm
helps in escaping local optimisation and achieving global optimisation.
Combining positive and negative feedback from the different pheromones can
incorporate the benefits of successful solutions while mitigating the negative
effects of poor solutions. SMC-EA algorithm fully uses all the information on
the solution space, avoiding unnecessary waste or duplication of information.
The EA algorithm updates the positions of all trees by using a combination of
their current and past positions within the population. This allows the
algorithm to maintain a history of information, preventing rapid jumps and
leading to a smooth algorithm evolution.

Table 2: SCADI dataset Chains_Trees | Iterations | MCMC | SMC-EA
---|---|---|---
10 | 1000 | 85 | 90.2
100 | 100 | 63 | 94.2
1000 | 10 | 57 | 96.6
Table 3: Heart Disease dataset Chains_Trees | Iterations | MCMC | SMC-EA
---|---|---|---
10 | 1000 | 75.7 | 78.9
100 | 100 | 73.5 | 79.2
1000 | 10 | 67.7 | 82.7
Table 4: Lung Cancer dataset Chains_Trees | Iterations | MCMC | SMC-EA
---|---|---|---
10 | 1000 | 70.1 | 73.9
100 | 100 | 69.7 | 79.1
1000 | 10 | 69.1 | 87.2

Due to the small size of the datasets we are using to conduct this study, we
can only show the effectiveness of our method in exploring the solutions space
faster. However, previous studies[12] have shown that the SMC DT algorithm can
improve the runtime compared to MCMC DT by up to a factor of $343$. We are
optimistic that we can achieve the same results, as SMC and EA are inherently
parallel algorithms. The main bottleneck of the SMC-EA algorithm is when we
evaluate a big number of trees; for example, see the test case of $10$
iterations and $1000$ $T$. We can overcome this problem by the distributed
implementation, as we can distribute the trees on many nodes and evaluate the
trees concurrently.

##  6 Conclusion

Our study has shown that by combining two novel algorithms, the SMC and EA can
tackle major problems on Bayesian Decision Tress and open the space for more
research. According to our experimental results, our novel approach based on
Sequential Monte Carlo and Evolutionary Algorithms explores the solution space
with at least 100 times fewer iterations than a standard probabilistic
algorithm, for example, Markov Chain Monte Carlo. As discussed in section 5,
we managed to tackle the problem of the naive proposals, and we suggest a
method that takes advantage of the communication between the trees. We
proposed a sophisticated method to propose new samples based on EA and
minimised the burn-in period through SMC.

This study, though, has some limitations. First, as we mentioned, both SMC and
EA are inherently parallel algorithms with existing parallel
implementations[26, 21, 25]. We plan to extend this study by parallelising the
SMC-EA, adding more testing scenarios with larger datasets, and showing
improvement in run time.

## References

  * [1] Guillaume Basse, Aaron Smith, and Natesh Pillai.  Parallel markov chain monte carlo via spectral clustering.  In Arthur Gretton and Christian C. Robert, editors, Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, volume 51 of Proceedings of Machine Learning Research, pages 1318–1327, Cadiz, Spain, 09–11 May 2016. PMLR.
  * [2] S Binitha, S Siva Sathya, et al.  A survey of bio inspired optimization algorithms.  International journal of soft computing and engineering, 2(2):137–151, 2012.
  * [3] Sergio Blanes, Fernando Casas, and Jesús María Sanz-Serna.  Numerical integrators for the hybrid monte carlo method.  SIAM Journal on Scientific Computing, 36(4):A1556–A1580, 2014.
  * [4] Hugh A Chipman, Edward I George, and Robert E McCulloch.  Bayesian cart model search.  Journal of the American Statistical Association, 93(443):935–948, 1998.
  * [5] Hugh A Chipman, Edward I George, and Robert E McCulloch.  Bart: Bayesian additive regression trees.  The Annals of Applied Statistics, 2010.
  * [6] Pierre Del Moral, Arnaud Doucet, and Ajay Jasra.  Sequential monte carlo samplers.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(3):411–436, 2006.
  * [7] Marco Dorigo and Thomas Stützle.  Ant colony optimization: overview and recent advances.  Springer, 2019.
  * [8] Randal Douc, Arnaud Guillin, J-M Marin, and Christian P Robert.  Convergence of adaptive mixtures of importance sampling schemes.  The Annals of Statistics, 35(1):420–448, 2007.
  * [9] Efthyvoulos Drousiotis, Panagiotis Pentaliotis, Lei Shi, and Alexandra I Cristea.  Capturing fairness and uncertainty in student dropout prediction–a comparison study.  In Artificial Intelligence in Education: 22nd International Conference, AIED 2021, Utrecht, The Netherlands, June 14–18, 2021, Proceedings, Part II, pages 139–144. Springer, 2021.
  * [10] Efthyvoulos Drousiotis, Lei Shi, and Simon Maskell.  Early predictor for student success based on behavioural and demographical indicators.  In International Conference on Intelligent Tutoring Systems. Springer, 2021.
  * [11] Efthyvoulos Drousiotis and Paul G Spirakis.  Single mcmc chain parallelisation on decision trees.  In Learning and Intelligent Optimization: 16th International Conference, LION 16, Milos Island, Greece, June 5–10, 2022, Revised Selected Papers, pages 191–204. Springer, 2023.
  * [12] Efthyvoulos Drousiotis, Paul G Spirakis, and Simon Maskell.  Parallel approaches to accelerate bayesian decision trees.  arXiv preprint arXiv:2301.09090, 2023.
  * [13] Simon Duane, Anthony D Kennedy, Brian J Pendleton, and Duncan Roweth.  Hybrid monte carlo.  Physics letters B, 195(2):216–222, 1987.
  * [14] Joffrey Dumont Le Brazidec, Marc Bocquet, Olivier Saunier, and Yelva Roustan.  Quantification of uncertainties in the assessment of an atmospheric release source applied to the autumn 2017.  Atmospheric Chemistry and Physics, 2021.
  * [15] Vijay Kalivarapu, Jung-Leng Foo, and Eliot Winer.  Improving solution characteristics of particle swarm optimization using digital pheromones.  Structural and Multidisciplinary Optimization, 37(4):415–427, 2009\.
  * [16] Jun S Liu, Faming Liang, and Wing Hung Wong.  The multiple-try method and local optimization in metropolis sampling.  Journal of the American Statistical Association, 95(449):121–134, 2000.
  * [17] Enzo Marinari and Giorgio Parisi.  Simulated tempering: a new monte carlo scheme.  EPL (Europhysics Letters), 19(6):451, 1992.
  * [18] Per Mykland, Luke Tierney, and Bin Yu.  Regeneration in markov chain samplers.  Journal of the American Statistical Association, 90(429):233–241, 1995.
  * [19] Christian P Robert, Víctor Elvira, Nick Tawn, and Changye Wu.  Accelerating mcmc algorithms.  Wiley Interdisciplinary Reviews: Computational Statistics, 10(5):e1435, 2018.
  * [20] Veronika Ročková and Enakshi Saha.  On theory for bart.  In The 22nd International Conference on Artificial Intelligence and Statistics. PMLR, 2019.
  * [21] Urvashi Prakash Shukla and Satyasai Jagannath Nanda.  Parallel social spider clustering algorithm for high dimensional datasets.  Engineering Applications of Artificial Intelligence, 56:75–90, 2016\.
  * [22] Matthew A Taddy, Robert B Gramacy, and Nicholas G Polson.  Dynamic trees for learning and design.  Journal of the American Statistical Association, 106(493):109–123, 2011.
  * [23] Duncan Taylor, Jo-Anne Bright, and John Buckleton.  Interpreting forensic dna profiling evidence without specifying the number of contributors.  Forensic Science International: Genetics, 2014.
  * [24] Gloria I Valderrama-Bahamóndez and Holger Fröhlich.  Mcmc techniques for parameter estimation of ode based models in systems biology.  Frontiers in Applied Mathematics and Statistics, 2019.
  * [25] Alessandro Varsi, Simon Maskell, and Paul G Spirakis.  An o (log2n) fully-balanced resampling algorithm for particle filters on distributed memory architectures.  Algorithms, 14(12):342, 2021.
  * [26] Alessandro Varsi, Jack Taylor, Lykourgos Kekempanos, Edward Pyzer Knapp, and Simon Maskell.  A fast parallel particle filter for shared memory systems.  IEEE Signal Processing Letters, 27:1570–1574, 2020.
  * [27] Yuhong Wu, Håkon Tjelmeland, and Mike West.  Bayesian cart: Prior specification and posterior simulation.  Journal of Computational and Graphical Statistics, 16(1):44–66, 2007.