\section{Introduction}
Granular flows are ubiquitous in nature and in a large array of industrial processes.
Pharmaceutical powders, plastic granulate and rocks obtained by mining are just some examples of granular media that are used in industries and which are processed in a multitude of different flow states.
While attempts have been made to formulate governing equations for granular flows, like the Navier-Stokes equations for fluid flow, they have so far eluded a general framework \citep{faccanoni2013}.
More recently, research focus has been expanded
towards applications of Deep Learning in the simulation of
physical domains, such as fluid dynamics, deformable materials, or aerodynamics~\citep{SanchezGonzalez2020, Pfaff2020}.
One key component of the recent progress in deep neural network simulation is the usage of
graph neural networks~\citep[][GNNs]{Scarselli2008}.
GNNs have demonstrated effectiveness in settings that involve interactions between many entities via forward dynamics~\citep{Battaglia2018}.
In~\citet{SanchezGonzalez2020}, GNNs
are used to obtain models that generalize fluid dynamics and 2D interactions over many timesteps and different initial conditions.
In~\citet{Pfaff2020}, mesh-based simulations
are learned using GNNs to predict the dynamics of a wide range of physical
systems, including aerodynamics, or structural mechanics.

Compared to previous work \citep[i.e., ][]{SanchezGonzalez2020, Pfaff2020}, we focus on learning 3D simulations of granular particle flow with nontrivial geometric boundary conditions. These simulations are highly relevant for the design of industrial processes and allow to understand and improve particle flow dynamics of various given materials.
Since no underlying governing equation for general granular flows exists, we use simulation data originating from the Discrete Element Method \citep[][DEM]{Cundall1979} as the ground truth (see Appendix \ref{DEM} for further details on DEM).
We generate granular flow simulation data with the open-source DEM software LIGGGHTS \citep{Kloss2012}. LIGGGHTS allows the simulation of particulate flows with a wide range of materials and complex mesh-based wall geometries, and therefore enables the simulation of relevant industrial processes. Such complex mesh-based wall geometries are in contrast to \citet{SanchezGonzalez2020}, where for the 3D simulations static cuboids have been assumed and a material point method  \citep[][MPM]{sulsky1995application} based simulator \citep{Hu2018} has been applied.

Similar to~\citet{SanchezGonzalez2020}, we learn a time-transition model, consisting of an encoder, processor, and, a decoder,  to predict particle accelerations and build upon their relative positional encoding variant.
The time transition model from time $t_k$ to time $t_{k+1}$
is given by
\begin{align}
\dot{\Bp}^{t_{k+1}}=\dot{\Bp}^{t_{k}}+\Delta t\; \ddot{\Bp}^{t_{k+1}} \ , \label{ttmodel} \\
\Bp^{t_{k+1}}=\Bp^{t_{k}}+\Delta t\; \dot{\Bp}^{t_{k+1}} \ ,
\end{align}
where $\Bp$ is the particle location, and $\dot{\Bp}$ the
particle velocity. The new particle acceleration $\ddot{\Bp}$ at time $t_{k+1}$ needs to be predicted.
The contributions of this work towards an accurate modeling of complex 3D granular flows are three-fold:
\begin{itemize}
    \item Our model deals with triangularization of 3D objects, which is a very common representation in engineering applications. We achieve this by inserting virtual nodes if particles are close to boundaries and computing the proximity to the closest triangle.
    \item We include particle - boundary interactions by forcing the network to be invariant with respect to the orientation of the normal vectors.
    \item We compare and analyze relevant physical quantities between simulated processes and the output of our network.
\end{itemize}

\section{Accurate Modeling of Granular Flow Dynamics}\label{sec:GranuFlow}

\textbf{Learning simulations that are governed by complex geometries. }
Triangularization of geometric boundaries is very common in engineering applications.
One idea to tackle this challenge for 2D scenes is the  insertion of stationary, virtual particles into the scene to describe boundaries~\citep{SanchezGonzalez2020}. However, for 3D scenes and thus for many practically relevant engineering applications,
a computationally more efficient approach is required for the following reasons:
Firstly, triangle surface areas in 3D would require significantly more stationary particles for the representation than curves in 2D scenes. Secondly, for certain time frames, only some parts of the mesh are relevant, e.g. as long as particles are in the state of a free fall in a container and the bottom of a container is still far away, the mesh part describing the bottom of the container is irrelevant for the next time step.

\newcommand\BNabla{\bm{\nabla}}
\textbf{Modeling distances to triangular boundaries in 3D scenes. }
 The challenge is to correctly incorporate the modeling of the proximity of a particle to a triangular 3D boundary and at the same time avoid large computational cost.
Firstly, to avoid large computational cost, we dynamically insert boundary particles if a real particle is close to the corresponding boundary.
The inserted virtual nodes in the graph describe whole surface areas in contrast to single 3D points. Therefore, we indicate the particle type, virtual or real, by introducing additional node features, such that the neural network is able to distinguish them. In other terms, this allows the model to learn different dynamics for particle - particle and particle - boundary interactions.
The additional node features are: (i) type feature, i.e.,  a binary indicator of whether a node represents a particle that is real or virtual, and, in the latter case, (ii) the components of the normal vector of the triangular surface (null vectors for real particles).
The second challenge is to correctly model the proximity of a
particle to the boundary, which
is computed as the minimum distance between the particle center and the closest point of the mesh triangles (adopted from \citet{Eberly1999}).
We assume a triangle to be represented by a function $\BT$, which is parameterized by two scalar values $u$ and $v \in \mathbb{R}$:
\begin{equation}
    \BT\left(u,v\right)=\BB+u\;\BE_0+v\;\BE_1,
\end{equation}
where $u\geq0$, $v\geq0$, $u+v\leq1$, $\BB$ represents one of the corner points of the triangle, and, $\BE_0$ and $\BE_1$ are vectors from $\BB$ towards the other two corner points (see Figure~\ref{myTriangle}). We solve the optimization problem
\begin{align}
& \min_{u,v}\; \BQ\left(u,v\right)=\left|\BT\left(u,v\right)-\BP\right|^2 \\
s.t.\;\; & u\geq0,\;v\geq0,\;u+v\leq1 \nonumber
\end{align}
to obtain the minimizing parameters $u,v$ for retrieving the closest point $\BT\left(u,v\right)$ of the triangle from the point $\BP$ in an Euclidean sense. As indicated in Figure~\ref{myTriangle}, we have to distinguish seven cases: the simplest one is \textit{Case 0}, where $\BNabla \BQ\left(u,v\right)=\BZero$ is fulfilled in the inner of the triangle (i.e., $u>0$, $v>0$, $u+v<1$). If $\BNabla \BQ\left(u,v\right)=\BZero$ is fulfilled outside from that triangle, cases where the minimizing point from the triangle is at the triangle border, i.e., the lines $u=0, v=0, u+v=1$ have to be considered. In order to find the minimizing parameters in this case we have to look for stationary points of the one-parameter functions $\BQ\left(u,1-u\right)$, $\BQ\left(u,0\right)$ and  $\BQ\left(0,u\right)$.

\begin{figure}
\begin{center}
\includegraphics[trim=0 10 0 0, scale=0.4]{figures/trianglePlot.pdf}
\end{center}
\caption{Visualization of point - triangle distance calculations in 3D. The triangle is represented by a parameterized function $\BT\left(u,v\right)=\BB+u\;\BE_0+v\;\BE_1$ with $u\geq0$, $v\geq0$, $u+v\leq1$ (indicated by the orange area). Level sets of $\BQ\left(u,v\right)$ are indicated by ellipses and describe the squared Euclidean distance of a triangle point $\BT\left(u,v\right)$ to the point $\BP$, for which we compute the minimum distance.}
\label{myTriangle}
\end{figure}

\textbf{Learning particle - boundary interactions. } In granular flow simulations, usually particle - particle interaction data outweigh particle - boundary interaction data, which makes learning particle - boundary interactions difficult.
We therefore
put more emphasis on the type and normal vector features by introducing hyperparameters for the features, that allow discrimination between real particles and virtual boundary particles.
The GNN predictions should be independent of the orientation of the normal vectors representing the planes of the triangle walls. Since there is a positive and a negative choice, we use both a positively and a negatively oriented version of the respective normal vectors as input features. However, as the prediction of a network should be invariant with respect to the orientation of the normal vectors, we define a partial ordering to be able to sort the normal vectors with respect to their orientations. For a given normal vector $\Bn=\PAR{n_1,n_2,n_3} \in \mathbb{R}^3$, we use the following partial order function
\begin{equation}
f_o\PAR{\Bn}=\sum_{i=1}^{3} 3^{i-1}\;\PAR{\sgn\PAR{n_i}+1}
\end{equation}
to retrieve the scalar values $o_1=f_o\PAR{\Bn}$ and $o_2 = f_o\PAR{-\Bn}$ and sort the two vectors according to their corresponding mapped values.
We choose a base of $3$ since zero entries of vectors are  possible. Different sign combinations of the normal vectors are shown in Figure~\ref{signBalls}.

\begin{figure}[ht]
\begin{center}
\includegraphics[trim=130 0 0 0, scale=0.27]{ballPlots/plot0.pdf}
\includegraphics[trim=15 0 15 0, scale=0.27]{ballPlots/plot1.pdf}
\includegraphics[scale=0.27, trim=0 0 130 0]{ballPlots/plot2.pdf}
\end{center}
\caption{Partial ordering of normal vectors. The numbers indicate the signs of normal vector components, which are used in the partial order function. The figures from left to right visualize different (ordered) sign combinations. Sets of sign combinations without zero values form volumes (left), sets with one zero value form planar areas (middle), and sets with two zero values form line sections (right). The 0/0/0 combination forms a point at the origin. }
\label{signBalls}
\end{figure}

To test the performance of our approach, we conduct a toy experiment as well as a simulation experiment with different representations of normal vectors. We describe both experiments in the appendix (Section \ref{toyExp} and Section \ref{simExp}).

\section{Results}

We focus on two typical applications of granular flow simulations in the design of industrial machinery: the particle flow through a hopper and the particle dynamics in a rotating drum.
Figure \ref{paravisplot} visualizes
the particle positions for both problem settings at different
time steps in the ground truth simulation and in the prediction of our model.
For both applications gravitation acts along the z-direction.
The upper part of the hopper is delimited along the y-axis by two planes, which are parallel to the x-z plane. The x-axis is delimited by two planes, that are inclined at certain angles $\alpha,\, 180^{\circ}-\alpha$ to the x-y plane and at corresponding angles $\alpha-90^{\circ}, 90^{\circ}-\alpha$ to the y-z plane. At the bottom of the hopper there is a hole with an adjustable radius, that is initially closed. Our generated training data consists of 30 simulation trajectories with different angles $\alpha$ and different hole sizes. Moreover, the initial filling distribution is varied.
The rotation axis of the drum is the y-axis. The initial filling is obtained by rotating the filling of a resting drum around the x-axis.

In Figure \ref{cmpplot} we compare ground truth trajectories to the learned trajectories.
The machine learning model does not exactly reproduce the ground truth trajectories due to chaotic behavior in the long run, but aggregated quantities like particle-averaged positions (upper left plot) show good qualitative agreement. The upper right and the lower left plots of the figure show that also the time- and particle-averaged particle flows are well reproduced. In the lower right plot we analyse the machine learning model in terms of mixing behavior. We quantify the extend of particle mixing via the mixing entropy \citep{fang1975appl}, which is computed by splitting particles into two classes +1, -1 at a certain time step $t_0=30$ and then summing over weighted local entropies $s(\Bx_{klm}, t)$ located at grid cells $\Bx_{klm}$ (see Appendix \ref{mixent} for details).  The lower right plot shows that the expected increase of the mixing entropy over time is well described by the GNN.
Overall, we conclude that our GNN approach is able to model 3D granular flow simulations accurately.

\section{Outlook}

Although, we have already used comparably large time steps in our machine learning model (i.e. multiples
of the ground truth simulation training data), we have not fully investigated the potential
of machine learning models for granular simulation data yet: Firstly, the simulation time of the ground truth simulation is strongly material-dependent. As a downside, the true Young‘s Modulus, which describes a stress-strain relationship in the linear elastic region of a material, is usually not used in the computation of simulations, but a much smaller value.
Machine learning models therefore could allow the simulation of granular flows with more realistic material parameters.
Secondly, larger time steps would allow to model physical phenomena occurring on larger time scales, which cannot be modeled accurately with current DEM approaches.
Further, we plan to leverage the
multiple symmetries occurring in the geometry of granular flow problems.
Another research direction might be the usage of global parameters
such as the angles of the side walls of the hopper as additional input features to the neural network.

\subsubsection*{Acknowledgments}
This research was supported by FFG grant 871302 (DL for granular flow).

The ELLIS Unit Linz, the LIT AI Lab, the Institute for Machine Learning, are supported by the Federal State Upper Austria. IARAI is supported by Here Technologies. We thank the projects AI-MOTION (LIT-2018-6-YOU-212), DeepToxGen (LIT-2017-3-YOU-003), AI-SNN (LIT-2018-6-YOU-214), DeepFlood (LIT-2019-8-YOU-213), Medical Cognitive Computing Center (MC3), PRIMAL (FFG-873979), S3AI (FFG-872172), ELISE (H2020-ICT-2019-3 ID: 951847), AIDD (MSCA-ITN-2020 ID: 956832).
We thank Janssen Pharmaceutica (MaDeSMart, HBC.2018.2287), Audi.JKU Deep Learning Center, TGW LOGISTICS GROUP GMBH, Silicon Austria Labs (SAL), FILL Gesellschaft mbH, Anyline GmbH, Google, ZF Friedrichshafen AG, Robert Bosch GmbH, UCB Biopharma SRL, Merck Healthcare KGaA, Software Competence Center Hagenberg GmbH, T\"{U}V Austria, and the NVIDIA Corporation.

\clearpage

\vspace*{-2cm}

\imagelabelset{
coarse grid color = red,
fine grid color = gray,
image label font = \sffamily\small,
image label distance = 2mm,
image label back = black,
image label text = white,
coordinate label font = \normalfont,
coordinate label distance = 2mm,
coordinate label back = white,
coordinate label text = black,
annotation font = \normalfont\small,
arrow distance = 1.5mm,
border thickness = 0.6pt,
arrow thickness = 0.4pt,
tip size = 1.2mm,
outer dist = 0.5cm,
}

\newsavebox\imagebox

\begin{figure}[h]
	\centering
	\makebox[ \textwidth ]{
	\centering
	\sbox{\imagebox}{
	\begin{annotationimage}[]{trim = 0mm 0mm 0mm 10mm, clip=True, width=0.9\textwidth}{figures/all1.jpeg}

\draw[coordinate label = {Hopper at (0.25,1.08)}];
\draw[coordinate label = {Rotating Drum at (0.75,1.08)}];
\draw[coordinate label = {Prediction at (0.125,1.03)}];
\draw[coordinate label = {Ground Truth at (0.375,1.03)}];
\draw[coordinate label = {Prediction at (0.625,1.03)}];
\draw[coordinate label = {Ground Truth at (0.875,1.03)}];
    \draw [thick,->] (-0.04, 1)  -- node[anchor=south, rotate=90] {Time} (-0.04, 0);
	\end{annotationimage}
	}
	\usebox{\imagebox}
	}

	\caption{Particle distributions for Hopper and Drum dynamics. Data obtained by the particle simulator LIGGGHTS (Ground Truth) and our trained graph neural network (Prediction) are compared. Particles are indicated by green spheres, triangular wall areas are yellow, the edges of these triangles are indicated by grey lines. The circular arrow indicates the rotation direction of the Drum. }
	\label{paravisplot}
\end{figure}

\begin{figure}[H]
	\centering
	\makebox[ \textwidth ]{
	\includegraphics[trim=0 100 0 0,clip, height=130px]{figuresSIMDL/avgPos_hop2.pdf}
	\includegraphics[trim=0 100 0 0,clip,height=130px]{figuresSIMDL/flowProfileZ_tsStart30.pdf}
	}
	\makebox[ \textwidth ]{
    \includegraphics[trim=0 100 0 0,clip,height=130px]{figuresSIMDL/flow_hop.pdf}
	\includegraphics[trim=0 100 0 0,clip,height=130px]{figuresSIMDL/entropy.pdf}
	}
	\makebox[ \textwidth ]{
    \includegraphics[trim=0 0 0 990,clip,height=12px]{figuresSIMDL/flow_hop.pdf}
    }
	\caption{Position (upper left) and flow profile (lower left) plots for the Hopper, and, flow profile (upper right) and entropy plot (lower right) for the Drum. The plots visualize ground truth (solid line) vs. predictions (dashed line) in dependence of the time step or a coordinate. The flow profiles are average velocities of of particles at a given time step or $z$ coordinate. The mixing entropies are obtained by splitting particles into two partitions according to a threshold on the respective $x$ or $z$ coordinate at time step 30.}
	\label{cmpplot}
\end{figure}

\clearpage


\section{Introduction}
\input{intro.tex}

\section{Preliminaries}
\input{prelim.tex}

\section{Reduction to the $3$-edge-connected case}
\label{Smu1}
\input{reduction.tex}

\section{Simple randomized algorithm}
\input{randomized.tex}

\section{Deterministic algorithm}
\input{deterministic.tex}

\section{Reconstructing the structure of $4$-edge-connected components}
\label{Smu2}
\input{3cuttree.tex}

\section{Open problems}
\label{open-problems}
\input{open-problems.tex}

\bibliographystyle{abbrv}

\subsection{One tree edge}
\label{ssec:deterministic_1edge}

Recall that in this case, we are to find all $3$-edge-cuts intersecting a~fixed depth-first search tree ${\ensuremath{\mathcal{T}}}$ in a~single edge.
This is fairly straightforward: if some $3$-edge-cut $C$ contains exactly one tree edge $e$ of ${\ensuremath{\mathcal{T}}}$, then the border of the cut is exactly ${\ensuremath{\mathcal{T}}}_e$; hence, this cut must include all back edges leaping over $e$.
Therefore, $e$ belongs to a~$3$-edge-cut of this kind if and only if $\low[3](e) = \bot$, i.e. if there only exist two back edges leaping over $e$; in this case, $e$ and these two edges form a $3$-edge-cut.
Since this check can be easily performed in ${\ensuremath{\mathcal{O}}}(n)$ time for all tree edges in ${\ensuremath{\mathcal{T}}}$, the whole subroutine runs in ${\ensuremath{\mathcal{O}}}(n + m)$ time complexity.

\subsection{Two tree edges, lower case}
\label{ssec:deterministic_2edges_down}

Recall that this case requires us to find all $3$-edge-cuts intersecting a~fixed depth-first search tree ${\ensuremath{\mathcal{T}}}$ in two edges, say $e$ and $f$, such that $e < f$ and the remaining back edge $g$ connects ${\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}}_e \setminus {\ensuremath{\mathcal{T}}}_f$ (Figure \ref{fig:randomized-two-tree-edge-cases}).
Hence, $g$ is the only back edge connecting ${\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}}_e \setminus {\ensuremath{\mathcal{T}}}_f$, no back edges connect ${\ensuremath{\mathcal{T}}}_e \setminus {\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$, but there may be multiple back edges connecting ${\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$.

We shall exploit the fact that in a~valid $3$-edge-cut of this kind, there are no back edges originating from ${\ensuremath{\mathcal{T}}}_e \setminus {\ensuremath{\mathcal{T}}}_f$ and leaping over $e$.
For a~fixed edge $e$, this severely limits the set of possible tree edges $f$ below $e$ in ${\ensuremath{\mathcal{T}}}$, since all back edges leaping over $e$ must originate from ${\ensuremath{\mathcal{T}}}_f$.
This is formalized by the notion of the \emph{deepest down cut} of~$e$:

\begin{definition}[deepest down cut]
  \label{def:deterministic_deepest_cut}
  For a~tree edge $e$ in ${\ensuremath{\mathcal{T}}}$, we define the \emph{deepest down cut} of~$e$, denoted $\DeepestDnCut(e)$, as the deepest tree edge $f \geq e$ for which all back edges leaping over $e$ originate from ${\ensuremath{\mathcal{T}}}_f$.
\end{definition}

\begin{lemma}
  \label{lem:deterministic_deepest_cut_def}
  For every tree edge $e$, $\DeepestDnCut(e)$ is defined correctly and uniquely. Moreover, every $3$-edge cut $\{e, f, g\}$ of the considered kind satisfies $e < f \leq \DeepestDnCut(e)$.
  \begin{proof}
    Let ${\ensuremath{\mathcal{S}}}_e$ be the set of the tails of all the edges leaping over~$e$.
    Then, $\DeepestDnCut(e)$ is the deepest tree edge $f=uv$ such that ${\ensuremath{\mathcal{S}}}_e \subseteq V({\ensuremath{\mathcal{T}}}_v)$.
    Equivalently, $v$ must be an ancestor of all the vertices in ${\ensuremath{\mathcal{S}}}_e$; hence, the deepest such vertex, the \emph{lowest common ancestor} of ${\ensuremath{\mathcal{S}}}_e$, is defined uniquely.
    The correctness of the definition of $\DeepestDnCut(e)$ follows.
    For the latter part of the lemma, observe that $f \leq \DeepestDnCut(e)$ if and only if every edge leaping over~$e$ originates from ${\ensuremath{\mathcal{T}}}_f$.
  \end{proof}
\end{lemma}

\begin{lemma}
  \label{lem:deterministic_deepest_cut_lca}
  For every tree edge $e$, $\DeepestDnCut(e) = uv$ is the tree edge whose head $v$ is the lowest common ancestor of two vertices: the tails of $\mindn[1](e)$ and $\maxdn[1](e)$.
  \begin{proof}
  Let ${\ensuremath{\mathcal{S}}}_e$ be the set defined in the proof of Lemma~\ref{lem:deterministic_deepest_cut_def}.
Then, $\DeepestDnCut(e) = uv$ is the tree edge whose head $v$ is the lowest common ancestor of ${\ensuremath{\mathcal{S}}}_e$.
  Equivalently, $v$ is the lowest common ancestor of two vertices: one of minimum preorder in ${\ensuremath{\mathcal{S}}}_e$, and the other of maximum preorder in ${\ensuremath{\mathcal{S}}}_e$.
  By definition, those two vertices are exactly the tails of $\mindn[1](e)$ and $\maxdn[1](e)$, respectively.
  \end{proof}
\end{lemma}

Thanks to Lemma~\ref{lem:deterministic_deepest_cut_lca}, we can compute $\DeepestDnCut(e)$ in constant time for each tree edge $e$ as the lowest common ancestor of two vertices can be computed in constant time after linear preprocessing~\cite{DBLP:journals/siamcomp/HarelT84}.

We now shift our focus to the lower tree edge $f$.
As in Section~\ref{sec:randomized_3cut}, the only possible candidate $g$ for a~back edge of the cut leaping over $f$ is given by $\maxup[1](f)$.
The only problematic part is locating the remaining tree edge $e$.
Previously, we utilized randomness in order to calculate the compressed hash of $e$ given the compressed hashes of $f$ and $g$.
Here, we instead use the following fact:

\begin{lemma}
  \label{lem:deterministic_deepest_edge_only}
  If $\{e, f, g\}$ is a~$3$-edge-cut of the considered kind for some tree edges $e$ and $f > e$, then $e$ is the deepest tree edge satisfying $\DeepestDnCut(e) \geq f$.
  \begin{proof}
    Naturally, the condition $\DeepestDnCut(e) \geq f$ must be satisfied if $\{e, f, g\}$ is a~$3$-edge-cut (Lemma~\ref{lem:deterministic_deepest_cut_def}).
    It only remains to prove that $e$ is the deepest edge with this property.
    Pick two tree edges $e_1$, $e_2$ such that $e_1 < e_2 < f$, both satisfying the condition regarding their deepest down cuts. We then partition the tree into four parts: ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_{e_1}$, ${\ensuremath{\mathcal{T}}}_{e_1} \setminus {\ensuremath{\mathcal{T}}}_{e_2}$, ${\ensuremath{\mathcal{T}}}_{e_2} \setminus {\ensuremath{\mathcal{T}}}_f$, and ${\ensuremath{\mathcal{T}}}_f$. By the assumptions, all back edges leaping over $e_1$ or over $e_2$ must originate from ${\ensuremath{\mathcal{T}}}_f$ (Figure~\ref{fig:deterministic_only_deepest_dn_cut}).

    Since $G$ is $3$-edge-connected, $\{e_2, f\}$ cannot be its edge cut; there must exist a~back edge $x_1$ connecting ${\ensuremath{\mathcal{T}}}_{e_2} \setminus {\ensuremath{\mathcal{T}}}_f$ with the rest of the graph.
    It cannot leap over $e_2$ as it would also need to originate from ${\ensuremath{\mathcal{T}}}_f$; hence, $x_1$ connects ${\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}}_{e_2} \setminus {\ensuremath{\mathcal{T}}}_f$.
    Analogously, $\{e_1, e_2\}$ is not an~edge cut of $G$, from which we infer that there exists another back edge $x_2$ connecting ${\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}}_{e_1} \setminus {\ensuremath{\mathcal{T}}}_{e_2}$.
    Since we produced two separate back edges $x_1$, $x_2$ connecting ${\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}}_{e_1} \setminus {\ensuremath{\mathcal{T}}}_f$, there cannot exist any $3$-edge-cut containing both $e_1$ and $f$.
    The statement of the lemma follows.
  \begin{figure}[h]
    \centering
    \includegraphics[scale=0.4]{figures/deterministic-only-deepest-dn-cut-inkscape.pdf}
    \caption{The setting in Lemma~\ref{lem:deterministic_deepest_edge_only}. Given that $\DeepestDnCut(e_1) \geq f$ and $\DeepestDnCut(e_2) \geq f$, the red back edges are forbidden from occurring in $G$, while the blue back edges can appear in $G$. We prove that both $x_1$ and $x_2$ appear in the graph.}
    \label{fig:deterministic_only_deepest_dn_cut}
  \end{figure}
  \end{proof}
\end{lemma}

Lemmas~\ref{lem:deterministic_deepest_cut_def} and~\ref{lem:deterministic_deepest_edge_only} naturally lead to the following idea: given a~tree edge $e$, the set of edges $f$ satisfying $\DeepestDnCut(e) \geq f$ is a~path $P_e$ connecting the head of $e$ with the head of $\DeepestDnCut(e)$; we assign this path a~weight $w_e$ equal to the depth of $e$ in~${\ensuremath{\mathcal{T}}}$.
We then invoke Theorem~\ref{thm:min_edge_path} with the tree ${\ensuremath{\mathcal{T}}}$, the weighted paths $\{P_e\,\mid\,e \in E({\ensuremath{\mathcal{T}}})\}$, and $k = 1$.
This allows us to find, in ${\ensuremath{\mathcal{O}}}(n)$ time, for each tree edge $f \in E({\ensuremath{\mathcal{T}}})$, the path $P_e$ of the maximum weight containing $f$ as an~edge.
This path naturally corresponds to the tree edge~$e$ from Lemma~\ref{lem:deterministic_deepest_edge_only}.
This way, for every tree edge $f$, we have uniquely identified a~back edge $g$ and a~tree edge $e$ such that the only possible $3$-edge-cut containing $f$ as a~deeper tree edge, if it exists, is $\{e, f, g\}$.

It only remains to verify that $\{e, f, g\}$ is a $3$-edge-cut.
We need to make sure that:
\begin{itemize}
  \item All back edges leaping over $e$ originate from ${\ensuremath{\mathcal{T}}}_f$.
    The equivalent condition $\DeepestDnCut(e) \geq f$ is guaranteed by the algorithm, so we do not need to check it again.
  \item Exactly one edge ($g$) connects ${\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}}_e \setminus {\ensuremath{\mathcal{T}}}_f$.
    Since $g = \maxup[1](f)$, it suffices to verify that the edge $\maxup[2](f)$, which is a~back edge leaping over $f$ whose head is the deepest apart from $g$, also leaps over $e$ (i.e., it does not terminate in ${\ensuremath{\mathcal{T}}}_e$).
\end{itemize}

It can be easily verified that the implementation of the subroutine is deterministic and runs in linear time with respect to the size of $G$.

\subsection{Two tree edges, upper case}
\label{ssec:deterministic_2edges_up}

Recall that in this case, we are required to find all $3$-edge-cuts intersecting ${\ensuremath{\mathcal{T}}}$ in two edges $e$ and $f$, such that $e < f$ and the remaining back edge $g$ connects ${\ensuremath{\mathcal{T}}}_e \setminus {\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$ (Figure \ref{fig:randomized-two-tree-edge-cases}).
Similarly to the randomized case, we use the fact that given a~tree edge $e$, the tail of $g$ has either the smallest preorder (i.e., $g = \mindn[1](e)$) or the largest preorder (i.e., $g = \maxdn[1](e)$) among all the back edges leaping over $e$.
Without loss of generality, assume that $g = \mindn[1](e)$; the latter case is analogous.

In a~similar vein to previous case, observe that if $\{e, f, g\}$ is a~$3$-edge-cut for some $f > e$, then all back edges leaping over $e$ other than $g$ must originate from ${\ensuremath{\mathcal{T}}}_f$.

This leads to the following slight generalization of $\DeepestDnCut$ (Definition~\ref{def:deterministic_deepest_cut}):
\begin{definition}
  \label{def:deterministic_deepest_nomin_cut}
  For a~tree edge $e$ in ${\ensuremath{\mathcal{T}}}$, we define the value $\DeepestDnCutNoMin(e)$ as the deepest tree edge $f \geq e$ for which all back edges leaping over $e$ other than $\mindn[1](e)$ originate from ${\ensuremath{\mathcal{T}}}_f$.
\end{definition}
The process of computation of $\DeepestDnCut(e)$ asserted by Lemma~\ref{lem:deterministic_deepest_cut_lca} can be easily modified to match our needs: recall that $\DeepestDnCut(e)$ is the tree edge whose head was the lowest common ancestor of the origins of $\mindn[1](e)$ and $\maxdn[1](e)$.
Since we exclude $\mindn[1](e)$ from the condition in Definition~\ref{def:deterministic_deepest_nomin_cut}, we use the edge $\mindn[2](e)$ instead: that is, the back edge leaping over $e$, with the next lowest preorder of the tail.
Therefore, we compute $\DeepestDnCutNoMin(e)$ as the tree edge whose head is the lowest common ancestor of the tails of $\mindn[2](e)$ and $\maxdn[1](e)$ (Figure~\ref{fig:deterministic_upcase_f0}).

\begin{figure}[h]
\centering
  \includegraphics[scale=0.4]{figures/deterministic-upcase-dncut-inkscape.pdf}
\caption{The computation of $f_0 = \DeepestDnCutNoMin(e)$ excludes $g = \mindn[1](e)$, and computes the lowest common ancestor of the remaining back edges leaping over $e$.}
\label{fig:deterministic_upcase_f0}
\end{figure}

Now, fix the shallower tree edge $e$ of the cut.
Let $f_0 := \DeepestDnCutNoMin(e)$.
Then, if some tree edge $f > e$ belongs to the $3$-edge-cut $\{e, f, g\}$, then $f \leq f_0$ (otherwise, there would be multiple edges connecting ${\ensuremath{\mathcal{T}}}_e \setminus {\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$).
The natural question is then: is $\{e, f_0, g\}$ a~$3$-edge-cut?
The only possible problem is that there may exist back edges connecting ${\ensuremath{\mathcal{T}}}_{f_0}$ with ${\ensuremath{\mathcal{T}}}_e \setminus {\ensuremath{\mathcal{T}}}_{f_0}$.
Fortunately, if this is the case, then the back edge $\maxup1(f_0)$, dependent only on $f_0$, is one of these back edges.

Hence, let $h_0 := \maxup[1](f_0)$.
If $h_0$ leaps over $e$, we are done, and $\{e, f_0, g\}$ is an~edge cut.
Otherwise, the subtree ${\ensuremath{\mathcal{T}}}_f$ for the sought $3$-edge-cut $\{e, f, g\}$ must contain the head of $h_0$ (or else $h_0$ would connect ${\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}}_e \setminus {\ensuremath{\mathcal{T}}}_f$).
Let then $f_1$, $e \leq f_1 < f_0$ be the deepest tree edge containing the head of $h_0$, that is, the edge whose head coincides with the head of $h_0$.
Then, the edge $f$ of the $3$-edge-cut $\{e, f, g\}$ must satisfy $e < f \leq f_1$.
We can then repeat this procedure: given $f_1$, we compute $h_1 := \maxup[1](f_1)$, and either $h_1$ leaps over $e$ and we are done, or we calculate another edge $f_2 < f_1$ supplying a~better bound on the depth of $f$ (Figure~\ref{fig:deterministic_upcase_iter}).
Since the graph is finite, this process terminates in $k = {\ensuremath{\mathcal{O}}}(n)$ steps, producing the lowest tree edge $f_k \geq e$ such that no back edge connects ${\ensuremath{\mathcal{T}}}_{f_k}$ with ${\ensuremath{\mathcal{T}}}_e \setminus {\ensuremath{\mathcal{T}}}_{f_k}$, and no back edge other than $g$ connects ${\ensuremath{\mathcal{T}}}_e \setminus {\ensuremath{\mathcal{T}}}_{f_k}$ with ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$.
If $f_k = e$, then no $3$-edge-cut $\{e, f, g\}$ exists for $f > e$.
Otherwise, $\{e, f_k, g\}$ is naturally a~correct $3$-edge-cut.
Since two edges of a $3$-edge-cut uniquely identify the third edge (Lemma~\ref{lem:randomized_unique_remaining_edge}), we conclude that this is the only $3$-edge-cut containing $e$ and $g$.

\begin{figure}[h]
\centering
  \includegraphics[scale=0.4]{figures/deterministic-upcase-iteration-inkscape.pdf}
\caption{The main iteration of the algorithm. We set $f_0 = \DeepestDnCutNoMin(e)$, $h_i = \maxup[1](f_i)$, and $f_{i+1}$ is the tree edge with the same head as $h_i$.}
\label{fig:deterministic_upcase_iter}
\end{figure}

We can now sum up the inefficient deterministic implementation of the algorithm for the current case (Algorithm~\ref{alg:deterministic_2edges_up_bad}).

\begin{algorithm}[H]
  \begin{algorithmic}
  \caption{Inefficient deterministic solution for ``two tree edges, upper case''}
  \label{alg:deterministic_2edges_up_bad}
  \Function {SolveTwoTreeEdgesUpper} {}
  \For{$e \in E({\ensuremath{\mathcal{T}}})$}
    \State $f_0 \gets \DeepestDnCutNoMin(e)$
    \State $g \gets \mindn[1](e)$
    \State $i \gets 0$
    \While {$f_i \neq e$ and $\maxup[1](e).\mathrm{head} \in {\ensuremath{\mathcal{T}}}_e$}
      \State $h_i \gets \maxup[1](e)$
      \State $f_{i+1} \gets$ the tree edge with the same head as $h_i$
      \State $i \gets i + 1$
    \EndWhile
    \If {$f_i \neq e$}
      \State add $\{e, f_i, g\}$ to the list of $3$-edge-cuts
    \EndIf
  \EndFor

  \Comment {{\small Run the analogous algorithm for $\DeepestDnCutNoMax$ instead of $\DeepestDnCutNoMin$}.}
  \EndFunction

  \end{algorithmic}
\end{algorithm}

In order to optimize the algorithm, we seek to eliminate the $\mathsf{while}$ loop.
Indeed, for each $e \in E({\ensuremath{\mathcal{T}}})$, the process is very similar: start with some edge $f > e$, and then repeatedly replace $f$ with a~higher edge, until a~replacement would result in an~edge closer to the root than $e$.
We can model this process with a~rooted tree ${\ensuremath{\mathcal{U}}}$ defined as~follows:
\begin{itemize}
  \item the vertices of ${\ensuremath{\mathcal{U}}}$ are the tree edges of ${\ensuremath{\mathcal{T}}}$, and $\bot$,
  \item $\bot$ is the root of ${\ensuremath{\mathcal{U}}}$,
  \item for a~non-root vertex $e$ of ${\ensuremath{\mathcal{U}}}$, its parent is the tree edge with the same head as $\maxup[1](e)$.
    If the head of $\maxup[1](e)$ coincides with the root of ${\ensuremath{\mathcal{T}}}$, then the parent of $e$ in ${\ensuremath{\mathcal{U}}}$ is $\bot$.
\end{itemize}
Naturally, ${\ensuremath{\mathcal{U}}}$ can be constructed in linear time with respect to the size of ${\ensuremath{\mathcal{T}}}$; moreover, if an~edge $p$ is a~parent of another edge $q$ in ${\ensuremath{\mathcal{U}}}$, then $p <_{\ensuremath{\mathcal{T}}} q$.

Now, the $\mathsf{while}$ loop is equivalent to repeated replacement of $f_0$ with its parent in ${\ensuremath{\mathcal{U}}}$ as long as the parent is greater or equal than $e$ with respect to $<_{\ensuremath{\mathcal{T}}}$.
In other words, the final edge $f' := \mathsf{final}(f_0, e)$ is taken as the shallowest ancestor of $f_0$ in ${\ensuremath{\mathcal{U}}}$ for which $f' \geq_{\ensuremath{\mathcal{T}}} e$.

This leads to the final idea: we simulate a~forest of rooted subtrees of ${\ensuremath{\mathcal{U}}}$ using a~disjoint set union data structure $F_{\ensuremath{\mathcal{U}}}$ (Theorem~\ref{thm:linear_fu}).
Each subtree additionally keeps its root, which can be retrieved from $F_{\ensuremath{\mathcal{U}}}$ (Lemma~\ref{enrich-parent}).
Initially, each subtree of ${\ensuremath{\mathcal{U}}}$ contains a~single vertex.

After the initialization of $F_{\ensuremath{\mathcal{U}}}$, we iterate $e$ over the tree edges of ${\ensuremath{\mathcal{T}}}$ in the decreasing order of depth in ${\ensuremath{\mathcal{T}}}$.
Throughout the process, we maintain the following invariant on $F_{\ensuremath{\mathcal{U}}}$: an edge $ef \in E({\ensuremath{\mathcal{U}}})$ for $e < f$ has been added to the forest if and only if $e$ has been considered at any previous iteration of the main loop.
Hence, at the beginning of the iteration for a~given edge $e$, we add to $F_{\ensuremath{\mathcal{U}}}$ all tree edges of ${\ensuremath{\mathcal{U}}}$ originating from $e$.
At this point of time, for every tree edge $f$ such that $f >_{\ensuremath{\mathcal{T}}} e$, the edge $\mathsf{final}(f, e)$ is given by $F_{\ensuremath{\mathcal{U}}}.\mathsf{lowest}(f)$.
This reduces the entire $\mathsf{while}$ loop to a~single $\mathsf{lowest}$ query on $F_{\ensuremath{\mathcal{U}}}$ (Algorithm~\ref{alg:deterministic_2edges_up_good}).
\begin{algorithm}[H]
  \begin{algorithmic}
  \caption{Efficient deterministic solution for ``two tree edges, upper case''}
  \label{alg:deterministic_2edges_up_good}
  \Function {SolveTwoTreeEdgesUpper} {}
  \State ${\ensuremath{\mathcal{U}}} \gets$ the rooted tree on $E({\ensuremath{\mathcal{T}}}) \cup \{\bot\}$ defined above
  \State $F_{\ensuremath{\mathcal{U}}} \gets$ the disjoint set union data structure built on ${\ensuremath{\mathcal{U}}}$ (Lemma~\ref{lem:linear_fu_with_lowest})
  \For{$e \in E({\ensuremath{\mathcal{T}}})$, in order from the deepest to the shallowest}
    \For{$c$ -- child of $e$ in ${\ensuremath{\mathcal{U}}}$}
      \State $F_{\ensuremath{\mathcal{U}}}.\mathsf{union}(e, c)$
    \EndFor
    \State $g \gets \mindn[1](e)$
    \State $f_0 \gets \DeepestDnCutNoMin(e)$
    \State $f' \gets F_{\ensuremath{\mathcal{U}}}.\mathsf{lowest}(f_0)$
    \If {$f' \neq e$}
      \State add $\{e, f', g\}$ to the list of $3$-edge-cuts
    \EndIf
  \EndFor

  \Comment {{\small Run the analogous algorithm for $\DeepestDnCutNoMax$ instead of $\DeepestDnCutNoMin$}.}
  \EndFunction

  \end{algorithmic}
\end{algorithm}

It can be easily seen that we initialize $F_{\ensuremath{\mathcal{U}}}$ on a~tree with $n$ vertices, and we issue ${\ensuremath{\mathcal{O}}}(n)$ queries to it in total.
Therefore, the whole subroutine runs in time linear with respect to the size of $G$.

Summing up, we replaced each randomized subroutine with its deterministic counterpart, preserving the linear guarantee on the runtime of the algorithm.
We conclude that there exists a deterministic linear-time algorithm listing $3$-edge-cuts in $3$-edge-connected graphs.
\subsection{Previous work}

There has been a~plethora of research into the algorithms deciding the $k$-vertex- and $k$-edge-connectedness of graphs, and decomposing the graphs into $k$-vertex- or $k$-edge-connected components.
However, while classical, elementary, and efficient algorithms exist for $k = 1$ and $k = 2$, these problems become increasingly more difficult for the larger values of $k$.
In fact, even for $k = 4$, there were no known linear time algorithms to any of the considered problems.
The following description presents the previous work in this area for $k \in \{1, 2, 3, 4\}$, and exhibits the related work for the larger values of $k$:

\paragraph{$\mathbf{k = 1}$.}
Here, the notions of $k$-vertex-connectedness and $k$-edge-connectedness reduce to that of connectivity and the connected components of a~graph. In the static setting, determining the connected components in linear time is trivial.
As a~consequence, more focus is being laid on dynamic algorithms maintaining the connected components of graphs.
In the incremental setting, where the edges can only be added to the dynamic graph, the optimal solution is provided by disjoint-set data structures \cite{FU}, which solve the problem in the amortized ${\ensuremath{\mathcal{O}}}(\alpha(n))$ time per query, where $\alpha(n)$ denotes the inverse of the Ackermann's fast-growing function. The fully dynamic data structures are also considered
 \cite{Wulff, Sparsification, Cos2, Cos3, Cos4, ThorupDuzo, Cos6, DynConnWorst, Bridge, DynConnSODA17}.

\paragraph{$\mathbf{k = 2}$.}
One step further are the notions of 2-vertex-connectivity (biconnectivity)
and 2-edge-connectivity.
In the static setting, partitioning of a~graph into $2$-vertex-connected or $2$-edge-connected components are classical problems, both solved in linear time by exploiting the properties of the low function~\cite{2ConnStatic}.
The incremental versions of both problems are again solved optimally in the amortized ${\ensuremath{\mathcal{O}}}(\alpha(n))$ time per query \cite{Westbrook1992}.
Significant research has been done in the dynamic setting as well \cite{DBLP:conf/wads/PengSS19, ThorupDuzo, Sparsification, Bridge, Dyn3EConn, Henzinger1995, Dyn2ConnBack}.

\paragraph{$\mathbf{k = 3}$.} As a next step, we consider 3-vertex-connectivity (triconnectivity) and 3-edge-connectivity.
An~optimal, linear time algorithm detecting the $3$-vertex-connected components was first given by Hopcroft and Tarjan~\cite{TriConn}.
The first linear algorithm for $3$-edge-connectivity was discovered much later by Galil and Italiano~\cite{EdgeToVertex}, where they present a~linear time reduction from the $k$-edge-connectivity problem to $k$-vertex-connectivity for $k \geq 3$, showing that in the static setting, the former problem is the easier of the two.
This was later followed by a~series of works simplifying the solution for $3$-edge-connectivity~\cite{Tsin1, Tsin2, 3E1, Tsin3}.
The incremental setting \cite{3Incr, IncrSPQR} and the dynamic setting \cite{Dyn3EConn, Sparsification} were also considered.

We also mention that in the case of 3-vertex-connectivity, there exists a structure called SPQR-tree which succinctly captures the structure of 2-vertex-cuts in graphs~\cite{SPQR, IncrSPQR}.
Its edge-connectivity analogue also exists, but we defer its introduction to the general setting.

\paragraph{$\mathbf{k = 4}$.}
We move on to the problems of 4-vertex-connectivity and 4-edge-connectivity.
A~notable result by Kanevsky et al.~\cite{QuadConn} supports maintaining 4-vertex-connected components in incremental graphs, with an~optimal ${\ensuremath{\mathcal{O}}}(\alpha(n))$ amortized time per query. Their result also yields the solution for static graphs in ${\ensuremath{\mathcal{O}}}(m + n \alpha(n))$ time complexity.
By applying the result of Galil and Italiano~\cite{EdgeToVertex}, we derive a~static algorithm determining the $4$-edge-connected components in the same time complexity.
This algorithm is optimal for $m = \Omega(n \alpha(n))$.

Another result by Dinitz and Westbrook~\cite{Dinitz1998} supports maintaining the $4$-edge-connected components in the incremental setting. Their algorithm processes any sequence of queries in ${\ensuremath{\mathcal{O}}}(q + m + n \log n)$ time where $q$ is the number of queries, and $m$ is the total number of inserted edges to the graph.

However, it is striking that the fastest solutions for $4$-edge-connectivity and $4$-vertex-connectivity for static graphs were derived from the on-line algorithms working in the incremental setting.
In particular, no linear time algorithms for $k = 4$ were known before.

\paragraph{$\mathbf{k \geq 5}$.}
As a side note, we also present the current knowledge on the general problems of $k$-vertex-connectivity and $k$-edge-connectivity. A~series of results~\cite{Karger, DetMinCut, HenzingerCut, AnotherMinCut} show that it is possible to compute the minimum edge cut of a~graph (i.e., determine the edge-connectivity of a~graph) in near-linear time.
The previously mentioned work by Dinitz and Westbrook~\cite{Dinitz1998} maintains the $k$-edge-connected components of an~incremental graph which is assumed to already have been initialized with a~$(k-1)$-edge-connected graph.
The data structure answers any sequence of on-line queries in ${\ensuremath{\mathcal{O}}}(q + m + k^2 n \log{(n / k)})$ time, where $q$ is the number of queries, and $m$ is the number of edges in the initial graph.

Gomory and Hu \cite{GomoryHu} proved that for any weighted, undirected graph $G$ there exists a weighted, undirected tree $T$ on the same vertex set such that for any two vertices $s, t \in V(G)$, the value of the minimum $s$-$t$ edge cut in $T$ is equal to the value of the minimum $s$-$t$ edge cut in $G$.
Moreover, such a~tree can be constructed using $n-1$ invocations of the maximum flow algorithm.
In an interesting result by Hariharan et al.~\cite{PartialGomory-Hu}, the decomposition of any graph into $k$-edge-connected components is constructed in ${\ensuremath{\mathcal{O}}}((m + nk^3) \cdot \polylog(n))$ time, producing a~partial Gomory-Hu tree as its result.

Dinitz et el.~\cite{Cactus} showed that the set of all minimum edge cuts can be succinctly represented with a~\textit{cactus graph}. When the minimum edge cut is odd, this cactus simplifies to a tree (see \cite[Corollary 8]{SimplerCactus}). These results imply that if the size of the minimum cut is odd, then the number of minimum cuts is ${\ensuremath{\mathcal{O}}}(n)$ and if it is even, then the number of minimum cuts is ${\ensuremath{\mathcal{O}}}(n^2)$. The structure of $k$-vertex-cuts was also investigated \cite{Longhui}.

\subsection{Our results}
In this work, we present a linear time, deterministic algorithm partitioning static, undirected graphs into $4$-edge-connected components. Even though the area of the dynamic versions of the algorithms for $k$-edge-connectivity is still thriving, the progress in static variants appears to have plateaued. In particular, both subquadratic algorithms determining the $4$-edge-connected components \cite{QuadConn, Dinitz1998} originate from their dynamic incremental equivalents and are almost thirty years old, yet they did not achieve the optimal linear running time. Hence, our work
constitutes the first progress in the static setting of $4$-edge-connectivity in a~long time. As a~side result, our algorithm also produces the tree representation of $3$-edge-cuts as explained in~\cite{SimplerCactus}.

\subsection{Organization of the work}
The paper is organized as follows.
In Section~\ref{Smu1}, we show how to reduce the problem of determining $4$-edge-connected components to the problem of determining $4$-edge-connected components in $3$-edge-connected graphs.
In Section~\ref{sec:randomized_3cut}, we show a~linear time, randomized Monte Carlo algorithm for listing all 3-edge-cuts in 3-edge-connected graphs.
In Section \ref{sec:deterministic_3cut}, we show how to remove the dependency on the randomness in the algorithm from the previous section, producing a~linear time, deterministic algorithm listing all 3-edge-cuts in 3-edge-connected graphs.
Then, in Section \ref{Smu2}, we construct a tree of 3-edge-cuts in a~3-connected graph, given the list of all its 3-edge cuts. This tree is then used to determine the 4-edge-connected components of the graph.
Finally, in Section \ref{open-problems}, we present open problems related to this work.

\subsection{Zero tree edges}\label{0t}
\label{ssec:randomized_0tree}
As ${\ensuremath{\mathcal{T}}}$ is a spanning connected subgraph, 3-edge-cuts not containing any tree edges simply do not exist.
\subsection{One tree edge}\label{1t}
\label{ssec:randomized_1tree}
Let $e$ be an edge of ${\ensuremath{\mathcal{T}}}$ that is the only tree edge of some 3-edge-cut. The two resulting connected components after the removal of such a cut are ${\ensuremath{\mathcal{T}}}_e$ and ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$ (we are slightly abusing the notation here; we actually mean that the connected components are the subgraphs induced by the vertex sets of ${\ensuremath{\mathcal{T}}}_e$ and ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$, respectively). Since $e$ is not a bridge, $\low(e)$ is well-defined and it connects ${\ensuremath{\mathcal{T}}}_e$ with ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$, so it belongs to this cut as well. This uniquely determines the third edge of this cut as it has to be the unique edge whose hash is $H(e) \oplus H(\low(e))$.

In order to detect all such cuts, we iterate over all tree edges $e$ and for each of them we look up in $M$ whether there exists an edge whose compressed hash is $CH(e) \oplus CH(\low(e))$. If it exists and if it turns out to be a non-tree edge, then we call it $g$ and output the triple $\{e, \low(e), g\}$ as a 3-edge-cut.
\subsection{Two tree edges}\label{2t}

Let $e$ and $f$ be the edges of ${\ensuremath{\mathcal{T}}}$ that form a 3-edge-cut together with some back edge $g$.

\begin{lemma}
	The edges $e$ and $f$ are comparable with respect to $\le_{{\ensuremath{\mathcal{T}}}}$.
\end{lemma}
\begin{proof} Assume otherwise. Then ${\ensuremath{\mathcal{T}}}_e$ and ${\ensuremath{\mathcal{T}}}_f$ are disjoint, and $L \coloneqq {\ensuremath{\mathcal{T}}}_e \cup {\ensuremath{\mathcal{T}}}_f$ and $R \coloneqq {\ensuremath{\mathcal{T}}} \setminus L$ are two connected components resulting from the removal of the 3-edge-cut containing $e$ and $f$. However, $\low(e)$ and $\low(f)$ are well-defined different edges connecting $L$ and $R$, so any cut between $L$ and $R$ would need to contain $e, f, \low(e)$ and $\low(f)$ at the same time---a~contradiction. This proves that $e$ and $f$ are comparable.
	\end{proof}

Without loss of generality, assume that $e <_{\ensuremath{\mathcal{T}}} f$. Then, the two resulting connected components
 are $L \coloneqq {\ensuremath{\mathcal{T}}}_e \setminus {\ensuremath{\mathcal{T}}}_f$ and $R \coloneqq {\ensuremath{\mathcal{T}}} \setminus L$. The remaining back edge $g$ has to connect $L$ and $R$. We will distinguish two cases, differing in the location of the remaining back edge:
 \textit{two tree edges, lower case}, where $g$ connects $L$ with ${\ensuremath{\mathcal{T}}}_f$, and \textit{two tree edges, upper case}, where $g$ connects $L$ with ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$.

\begin{figure}[h]
	\centering
	\includegraphics[scale=0.4]{figures/randomized-two-tree-edges-cases.pdf}
	\caption{The settings in \ref{2tl} and \ref{2tu}. Cuts that we are looking for are formed by edges $\{e, f, g\}$ in both cases.}
	\label{fig:randomized-two-tree-edge-cases}
\end{figure}

\subsubsection{Two tree edges, lower case}\label{2tl}
\label{ssec:randomized_2tree_lo}

Let $A$ be the set of back edges between ${\ensuremath{\mathcal{T}}}_f$ and ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_f$. It consists of the edge $g$, which connects $L$ with ${\ensuremath{\mathcal{T}}}_f$, and of several edges connecting ${\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$. Let $B$ be the set of heads of edges from $A$ (we remind that back edges are directed towards the root $r$). All elements of $B$ lie on the path from $f$ to the root $r$, and the head of $g$ is the deepest element of $B$.

Based on this observation, we are going to use the data structure from Theorem \ref{thm:min_edge_path}. We initialize an instance of it with the tree ${\ensuremath{\mathcal{T}}}$, $k=1$, $C=2n$ and paths $\{P_e \ | \ e$ is a back edge$\}$, where for each non-tree edge $e=xy$, where $y$ is its head, we create one input path $P_e$ from $x$ to $y$, and we set $w_e \coloneqq 2n - {\rm pre}(y)$. This data structure lets us for each tree edge $f$ determine the back edge leaping over $f$ with the biggest value of preorder of its head. We call this edge $\maxup(f)$. This back edge is the only candidate for the edge $g$ for a~fixed edge $f$.

Hence, we can find all such cuts by firstly initializing the data structure, and then iterating over all tree edges $f$.
For each $f$, we take $g=\maxup(f)$. Knowing $f$ and $g$, we can look up in $M$ whether there exists an edge whose compressed hash is $CH(f) \oplus CH(g)$ (Lemma~\ref{lem:randomized_unique_remaining_edge}). If it exists and if it is a tree edge, then we call it $e$ and output the triple $\{e, f, g\}$ as a 3-edge-cut.

\subsubsection{Two tree edges, upper case}\label{2tu}
\label{ssec:randomized_2tree_hi}

Let $A$ be the set of non-tree edges between ${\ensuremath{\mathcal{T}}}_e$ and ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$. It consists of the back edge $g$, which connects $L$ with ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$, and of several edges connecting ${\ensuremath{\mathcal{T}}}_f$ with ${\ensuremath{\mathcal{T}}} \setminus {\ensuremath{\mathcal{T}}}_e$. Let $B$ be the set of the tails of the edges from $A$.
Let $v \in L$ be the tail of $g$. As $v \notin {\ensuremath{\mathcal{T}}}_f$, either for all vertices $u \in {\ensuremath{\mathcal{T}}}_f$ it holds that ${\rm pre}(v) < {\rm pre}(u)$, or for all vertices $u \in {\ensuremath{\mathcal{T}}}_f$ it holds that ${\rm pre}(u) < {\rm pre}(v)$.
In the first case, $v$ is the vertex with the smallest preorder which is a~tail of some edge leaping over $e$; while in the second case $v$ is the analogous vertex with largest preorder.

Based on this observation, we are going to use the data structure from Theorem \ref{thm:min_edge_path} again.
We initialize one instance of it with the tree ${\ensuremath{\mathcal{T}}}$, $k=1$, $C=2n$ and the paths $\{P_e \,\mid\, e\text{ is a back edge}\}$, where for each back edge $e=xy$, we create one input path $P_e$ from $x$ to $y$, and set $w_e \coloneqq {\rm pre}(x)$.
We also initialize another instance of this data structure in the same way, with the only difference that we set $w_e \coloneqq 2n - {\rm pre}(x)$ instead.
The first instance lets us for each tree edge $e$ determine the edge $\mindn(e)$ leaping over $e$ with the smallest preorder of its tail; while the second instance determines the analogous edge $\maxdn(e)$, only with the largest preorder.
By our considerations above, if any~desired $\{e, f, g\}$ cut exists, then $g \in \{\mindn(e), \maxdn(e)\}$.

Hence, we can find all such cuts by firstly initializing both instances of the data structure. Then, we iterate over all tree edges $e$, and for each $e$ we check both candidates for $g$. If our hash table contains an~edge whose compressed hash is $CH(e) \oplus CH(g)$ and it is a tree edge, then we call it $f$ and output the triple $\{e, f, g\}$ as a 3-edge-cut.

\subsection{Three tree edges}\label{3t}
\label{ssec:randomized_3tree}

Solving this case in a way similar to the previous cases seems intractable. We consider this subsection, together with the time analysis following it, as one of the key ideas of this work.

In this case, we assume that no non-tree edges belong to the cut.
Therefore, we can contract all of them simultaneously and recursively list all 3-edge-cuts in the resulting graph.
Since 3-edge-cuts in the contracted graph exactly correspond to the 3-edge-cuts consisting solely of tree edges in the original graph, this reduction is sound.
The contraction is performed in ${\ensuremath{\mathcal{O}}}(n + m)$ time by identifying the vertices within the connected components of $G - E({\ensuremath{\mathcal{T}}})$. Let $G'$ be the graph after these contractions. Since $G$ is 3-edge-connected, $G'$ has to be 3-edge connected as well, so the assumption about the input graph being 3-edge-connected is preserved. We do not modify the value of $b$ in the subsequent recursive calls, even though the value of $m$ will decrease.

\subsection{Time analysis}
\label{ssec:randomized_time_analysis}
We will now compute the expected time complexity of our algorithm. The subroutines detecting the cuts with at most two tree edges clearly take expected linear time. The only nontrivial part of the time analysis is the recursion in Subsection \ref{3t}. Let $T(m)$ denote the maximum expected time our algorithm needs to solve any graph with at most $m$ edges.
Since our graph is 3-edge-connected, the degree of each vertex in $G$ is at least three, so $|E(G)| \ge \frac{3}{2} |V(G)| \ge \frac{3}{2} |E(G')| \Rightarrow |E(G')| \le \frac{2}{3} |E(G)|$. The bound on $T(m)$ follows:
$$T(m) \le {\ensuremath{\mathcal{O}}}(m) + T\left(\frac{2}{3}m\right).$$
The solution to this recurrence is $T(m) = {\ensuremath{\mathcal{O}}}(m)$, hence the whole algorithm runs in expected linear time.

\subsection{Correctness analysis}
In this subsection, we are going to prove that our algorithm works with sufficiently high probability. The only reason it can output wrong result comes from the compression of hashes---if we used their uncompressed version instead, the algorithm would clearly be correct.

The maximum number of queries to our hash table is linear in terms of $n$, which follows from a similar argument to the one presented in Subsection~\ref{ssec:randomized_time_analysis}.
Hence, there exists some absolute constant $c$ such that the number of queries is bounded by $cn$.
For each query with value $q$, and for each edge whose hash is not equal to the hash whose compressed version we ask about, there is $2^{-b}$ probability that compressed version of this hash is equal to $q$.
Hence, the probability that we ever get a~false positive is bounded from above by $cnm2^{-b}$.
If we never get a~false positive, then our algorithm returns the correct output. Since $b = \ceil{3 \log_2(m)}$ and $n \le m$, we infer that $cnm2^{-b} \le \frac{c}{m}$.
Therefore, our algorithm works correctly with probability at least $1 - \frac{c}{m}$.

\subsection{Removing hash tables}
\label{ssec:removing_hash_tables}
As mentioned earlier, there is a way to avoid using hash tables in our algorithm. Since the values of compressed hashes of edges are polynomial in $m$ (at most $2m^3$), we are able to sort them in ${\ensuremath{\mathcal{O}}}(m)$ time using the radix sort.
Now, instead of requiring a data structure that allows us to query the existence of an edge with a particular value of its compressed hash, we create an object representing such query, which we are going to answer in the future in an offline manner. After gathering all such objects, we may sort them together with edges, where as keys for comparisons we use values of compressed hashes for edges and values of required compressed hashes for queries. All the keys for comparisons are polynomial, so we are able to sort them together in ${\ensuremath{\mathcal{O}}}(m)$ time using the radix sort. After the sorting, we are able to easily answer all queries offline.

This concludes the description of the randomized variant of the algorithm.


\section{Introduction}

Malaria is a mosquito-borne disease that continues to pose a heavy burden on South Sahara Africa~(SSA)~\cite{moran2007malaria}. Recently, there has been significant progress in improving treatment efficiency and reducing the mortality rate of malaria. Unfortunately, due to financial constraints, the policymakers face the challenge of ensuring continued success in disease control with insufficient resources. To make intelligent decisions, learning control policies over the years have been formulated as Reinforcement Learning~(RL) problems~\cite{bent2018novel}. Nevertheless, applying RL to malaria control seems to be a tricky issue since RL usually requires numerous trial-and-error searches to learn from scratch. Unlike simulation-based games, \emph{e.g., } Atari games~\cite{mnih2013playing} and game GO~\cite{silver2017mastering}, the endless intervention trial is unacceptable to regions over the years since the actual cost of life and money is enormous. Hence, as in many human-in-loop systems~\cite{zou2019longterm,zou2020neural,zou2020pseudo}, it is too expensive to apply RL to learn malaria intervention policy from scratch directly.

Therefore, to reduce the heavy burden of malaria in SSA, it is urgent to improve the data efficiency of policy learning. In~\cite{bent2018novel}, novel exploration techniques, such as Genetic Algorithm~\cite{holland1992genetic}, Batch Policy Gradient~\cite{sutton2000policy} and Upper/Lower Confidence Bound~\cite{auer2010ucb}, have been firstly applied to learn malaria control policies from scratch. However, these solutions are introduced under the Stochastic Multi-Armed Bandit~(SMAB) setting, which myopically ignores the delayed impact of the interventions in the future and might result in serious problems. For example, the large-scale use of spraying may lead to mosquito resistance and bring about the uncontrolled spread of malaria in the coming years. Hence, it requires us to optimize disease control policies in the long run, which is far more challenging than the exploration in SMAB.

Considering the long-term effects, the finite horizon continuous-space Markov Decision Process is employed to model the disease control in this work. Under this setting, we propose a framework named Variance-Bonus Monte Carlo Tree Search~(VB-MCTS) for data-efficient policy searching, illustrated in Figure~\ref{fig:solution}. Particularly, it is a model-based training framework, which iterates between updating the world model and collecting data. In model training, Gaussian Process~(GP) is used to approximate the state transition function with collected rollouts. As a non-parametric probabilistic model, GP can avoid the model bias and explicitly model the uncertainty about the transitions, \emph{i.e., }~the variance of the state. In data collection, we propose to employ MCTS for generating the policy with the mean MDP plus variance-bonus reward. The variance-bonus reward can decrease the uncertainty at the state-action pairs with high potential reward by explicitly motivating the agent to sample the state-actions with the highest upper-bounded reward. Furthermore, the sample complexity of the proposed method indicates that it is a PAC optimal exploration solution for malaria control. Finally, to verify the effectiveness of our policy search solution, extensive experiments are conducted on the malaria control simulators\footnote{\url{https://github.com/IBM/ushiriki-policy-engine-library}}~\cite{bent2018novel}, which are gym-like\footnote{\url{https://gym.openai.com}} environments for KDD Cup 2019~\cite{zhou2020kdd}. The outstanding performance on the competition and extensive experimental results demonstrated that our approach could achieve unprecedented data efficiency on malaria control compared to the state-of-the-art methods.

\begin{figure}
\centering
\includegraphics[width=2.7in]{figure/model.png}
    \caption{Overview of the proposed learning system on Malaria Control. The system
    alternates model training and data collection.}
    \label{fig:solution}
\end{figure}

Our main contributions are: \textbf{(1)} We propose a highly data-efficient learning framework for malaria control. Under the framework, the policymakers can successfully learn control policies from scratch within 20 rollouts. \textbf{(2)} We derive the sample complexity of the proposed method and verify that VB-MCTS is an efficient PAC-MDP algorithm. \textbf{(3)} Extensive experiments conducted on malaria control demonstrate that our solution can outperform the state-of-the-art methods.

\section{Related Work}
As a highly pathogenic disease, malaria has been widely studied from the perspective of predicting the disease spread, diagnosis, and personalized care planning. However, rarely work focuses on applying RL to learn the cost-effectiveness intervention strategies, which plays a crucial role in controlling the spread of malaria~\cite{moran2007malaria}. In~\cite{bent2018novel}, malaria control has firstly been formulated as a stochastic multi-armed bandit (SMAB) problem and solved with novel exploration techniques. Nevertheless, SMAB based solutions only myopically maximize instant rewards, and the ignorance of delayed influences might result in disease outbreaks in the future. Therefore, in this work, a comprehensive solution has been proposed to facilitate policy learning in a few trials under the setting of finite-horizon MDP.

Another topic is data-efficient RL. To increase the data efficiency, we are required to extract more information from available trials~\cite{deisenroth2011pilco}, which involves utilizing the samples in the most efficient way (\emph{e.g., } exploitation) and choosing the samples with more information (\emph{e.g., } exploration). Generally, for exploitation, model-based methods~\cite{ha2018recurrent,kamthe2017data} are more sample efficient but require more computation time for the planning. Model-free~\cite{szita2006learning,krause2016cma,van2009theoretical} methods are generally computationally light and can be applied without a planner, but need (sometimes exponentially) more samples, and are usually not efficient PAC-MDP algorithm~\cite{strehl2009reinforcement}. For exploration, there are two options: \textbf{(1)} Bayesian approaches, considering a distribution over possible models and acting to maximize expected reward; unluckily, it is intractable for all but very restricted cases, such as the linear policy assumption in PILCO~\cite{deisenroth2011pilco}. \textbf{(2)} intrinsically motivated exploration, implicitly negotiating the exploration/exploitation dilemma by always exploiting a modified reward for directly accomplishing exploration. However, on the one hand, the vast majority of papers only address the discrete state case, providing incremental improvements on the complexity bounds, such as ~MMDP-RB~\cite{sorg2012variance}, metric-E$^3$~\cite{kakade2003exploration}, and  BED~\cite{kolter2009near}. On the other hand, for more realistic continuous state space MDP, over-exploration has been introduced for achieving polynomial sample complexity in many work, such as KWIK~\cite{li2011knows}, and GP-Rmax~\cite{grande2014sample}. These methods will explore all regions equally until the reward function is highly accurate everywhere. By drawing on the strength of existing methods, our solution is a model-based RL framework, which efficiently plans with MCTS and trade-off exploitation and exploration by exploiting a variance-bonus reward.

\section{Proposed Method: VB-MCTS}
\subsection{Malaria Control as MDP}
Finding an optimal malaria control policy can be posted as a reinforcement learning task by illustrating it as a Markov Decision Process~(MDP).
Specifically, we formulate the task as a finite-horizon MDP, defined by the tuple $\langle\mathcal{S},\mathcal{A},P,R,\gamma\rangle$ with $\mathcal{S}$ as the potential infinite state space, $\mathcal{A}$ as the finite set of actions, $P:\mathcal{S}\times\mathcal{A} \rightarrow \mathcal{S}$ as the deterministic transition function,  $R:\mathcal{S}\rightarrow \mathbb{R}$ as the reward function, and $\gamma\in (0,1]$ as the discount factor.
In this case, we face the challenge of developing an efficient policy for a population over a 5 year intervention time frame.
As shown in Figure \ref{fig:mdp}, the corresponding components in malaria control are defined as,

\begin{figure}
\centering
\includegraphics[width=2.7in]{figure/mdp.png}
    \caption{Malaria intervention as MDP.}
    \label{fig:mdp}
\end{figure}

\paragraph{Action} The actions are the available means of interventions, including the mass-distribution of long-lasting insecticide-treated nets~(ITNs) and indoor residual spraying~(IRS) with pyrethroids in SSA~\cite{stuckey2014modeling}. In this work, the action space $\mathcal{A}$ is constructed through $a_i \in \mathcal{A} = \{(a_{ITN}, a_{IRS})\}$ with $a_{ITN},a_{IRS} \in (0,1]$, which represent the population coverage for ITNs and IRS of a specific area. Without significantly affecting performance, we discrete the action space with an accuracy of 0.1 for simplicity.

\paragraph{Reward} The reward is a scalar $r_t=R(s_{t+1})\in [R_{\min},R_{\max}]$, associated with next state $s_{t+1}$. In malaria control, it is determined through an economic cost-effectiveness analysis. In~\cite{bent2018novel}, an overview of the reward calculation is specified. Without loss of generality, the reward function $R(\ast)$ is assumed known to us since an MDP with unknown rewards and unknown transitions can be represented as an MDP with known rewards and unknown transitions by adding additional states to the system.

 \paragraph{State} The state contains important observations for decision making in every time step. In malaria control, it includes the number of life with disability, life expectancy, life lost, and treatment expenses~\cite{bent2018novel}. We set the state in the form $s_t = (r_{t-1},a_{t-1},t)$, including current reward $r_{t-1}$, previous action $a_{t-1}$ and the current intervention timestamp, which covers the most crucial and useful observations for malaria control, as shown in Figure \ref{fig:mdp}. For the start state $s_1$, the reward and action are initialized with 0.

Let $\pi:\mathcal{S}\rightarrow\mathcal{A}$ denote a deterministic mapping from states to actions, and let $V_{\pi}(s_t)=\mathbb{E}_{a_{i}\sim\pi(s_i)}\left[\sum_{i=t}^{T} \gamma^{i-t} r_i\right]$ denote the expected discounted reward by following policy $\pi$ in state $s_t$. The objective is to find a deterministic policy $\pi^\ast$ that maximizes the expected return at state $s_1$ as
\begin{eqnarray}\nonumber
    \pi^\ast = {\arg\max}_{\pi\in \Pi} V_{\pi}(s_1).
\end{eqnarray}

\subsection{Model-based Indirect Policy Search}
In the following, we detail the key components of the proposed framework VB-MCTS, including the world model, the planner, and variance-bonus reward with its sample complexity.

\paragraph{World Model Learning}
The probabilistic world model is implemented as a GP, where we use a predefined feature mapping $x_t = \phi(s_t,a_t) \in \mathbb{R}^{m}$ as training input and the target state $s_{t+1}$ as the training target. The GP yields one-step predictions
\begin{eqnarray}\nonumber
    p(s_{t+1}|s_t,a_t) &=& \mathcal{N} (s_{t+1}|\mu_t,\sigma^2_t)\\ \nonumber
    \mu_t &=& \bm{k}_{\ast}^\top(K+\omega^2_n I)^{-1}\bm{y},\\ \nonumber
    \sigma^2_t &=& k(x_{t},x_{t}) - \bm{k}_{\ast}^\top(K+\omega^2_n I)^{-1}\bm{k}_{\ast},
\end{eqnarray}
where $k$ is the kernel function, $\bm{k}_{\ast} \equiv k(\bm{X},x_{t})$ denotes the vector of covariances between the test point and all training points with $\bm{X}=[x_1,\cdots,x_n] $, and $\bm{y}=[s_1,\dots,s_n]^\top$ is the corresponding training targets. $\omega_n$ is the noise variance. $K$ is the Gram matrix with entries $K_{ij} = k(x_i,x_j)$.

Throughout this paper, we consider a prior mean function $m \equiv 0$ and a squared exponential~(SE) kernel with automatic relevance determination. The SE covariance function is defined as
\begin{eqnarray}\nonumber
     k(x,x') = \alpha^2 \exp (-\frac{1}{2}(x-x')^\top \mathbf{\Lambda}^{-1} (x-x')),
\end{eqnarray}
where $\alpha^2$ is the variance of state transition and $\bm{\Lambda} \equiv \operatorname{diag}\left(\left[l_{1}^{2}, \ldots, l_{m}^{2}\right]\right)$. The characteristic length-scale $l_i$ controls the importance of $i$-th feature. Given $n$ training inputs $\bm{X}$ and the corresponding targets $\bm{y}$, the posterior hyper-parameters of GP~(length-scales $l_i$ and signal variance $\alpha^2$) are determined through evidence maximization technique~\cite{williams2006gaussian}.

\paragraph{Exploration with Variance-Bonus Reward}

The algorithm we propose is itself very straightforward and similar to many previously proposed exploration
heuristics~\cite{kolter2009near,srinivas2009gaussian,sorg2012variance,grande2014computationally}. We call the algorithm Variance-Bonus Reward, since it chooses action according to the current mean estimation of the reward plus an additional variance-based reward bonus for state-actions that have not been well explored as
{\small
\begin{eqnarray} \nonumber
&& \tilde{R}(s_{t+1}|s_t,a_t) =
    R(s_{t+1})|_{s_{t+1} = \mathbb{E}_{s_{t+1}\sim GP(s_t,a_t)}[s_{t+1}]} + \\ \nonumber
&& \beta_1 \text{Var}_{s_{t+1}\sim GP(s_t,a_t)}[R(s_{t+1})]
    +	 \beta_2 \text{Var}_{s_{t+1}\sim GP(s_t,a_t)}[s_{t+1}],
\end{eqnarray}}
where $\beta_1$ and $\beta_2$ are the parameters that trade-off the balance of exploitation and exploration. $\text{Var}_{s_{t+1}\sim GP(s_t,a_t)}[R(s_{t+1})]$ and $\text{Var}_{s_{t+1}\sim GP(s_t,a_t)}[s_{t+1}]$ are the predicted variances for state and reward. The variance of reward can be exactly computed following the law of iterated variances.

\paragraph{Planning with Mean MDP + Reward Bonus}
MCTS is a strikingly successful planning algorithm~\cite{silver2017mastering}, which can find out the optimal solution with enough computation resources. Since disease control is not a real-time task, we propose to apply MCTS~(Figure~\ref{fig:solution}) as the planner for generating the policy to maximize the variance-bonus reward. In the executing process, the MCTS planner incrementally builds an asymmetric search tree guided to the most promising direction by a tree policy. This process usually consists of four phases --- \textbf{selection}, \textbf{expansion}, \textbf{evaluation}, and \textbf{backup}~(as shown in Figure~\ref{fig:mcts}).

Specifically, each edge $(s,a)$ of the search tree stores an average action value $\tilde{Q}(s, a)$ and visit count $N(s, a)$. In the \textbf{selection} phase, starting from the root state, the tree is traversed by simulation~(that is, descending the tree with the mean prediction of states without backup). At each time step $t$ of each simulation, an action $a_t$ is selected from state $s_t$
    \begin{eqnarray*}
        a_{t} = {\arg\max}_{a\in \mathcal{A}} ( \tilde{Q}(s_t,a)+ \frac{c_{puct}}{|\mathcal{A}|}\frac{\sqrt{\sum_{b\in \mathcal{A}}N(s_t,b)}}{1+N(s_t,a)} ),
    \end{eqnarray*}
so as maximize action value plus a bonus that decays with repeated visits to encourage exploration for tree search. Here, $c_{puct}$ is the constant determining the level of exploration. When the traversal reaches a leaf node $s_L$ at step $L$, the leaf node may be \textbf{expanded} with each edge initialized as $\tilde{Q}(s_L,a) = 0$, $N(s_L,a) = 0$, and the corresponding leaf nodes are initialized with the mean GP prediction $s_{L+1}=\mathbb{E}[s_{L+1}|s_L,a]$. Then, the leaf node is \textbf{evaluated} by the average outcome $v(s_L)$ of rollouts, which played out until terminal step $T$ using the fast rollout policies, such as random policy and greedy policy. At the end of the simulation, the action values and visit counts of all traversed edges are updated, \emph{i.e., } \textbf{backup}. Each edge accumulates the visit counts and means the evaluation of all simulations passing through that edge as
\begin{eqnarray*}
      \tilde{Q}(s_j,a_j) &\leftarrow & \frac{N(s_j,a_j)\times \tilde{Q}({s}_j,{a}_j)+v({{s}_j})}{N(s_j,a_j)+1},\\
  v({s}_j) &\leftarrow & v({{s}_{j+1}})+\tilde{R}(s_{j+1}|s_j,a_j),\\
  N(s_j,a_j) &\leftarrow &  N(s_j,a_j)+1,
\end{eqnarray*}
where $(s_j,a_j)$ with $j< L$ is the edge in the forward trace.

\begin{figure}
\includegraphics[width=2.7in]{figure/mcts.png}
    \caption{An illustration of MCTS for policy generation.}
    \label{fig:mcts}
\end{figure}

This cycle of selection, evaluation, expansion and backup is repeated until the maximum iteration number has been reached. At this point, the best action has been chosen by selecting the action that leads to the highest reward (max child) as follows
\begin{eqnarray}\label{equ:mcts_action}
    \pi(s_{t}) = \arg\max_{a\in\mathcal{A}} \tilde{Q}(s_t,a),
\end{eqnarray}
where $\pi(s_{t})$ is the policy generated by MCTS. To be noticed, the variance-bonus reward $\tilde{R}$ can be replaced with the $R$ for generating the best-known policy since it maximizes the expected reward based on the posterior so far.

\begin{algorithm}[tbh]
\small
\SetAlgoLined
\KwIn{Hype-parameters $\beta$, maximum trials number $N$.}
\KwOut{The malaria control policy $\bm{\pi}$.}
{\color{blue}\# Initialization.} \\
Generate random policy $\bm{\pi} = \left[a_1,\dots,a_T\right]$.\\
Apply policy $\bm{\pi}$ in real world and collect training samples $\mathcal{D}=\{(s_i,a_i,s_{i+1})\}_{i=1}^{T}$.\\
Update GP world model with $\mathcal{D}$.\\
Set number of trials $k=1$. \\
\For{$k<=N$}{
    Initialize time step $t=1$, and state $s_t$.\\
    \For{$t<=T$}{
    {\color{blue}\# Data Collection.} \\
    Generate action $a_t =\arg\max_{a\in\mathcal{A}} Q(s_t,a) $ by MCTS with \textbf{variance-bonus reward} $\tilde{R}$. \\
    Apply $a_t$ in real world and collect $s_{t+1}$.\\
    Update training samples $\mathcal{D} = \mathcal{D} \cup \{(s_t,a_t,s_{t+1})\}$.\\
    {\color{blue}\# Model Training.}\\
    Update GP world model with $\mathcal{D}$.\\
    Update $s_t=s_{t+1}$, and $t=t+1$.\\
    }
    Increase the training trials $k=k+1$.
}
{\color{blue}\# Generate final decision.}\\
Initialize time step $t=1$, state $s_t$, and policy $\bm{\pi}=\left[\right]$.\\
\For{$t<=T$}{
Generate action $a_t =\arg\max_{a\in\mathcal{A}} Q(s_t,a) $ by MCTS with \textbf{mean reward} $R$. \\
Predict $s_{t+1} = \mathbb{E}[s_{t+1}|s_t,a_{t}]$ with Gaussian Process.\\
Update $s_t=s_{t+1}$, and $t=t+1$.\\
Append action to final decision $\bm{\pi} = \bm{\pi}\oplus a_t$.
}
Return the malaria control policy $\bm{\pi}$.
\caption{Variance-Bonus Monte Carlo Tree Search}
\label{alg:training}
\end{algorithm}

Finally, we implement an iterative training procedure, as shown in Algorithm \ref{alg:training}, where we specify the order in which each component occurs within the iteration.

\subsection{Sample Complexity}
In the following, we derive sample complexity (\emph{e.g., } the required samples to learn near-optimal performances) for our proposed solution. In the worst case, when the reward function is equal everywhere, and all state-action pairs will be equally explored, VB-MCTS will have the same complexity bounds as over-exploration methods, such as GP-Rmax ~\cite{grande2014sample}. In practice, VB-MCTS will learn the optimal policy in fewer steps than the over-exploration methods because the high uncertainty region with low reward is not worth exploration. Following the general PAC-MDP theorem, Theorem 10 in~\cite{strehl2009reinforcement}, we derive the polynomial sample complexity of VB-MCTS.

\begin{theorem}\label{th:sample_complexity}
Assume that the feature space $\mathcal{X}\subset\mathbb{R}^{m}$ of state-action pairs is a compact domain and the target value is bounded $y\in [0, V_m]$. The reward and action value are Lipschitz continuous w.r.t the state-action pairs and let $L_q$ and $L_r$ be the Lipschitz constants for reward and action value respectively. If VB-MCTS is executed with $\beta_1 = \frac{L_r}{2 \sigma_{n}^{2}}$ and $\beta_2 = \frac{L_p}{2 \sigma_{n}^{2}}$ for any MDP $M$, then with probability $1-\delta$, VB-MCTS will follow a $4\epsilon$-optimal policy from its current state on all but
\begin{eqnarray}
    O\left(\frac{V_{m} \zeta(\epsilon, \delta)}{\epsilon(1-\gamma)} \log \frac{1}{\delta} \log \frac{1}{\epsilon(1-\gamma)}\right)
\end{eqnarray}
timesteps, with probability at least $1 - 2\delta$, where
$\zeta(\epsilon, \delta)=\left(\frac{4 V_{m}^{2}}{\epsilon^{2}(1-\gamma)^{2}} \log \left(\frac{2\left|\mathcal{N}_{\mathcal{X}}\left(\tau\left(\sigma_{\text {tol }}^{2}\right)\right)\right|}{\delta}\right)\right) \mathcal{N}_{\mathcal{X}}\left(\tau\left(\sigma_{\text {tol }}^{2}\right)\right)$,
$ \sigma_{\text{tol}}^{2} = \frac{2 \omega_{n}^{2} \epsilon_{1}^{2}}{V_{m}^{2} \log \left(\frac{2}{\delta_{1}}\right)}$. $\mathcal{N}_{\mathcal{X}}(\tau)$ is the covering number of domain $\mathcal{X}$, which defines the cardinality of the minimal set $C=\left\{c_{i}, \ldots, c_{N_{c}}\right\} \text { s.t. } \forall x \in \mathcal{X}, \exists c_{j} \in C \text { s.t. } d\left(x, c_{j}\right) \leq \tau,$ where $d(\cdot,\cdot)$ is a distance measure.
\end{theorem}

\begin{proof}[Sketch of Proof]
The proof is organized by showing the required key properties (optimism, accuracy, and learning complexity) for general PAC-MDP~\cite{strehl2009reinforcement} are satisfied.

\paragraph{Optimism} Assume that MCTS can find the optimal decision under every state $s$ with enough computation resources. Then, $\forall (s_t,a_t)$, the optimal action value $Q^\ast(s_t,a_t)$ satisfies
{\small
\begin{equation*}
\begin{aligned}
& Q^{*}\left(s_{t}, a_{t}\right)=\arg \max _{a_{t+1} \in \mathcal{A}} R\left(s_{t+1}\right)+\gamma Q^{*}\left(s_{t+1}, a_{t+1}\right) \\
\leq & R\left(s_{t+1}^{\prime}\right)+\frac{L_{r}}{\omega_{n}^{2}} \operatorname{Var}_{s_{t+1}^{\prime}\sim G P\left(s_{t}, a_{t}\right)} \left[R\left(s_{t+1}^{\prime}\right)\right] \\
+& \gamma \tilde{Q}^{*}\left(s_{t+1}^{\prime}, a_{t+1}\right)+\gamma \frac{L_{q}}{\omega_{n}^{2}} \operatorname{Var}_{s_{t+1}^{\prime}\sim G P\left(s_{t}, a_{t}\right)} \left[s_{t+1}^{\prime}\right] \\
\leq & \arg \max _{a_{t+1}^{\prime} \in \mathcal{A}} R\left(s_{t+1}^{\prime}\right)+\frac{L_{r}}{\omega_{n}^{2}} \operatorname{Var}_{s_{t+1}^{\prime}\sim G P\left(s_{t}, a_{t}\right)} \left[R\left(s_{t+1}^{\prime}\right)\right] \\
+& \gamma \tilde{Q}^{*}\left(s_{t+1}^{\prime}, a_{t+1}^{\prime}\right)+\gamma \frac{L_{q}}{\omega_{n}^{2}} \operatorname{Var}_{s_{t+1}^{\prime}\sim G P\left(s_{t}, a_{t}\right)} \left[s_{t+1}^{\prime}\right] \\
=& \tilde{Q}\left(s_{t}, a_{t}\right),
\end{aligned}
\end{equation*}
}
where $s_{t+1} {=} P(s_t,a_t)$, $s'_{t+1} {=}\mathbb{E}_{s_{t+1}\sim
GP(s_t,a_t)}[s_{t+1}|s_t,a_t]$.
The maximum error of propagating $Q(s_{t+1},a_{t+1})$ instead of $Q(s'_{t+1},a_{t+1})$ is given by $\|s'_{t+1}-s_{t+1}\|L_q$.
The norm $\|s'_{t+1}-s_{t+1}\|$ is upper bounded by the regularization error
$\frac{L_{q}}{\omega_{n}^{2}} \operatorname{Var}_{s_{t+1}^{\prime}\sim G P\left(s_{t}, a_{t}\right)} \left[s_{t+1}^{\prime}\right]$
with probability of $(1-\delta)$~(Section A.2 in ~\cite{grande2014computationally}).
In the same way, we have  $R(s'_{t+1}) + \frac{L_r}{\omega^2_n}\text{Var}_{s'_{t+1}\sim GP(s_t,a_t)}[R(s'_{t+1})] - R(s_{t+1})\geq 0$ with probability $1-\delta$.
It follows that Algorithm \ref{alg:training} remains optimistic with probability at least 1-2$\delta$.

\paragraph{Accuracy} Following the Lemma 1 from \cite{grande2014sample}, the prediction error at $x_i\in \mathcal{X}$ is bounded in probability $Pr\left\{\left|\hat{\mu}\left(x_{i}\right)-f\left(x_{i}\right)\right| \geq \epsilon_{1}\right\} \leq \delta_{1}$ if the predictive variance of the GP at $ \sigma_n^{2}\left(x_{i}\right) \leq \sigma_{\text{tol}}^{2}=\frac{2 \omega_{n}^{2} \epsilon_{1}^{2}}{V_{m}^{2} \log \left(\frac{2}{\delta_{1}}\right)}$.

\paragraph{Learning complexity} \textbf{(a)} For $\forall x\in\mathcal{X}$ with $\sigma_0^2 = 1$, let $\rho = \min\limits_{j=1,\dots,n} k(x,x_j)$, $d_{\max} = \max\limits_{j=1,\dots,n} d(x,x_j)$, then $d_{\max}$ can be rewritten as $d_{\max} = \sqrt{2\theta^2\log{\frac{1}{\rho}}}$ with $\rho =e^{-\frac{d_{\max}^2}{2\theta^2}}$. If $n\geq \frac{\omega_n^2}{\sigma_{\text{tol}}^2-1+\rho^2}$, the posterior variance of $x$ satisfies $\sigma_n^2(x)\leq 1-\frac{n\rho^2}{n+\omega^2_n}\leq\frac{n(1-p^2)+\omega^2_n}{n}\leq \sigma^2_{\text{tol}}$. \textbf{(b)} Since $\mathcal{X}$ is a compact domain with the length of $i$-th $\left(i\in [1,m]\right)$ dimension as $L_i$, we have $\mathcal{N}_{\mathcal{X}}\left(d_{\max }\right) \leq \mathcal{N}_{\mathcal{X}}\left(\frac{d_{\max }}{2}\right)$ and $\mathcal{X}$ can be covered by $2^{m} \prod_{j=1}^{m} L_{j} / d_{\max }^{m}$ balls $\left\{B\left(c_{j}, \frac{d_{\max }}{2}\right)\right\}_{j}$, which centers at $c_j$ with a radius of $\frac{d_{\max }}{2}$.\\

Given \textbf{(a)} and \textbf{(b)}, $\forall x\in \mathcal{X}$, $\exists c_j\in\mathcal{X}$ s.t. $B(x,d_{\max})\supset B(x,d_{\max}/2)$. If there are at least $\frac{\omega_n^2}{\sigma_{\text{tol}}^2-1+\rho^2}$ observations in $B(x,d_{\max}/2)$, then (2) satisfied. Combining the Lemma 8 from~\cite{strehl2009reinforcement} and denoting $\tau(\sigma^2_{\text{tol}}) = {d_{\max}}/2$, the total number of updates occurs will be bounded by $\zeta(\epsilon, \delta)=\left(\frac{4 V_{m}^{2}}{\epsilon^{2}(1-\gamma)^{2}} \log \left(\frac{2\left|\mathcal{N}_{\mathcal{X}}\left(\tau\left(\sigma_{\mathrm{tol}}^{2}\right)\right)\right|}{\delta}\right)\right) \mathcal{N}_{\mathcal{X}}\left(\tau\left(\sigma_{\mathrm{tol}}^{2}\right)\right)$ with probability of $1-\delta$.

Now that the key properties of optimism, accuracy, and learning complexity have been established, the general PAC-MDP Theorem 10 of~\cite{strehl2009reinforcement} is invoked.
\end{proof}

Theorem~\ref{th:sample_complexity} allows us to guarantee that the number of steps in which the performance of VB-MCTS is significantly worse than that of an optimal policy starting from the current state is at most log-linear in the covering number of the state-action space with probability $1-\delta$.

\section{Experiments}
In this section, we conduct extensive experiments on two different OpenMalaria~\cite{smith2008towards} based simulators: {\it SeqDecChallenge} and {\it ProveChallenge}, which are the testing environments used in KDD Cup 2019. In these simulators, the parameters of the simulator are hidden from the RL agents since the {\it“simulation parameters”} for SSA are unknown for policymakers. Additionally, to simulate the real disease control problem, only {\bf 20 trials} were allowed before generating the final decision, which is much more challenging than traditional RL tasks. These simulating environments are available at \url{https://github.com/IBM/ushiriki-policy-engine-library}.

\paragraph{Agents for Comparison}
To show the advantage of VB-MCTS, many benchmarking reinforcement learning methods and open source solutions have been deployed to verify its effectiveness: \textbf{Random Policy}: The random policy is executed in 20 trials and chooses the generated policy with the maximum reward as the final decision. \textbf{SMAB}: This kind of policy treats the problem as a Stochastic Multi-Armed Bandit problem and independently optimizes the policy every year with Thompson sampling~\cite{chapelle2011empirical}. \textbf{CEM}: Cross-Entropy Method is a simple gradient-free policy searching method~\cite{szita2006learning}. \textbf{CMA-ES}: CMA-ES is a gradient-free evolutionary approach to optimizing non-convex objective functions~\cite{krause2016cma}. \textbf{Q-learning-GA}: It learns the malaria control policy by combining Q-learning and Genetic Algorithm. \textbf{Expected-Sarsa}: It collects 13 random episodes and runs expected value SARSA~\cite{van2009theoretical} for 7 episodes to improve the best policy using the collected statistics. \textbf{GP-Rmax}: It uses GP learners to model $T$ and $R$ and replaces the value of any $Q^\ast (s, a)$ where $T (s, a)$ is {\it "unknown"} with a value of $\frac{R_{\max }}{1-\gamma}$~\cite{li2011knows,grande2014sample}. \textbf{GP-MC}: It employs Gaussian Process to regress the world model. The policy is generated by sampling from the posterior and choosing the max rewarded action. \textbf{VB-MCTS}: Our proposed method.

\paragraph{Implementation Details} We build a feature map of 14-dimension for this task, which includes the periodic feature and cross term feature. Specifically, the feature map is set as,
{\small\begin{eqnarray}
\nonumber
&&\phi(s_t,a_t) = [t,t\%2,t\%3,r_{t-1},a_{(t-1,{ITN})},a_{(t-1,{IRS})},\\ \nonumber
&& a_{(t,{ITN})},a_{(t,{IRS})},a_{(t,{ITN})}\times a_{(t,{IRS})},\\ \nonumber
&& a_{(t-1,{ITN})}\times a_{(t-1,{IRS})},a_{(t,{ITN})}\times a_{(t-1,{ITN})},  \\\nonumber
&& a_{(t,{IRS})}\times a_{(t-1,{IRS})},a_{(t,{ITN})}\times \left(1-a_{(t-1,{ITN})}\right),  \\ \nonumber
&& a_{(t,{IRS})}\times \left(1-a_{(t-1,{IRS})}\right)],
\end{eqnarray}}
where $t,t\%2,t\%3$ are the periodic features and $a_{(\ast,{\ast})}\times a_{(\ast,{\ast})}$ is the cross term feature. Since the predicted variance of state and reward are the same in our setting, we empirically set the sum of exploration/exploitation parameters as $\beta_1+\beta_2 = 3.5$ in the experiments. In MCTS, $c_{puct}$ is 5, and only the top 50 rewarded child nodes are expanded. The number of iterations does not exceed 100,000. For Gaussian Process, to avoid overfitting problems, 5-fold cross-validation is performed during the updates of the GP world model. Particularly, we use 1-fold for training and 4-fold for validation, which ensures the generalizability of the GP world model over different state-action pairs. Our implementation and all baseline codes are available at \url{https://github.com/zoulixin93/VB_MCTS}.

{\centering
\small
\begin{table}[!hpt]
\tabcolsep 0.04in
\caption{Performance comparisons between different agents.}
\label{tab:results}
\begin{tabular}{l|| c c c} \ChangeRT{1pt}
\centering
{\multirow{3}*{Agents}}
&\multicolumn{3}{c}{\textit{SeqDecChallenge}} \\
& Med. Reward & Max Reward & Min Reward  \\ \hline\hline
Random Policy & 167.79 & 193.24 & 135.06  \\
SMAB & 209.05 & 386.28 & -6.44 \\
\hline
CEM & 179.30 & 214.87 & 120.92 \\
CMA-ES & 185.34 & 246.12 & 108.18 \\
Q-Learning-GA & 247.75 & 332.40 & 171.33  \\
Expected-Sarsa & 462.76 & 495.03 & 423.93  \\
\hline
GP-Rmax & 233.95 & 292.99 & 200.35 \\
GP-MC & 475.99 & 499.60 & 435.51 \\
\rowcolor{lightgray}
VB-MCTS & \textbf{533.38} & \textbf{552.78} & \textbf{519.61}  \\
\ChangeRT{0.8pt}
{\multirow{3}*{Agents}}
&\multicolumn{3}{c}{\textit{ProveChallenge}} \\
& Med. Reward & Max Reward & Min Reward  \\ \hline\hline
Random Policy  & 248.25 & 464.92 & 55.24  \\
SMAB & 18.02 & 135.37 & -56.86 \\
\hline
CEM & 229.61 & 373.83 & 20.09 \\
CMA-ES & 289.03 & 314.57 & 92.95 \\
Q-Learning-GA & 242.97 & 325.24 & 88.70\\
Expected-Sarsa & 190.08 & 296.16 & 140.86 \\
\hline
GP-Rmax & 287.45 & 371.49 & 153.98 \\
GP-MC  & 300.37 & 447.15 & \textbf{263.96}\\
\rowcolor{lightgray}
VB-MCTS & \textbf{352.17} & \textbf{492.23} & 259.97 \\
\ChangeRT{0.8pt}
\end{tabular}
\end{table}}

\subsection{Results}
\paragraph{Main Results} In Table \ref{tab:results}, we report the median reward, the maximum reward, and the minimal reward over 10 independent repeat runs. The results are quite consistent with our intuition. We have the following observations: \textbf{(1)} For the finite-horizon decision-making problem, treating it as SMAB does not work, and the delayed influence of actions can not be ignored in malaria control. As presented in Table~\ref{tab:results}, SMAB's performances have large variance and are even worse than random policy in {\it ProveChallenge}. \textbf{(2)} Overall, two model-based methods~(GP-MC and VB-MCTS) consistently outperform the model-free methods~(CEM, CMA-ES, Q-learning-GA, and Expected-Sarsa) in {\it SeqDecChallenge} and {\it ProveChallenge}, which indicates that empirically model-based solutions are generally more data-efficient than model-free solutions. From the results in Table~\ref{tab:results}, performances of model-free methods are defeated by model-based methods with a large margin in {\it SeqDecChallenge}, and their performances are almost the same as the random policy in {\it ProveChallenge}. \textbf{(3)} The proposed VB-MCTS can outperform all the baselines in {\it SeqDecChallenge} and {\it ProveChallenge}. Compared with GP-MC and GP-Rmax, the advantage from the efficient MCTS with variance-bonus reward leads to the success on {\it SeqDecChallenge} and {\it ProveChallenge}.

{\begin{figure}[hbt]
\centering
\includegraphics[width=2.7in]{figure/learning_curves.png}
\caption{Malaria policy learning curves in five different agents.}
\label{fig:learning_curves}
\end{figure}}

\paragraph{Data Efficiency}
This paragraph compares data-efficiency~(required trials) of VB-MCTS with other RL methods that learn malaria policies from scratch. In Figure~\ref{fig:learning_curves}(a) and~\ref{fig:learning_curves}(b), we report agents' performances after collecting every trial episode in {\it SeqDecChallenge} and {\it ProveChallenge}. The horizontal axis indicates the number of trials. The vertical axis shows the average performance after collecting every episode. Figure \ref{fig:learning_curves}(a) and \ref{fig:learning_curves}(b) highlighted that our proposed VB-MCTS approach (brown) requires on average only 8 trials to achieve the best performances in {\it SeqDecChallenge} and {\it ProveChallenge}, including the first random trial. The results indicate that VB-MCTS can outperform the state-of-the-art method on both data efficiency and performance.

\section{Conclusion and Future Work}

We proposed a model-based approach employing the Gaussian Process to regress the state transition for data-efficient RL in malaria control. By planning with variance-bonus reward, our method can naturally deal with the dilemma of exploration and exploitation by efficiently MCTS planning. Extensive experiments conducted on the challenging malaria control task have demonstrated the advantage of VB-MCTS over state-of-the-arts both on performance and efficiency. However, the stationary setting of MDP may be unrealistic due to the development of disease control tools and the evolution of the disease. Therefore, data-efficient reinforcement learning under nonstationary settings will be more realistic and more challenging task.

\bibliographystyle{named}