# Non-greedy Tree-based Learning for Estimating Global Optimal Dynamic
Treatment Decision Rules with Continuous Treatment Dosage

Chang Wang, Lu Wang
Department of Biostatistics, University of Michigan

###### Abstract

Dynamic treatment regime (DTR) plays a critical role in precision medicine
when assigning patient-specific treatments at multiple stages and optimizing a
long term clinical outcome. However, most of existing work about DTRs have
been focused on categorical treatment scenarios, instead of continuous
treatment options. Also, the performances of regular black-box machine
learning methods and regular tree learning methods are lack of
interpretability and global optimality respectively. In this paper, we propose
a non-greedy global optimization method for dose search, namely Global Optimal
Dosage Tree-based learning method (GoDoTree), which combines a robust
estimation of the counterfactual outcome mean with an interpretable and non-
greedy decision tree for estimating the global optimal dynamic dosage
treatment regime in a multiple-stage setting. GoDoTree-Learning recursively
estimates how the counterfactual outcome mean depends on a continuous
treatment dosage using doubly robust estimators at each stage, and optimizes
the stage-specific decision tree in a non-greedy way. We conduct simulation
studies to evaluate the finite sample performance of the proposed method and
apply it to a real data application for optimal warfarin dose finding.

Keywords: dynamic treatment regime; dose finding; causal inference; non-greedy
tree-based learning; global optimality

##  1 Introduction

A dynamic treatment regime (DTR) is a sequence of decision rules that
determine the optimal treatment for individual patients at multiple stages.
Treatment decisions are based on each patient’s unique characteristics and
medical history to optimize their long-term clinical outcomes. With the
emergence of precision health care, DTRs with continuous dosage treatment have
become increasingly important, allowing for personalized optimal dosage
intervention. Examples include optimal dose finding in radiation oncology
therapy and drug trials with multiple stages.

Although black box learning methods like random forest and deep learning can
produce accurate predictions of the optimal treatments, their lack of
interpretability makes them difficult for medical experts to understand and
implement. Conventional tree learning methods such as CART are easy to
interpret and predict, but they often use greedy algorithms that can fail to
converge to the global optimum or achieve high performance under certain
circumstances. Additionally, while most research on dynamic treatment regimes
has focused on selecting the optimal treatment type, there is limited
knowledge on DTRs with continuous dose finding. Another challenge is how to
maintain desirable properties like doubly robustness and asymptotic normality
for continuous dose finding in DTR settings.

The literature on evaluating dynamic treatment regimes is extensive and can be
divided into two categories for finding the optimal tree-based DTR. One
approach is to use supervised learning methods, which involve estimating the
optimal treatment first and then transforming it into tree versions (known as
the treatment-tree algorithm). Several algorithms are commonly used for
optimal treatment estimation, including Q-learning (e.g. Murphy,, 2005),
Marginal Structural Models (e.g. Robins et al.,, 2000), and Outcome Weighted
Learning (e.g. Zhao et al.,, 2012). However, most of these methods require a
correct specification of working models for the propensity model or the
conditional outcome mean model, which can be challenging to validate with
limited knowledge. Therefore, a method with doubly robustness, which can
ensure consistent estimation as long as one of two models is correctly
specified, is desirable.

The other tree-based learning approach is searching for trees with the best
performance, also known as creating a new tree and evaluating its performance
under the optimal treatment (denoted as the tree-treatment algorithm).
However, this approach often relies on either random/stochastic search, which
lacks efficiency (e.g. Sun and Wang,, 2021), or greedy search, which may
converge to a local optimum and fail to find the global optimum as it only
maximizes the current purity measure (e.g. Tao et al.,, 2018). In other words,
both approaches have limitations, with the treatment-tree algorithm relying on
stochastic or greedy search and most treatment-tree algorithms not being
doubly robust. Therefore, we aim to propose a treatment-tree algorithm with
doubly robustness that avoids stochastic/greedy tree learning and applies a
non-greedy algorithm.

The second challenge arises when developing an individualized causal inference
algorithm for continuous treatment dosage with doubly robustness. The majority
of published papers about DTRs deal with categorical treatment assignment.
While there are some algorithms that can estimate doubly robust causal effects
and find optimal continuous treatment at the population level, these methods
require the use of a tree-treatment algorithm, which necessitates the decision
of the study population before the algorithms can be applied. As a result,
stochastic or greedy searching cannot be avoided. Alternatively, some methods
can estimate the optimal continuous treatment at the individual level by
optimizing a well-designed objective function. However, these methods are not
robust as some depend on the correct specification of the conditional outcome
model, while others rely on the correct specification of the propensity model
(Laber and Zhao,, 2015; Chen et al., 2016a, ). Therefore, we aim to propose a
new method for causal inference that can achieve individualized optimal dose
with doubly robustness.

To overcome the aforementioned challenges, we introduce a new method called
Global Optimal Dosage Tree-based learning (GoDoTree), which is a non-greedy
tree-based optimization technique that estimates optimal DTRs in a multi-stage
continuous-treatment environment, using data from randomized trials or
observational studies to provide personalized intervention with patient-
specific medicine dosage. At each stage, GoDoTree constructs a decision tree
by first modeling an individual-level counterfactual treatment effect curve
via semiparametric regression models, and then performing a non-greedy tree-
based search to optimize the counterfactual treatment effect. The two main
components of GoDoTree, namely individual-level counterfactual outcome
estimation and tree learning, are implemented in a backward manner at every
stage to optimize the long-term objective function. The proposed GoDoTree has
several advantages, including great interpretability, doubly robustness, and
the ability to achieve global optimal DTRs. It also contributes to the
existing literature on the development of DTRs with continuous treatment. To
demonstrate the performance of GoDoTree, we conduct simulation studies to show
its global optimization ability and apply it to warfarin dosing data collected
from Klein et al.,, 2009.

This paper is structured as follows. Section 2 introduces the proposed method
for individualized counterfactual outcome estimation. In Section 3, we
describe non-greedy tree-based supervised learning, which allows for a global
optimal tree search with known counterfactual outcome for every sample. The
details of finding the optimal kernel function are presented in Section 4.
Section 5 formalizes the problem of estimating optimal DTRs with continuous
treatment and outlines the framework of our algorithm. In Sections 6 and 7, we
present the results of numerical studies and an application example,
respectively. Finally, we provide a concluding discussion in Section 8.

##  2 Individualized Counterfactual Outcome Estimation With Continuous Dosage

We will firstly introduce treatment regime optimization with single stage. In
single stage scenario, we denote $Z_{i}=(\bm{X}_{i},A_{i},Y_{i})$ as the
observed data for patient $i$, where $\bm{X}_{i}$ is a vector of covariates,
$A_{i}$ is a continuous treatment or exposure, and $Y_{i}$ is the outcome of
interest. $Y^{a}$ is the counterfactual outcome when patient takes the
treatment $a$ and $\mu(\bm{x},a)=\mathbb{E}[Y^{a}\mid\bm{X}=\bm{x}]$ is the
conditional expectation of $Y^{a}$.

Population level effect curve
$\theta_{p}(a)=\mathbb{E}Y^{a}=\int_{\mathscr{X}}\mu(\bm{x},a)dP(\bm{x})$ is
frequently used in causal inference with continuous treatment, which is the
potential outcome that would have been observed under continuous treatment
level $A=a$ for the whole population. However, to perform global tree
optimization, we define a new individual-level effect curve $\theta_{i}(a)$
for the $i$ th patient as

$\theta_{i}(a)=\frac{\int_{\mathscr{X}}\mu(\bm{x},a)K_{i}(\bm{x})dP(\bm{x})}{\int_{\mathscr{X}}K_{i}(\bm{x})dP(\bm{x})}=\frac{\mathbb{E}\mu(\bm{X},a)K_{i}(\bm{X})}{\mathbb{E}K_{i}(\bm{X})},$

where $K_{i}(\bm{x})$ is a pre-specified kernel function centered at
$\bm{X}_{i}$ and gives samples different weights according to their
similarities with sample $i$. In another word, the neighborhood of sample $i$
are detected by $K_{i}$ and this subgroup’s weighted population level effect
curves are defined as the individual level effect curve for sample $i$.

Instead of optimizing the population level outcome, we propose a surrogate
objective function and optimize the individual level potential outcome:
$g^{opt}=\operatorname*{arg\,max}\ \sum_{i=1}^{n}\theta_{i}(g(\bm{X}))$.

Algorithm 1 Individual level effect curve estimation

1:Fit a conditional outcome mean model of $\mu(\bm{x},a)$ and a propensity
model of $\pi(a\lvert\bm{x})$. Get estimates $\hat{\mu}(\bm{x},a)$ and
$\hat{\pi}(a\lvert\bm{x})$, with $\sup\lvert\hat{\pi}-\pi^{*}\rvert=o_{p}(1)$
and $\sup\lvert\hat{\mu}-\mu^{*}\rvert=o_{p}(1)$.

2:In order to estimate the effect curve $\theta_{i}(a)$, we propose a doubly
robust estimator $\hat{\xi}(\bm{Z})$:
$\hat{\xi}(\bm{Z}{;}\hat{\pi},\hat{\mu},\hat{\kappa})=\hat{\kappa}^{-1}\left(\frac{Y-\hat{\mu}(\bm{X},A)}{\hat{\pi}(A\mid\bm{X})}\hat{w}(A)K_{i}(\bm{X})+\hat{m}(A)\right)$
where
$\hat{w}(a)=\int_{\mathcal{X}}\hat{\pi}(a\lvert\bm{x})d\mathbb{P}_{n}(\bm{x}),\hat{m}(a)=\int_{\mathcal{X}}K_{i}(\bm{x})\hat{\mu}(\bm{x},a)d\mathbb{P}_{n}(\bm{x}),\hat{\kappa}=\int_{\mathcal{X}}K_{i}(\bm{x})d\mathbb{P}_{n}(\bm{x}).$

3:Estimate patient-specific effect curve $\hat{\theta}_{i}(a)$ for every
patient $i$:
$\hat{\theta}_{i}(a)=\gamma_{ba}(a)^{T}\hat{\bm{D}}^{-1}_{ba}\mathbb{P}_{n}\\{\gamma_{ba}(A)K_{ba}(A)\hat{\xi}(\bm{Z}{;}\hat{\pi},\hat{\mu},\hat{\kappa})\\},$
where
$\hat{\bm{D}}^{-1}_{ba}=\mathbb{P}_{n}\\{\gamma_{ba}(a)K_{ba}(A)\gamma_{ba}(a)^{T}\\}$,
$\gamma_{ba}(A)=(1,(A-a)/b)^{T}$, $K_{ba}(A)$ is a kernel function centered at
$a$ and scaled with $b$.

###  2.1 Working Model Estimation

The estimations in the conditional outcome mean model and the propensity model
are vital for estimating optimal DTRs. Researchers can construct parametric
models with subject-relevant knowledge or just use some nonparametric models
like random forest, generalized additive model. Specifically, Bayesian
additive regression trees (BART) values as a great candidate because it
requires minimal specification of the model and can approximate complex
functions well with a tree ensemble (Chipman et al.,, 2010). Nonparametric
Kernel Smoothing Methods are also frequently used to fit propensity model
(Hayfield and Racine,, 2008).

###  2.2 Doubly Robust Estimator of the Individualized Effect Curve

We simplify the notation of $\theta_{i}(a)$ to $\theta(a)$ when it causes no
confusions, as we estimate effect curves $\theta_{i}(a),i=1,2,...,N$ for all
samples. To derive doubly robust estimators for $\theta(a)$, we adapt
semiparametric theory in a novel way similar to the approach of Rubin and van
der Laan,, 2005 and Kennedy et al.,, 2016. Our goal is to find a function
$\xi(\bm{Z}{;}\pi,\mu,\kappa)$ of the observed data $\bm{Z}$ and nuisance
functions $(\pi,\mu)$ with doubly robustness, which means
$\mathbb{E}\\{\xi(\bm{Z}{;}\pi^{*},\mu^{*},\kappa)\lvert A=a\\}=\theta(a)$ if
either $\pi^{*}=\pi$ or $\mu^{*}=\mu$ (no necessarily both). Here
$\kappa=\mathbb{E}K(\bm{X})$ is a constant when the kernel function $K$ is
fixed.

According to the semiparametric theory, a doubly robust mapping is related to
the efficient influence function for a certain parameter. If
$\mathbb{E}\\{\xi(\bm{Z}{;}\pi^{*},\mu^{*},\kappa)\lvert A=a\\}=\theta(a)$,
then it follows $\mathbb{E}\\{\xi(\bm{Z}{;}\pi^{*},\mu^{*},\kappa)\\}=\psi$
for $\psi=\mathbb{E}\theta(a)=\int_{\mathscr{A}}\theta(a)w(a)da$, where $w(a)$
is the marginal probability density function of treatment $A$. This indicates
that a component of the efficient influence function for the parameter $\psi$
can be a candidate for the doubly robust mapping
$\xi(\bm{Z}{;}\pi,\mu,\kappa)$(Robins and Rotnitzky,, 2001; van der Laan and
Robins,, 2003). We derive the efficient influence function for $\psi$ as
follows:

######  Theorem 2.1 (Efficient Influence Function)

Under a semiparametric model,
$\xi-\psi+(\mathbb{E}K_{i}(\bm{X}))^{-1}\int_{\mathscr{A}}[\mu(\bm{X},a)K_{i}(\bm{X})-m(a)]w(a)da$
is the efficient influence function for
$\psi=\int_{\mathscr{A}}\theta_{i}(a)w(a)da$, where
$w(a)=\int_{\mathcal{X}}\pi(a\lvert
x)dP(x),m(a)=\int_{\mathcal{X}}\mu(x,a)K_{i}(\bm{x})dP(x)$ and

$\xi(\bm{Z}{;}\pi,\mu)=(\mathbb{E}K_{i}(\bm{X}))^{-1}[\frac{Y-\mu(\bm{X},A)}{\pi(A\mid\bm{X})}w(A)K_{i}(\bm{X})+m(A)].$

Based on this, we propose a new estimator $\xi$, which is a component of the
efficient influence function of $\psi$, and further prove its doubly
robustness.

######  Theorem 2.2 (Doubly Robustness)

The proposed estimator $\xi$ has the doubly robust property: if $\pi^{*}=\pi$
or $\mu^{*}=\mu$, then we have

$\mathbb{E}\\{\xi(\bm{Z}{;}\pi^{*},\mu^{*},\kappa)\lvert A=a\\}=\theta(a).$

So as long as one of conditional outcome mean model $\mu$ or propensity model
$\pi$ is correctly specified (not necessarily both), our estimator $\xi$ will
be unbiased for $\theta(a)$.

###  2.3 Individual Level Effect Curve Estimation

With the doubly robust mapping $\xi(\bm{Z}{;}\pi,\mu,\kappa)$ derived in the
previous subsection, for which
$\theta_{i}(a)=\mathbb{E}\\{\xi(\bm{Z}{;}\pi^{*},\mu^{*},\kappa)\lvert A=a\\}$
as long as $\pi^{*}=\pi$ or $\mu^{*}=\mu$, we can construct doubly robust
estimate $\hat{\xi}(\bm{Z}{;}\hat{\pi},\hat{\mu},\hat{\kappa})$ and regress on
treatment variable $A$. The local linear kernel version of the estimator is
$\hat{\theta}_{b}(a)=\gamma_{ba}(a)\hat{\beta}_{b}(a)$, where

$\hat{\beta}_{b}(a)=\operatorname*{arg\,max}_{\beta\in\mathbb{R}^{2}}\mathbb{P}_{n}[K_{ba}(A)\\{\hat{\xi}(\bm{Z}{;}\hat{\pi},\hat{\mu},\hat{\kappa})-\gamma_{ba}(a)^{T}\beta\\}^{2}],\gamma_{ba}(A)=(1,(A-a)/b)^{T},$

for $K_{ba}(t)=b^{-1}K\\{(t-a)/b\\}$, with K a standard kernel function and
$b$ a scalar bandwidth parameter. Then we have the close form for estimation
of $\theta_{i}(a)$:

$\hat{\theta}_{i}(a)=\gamma_{ba}(a)^{T}\hat{\bm{D}}^{-1}_{ba}\mathbb{P}_{n}\\{\gamma_{ba}(A)K_{ba}(A)\hat{\xi}(\bm{Z}{;}\hat{\pi},\hat{\mu},\hat{\kappa})\\},$

where
$\hat{\bm{D}}^{-1}_{ba}=\mathbb{P}_{n}\\{\gamma_{ba}(a)K_{ba}(A)\gamma_{ba}(a)^{T}\\}$.

We will show the asymptotic bias and normality of the proposed individual
level effect curve $\hat{\theta}(a)$:

######  Theorem 2.3 (Asymptotic Bias Analysis)

Let $\pi^{*}$ and $\mu^{*}$ denote fixed functions to which $\hat{\pi}$ and
$\hat{\mu}$ converge in the sense that
$\sup\lvert\hat{\pi}-\pi^{*}\rvert=o_{p}(1)$ and
$\sup\lvert\hat{\mu}-\mu^{*}\rvert=o_{p}(1)$, and let $a\in\mathcal{A}$ denote
a point in the interior of the compact support $\mathcal{A}$ of treatment $A$.
Assume
$[\mu(\bm{X},t)-\hat{\mu}(\bm{X},t)][\pi(t\lvert\bm{X})-\hat{\pi}(t\lvert\bm{X})]=O_{p}(r_{n})$,
along with several regularity assumptions (see Appendix), we have:

(a) Either $\pi^{*}=\pi$ or $\mu^{*}=\mu$, where $\mu$ and $\pi$ are the true
conditional outcome model and propensity model respectively.

(b) The bandwidth $b=b_{n}$ satisfies $b\to 0$ and $nb^{3}\to\infty$ as
$n\to\infty$.

(c) $K$ is a continuous symmetric probability density.

(d) $\theta(a)$ is twice continuously differentiable, and both $\pi(a)$ and
the conditional density of $\xi(Z{;}\pi,\mu)$ given $A=a$ are continuous as
functions of $a$.

(e) The estimators $(\hat{\pi},\hat{\mu},\hat{\kappa})$ and their limits
$(\pi,\mu,\kappa)$ are contained in uniformly bounded function classes with
finite uniform entropy integrals (as defined in Section 4 of the Appendix),
with $1/\hat{\pi}$, $1/\pi$, $1/\hat{\kappa}$ and $1/\kappa$ also uniformly
bounded.

$\hat{\theta}_{b}(a)-\theta(a)=O_{p}(\frac{1}{\sqrt{nb}}+b^{2}+r_{n}).$

Here $O_{p}(r_{n})$ denotes the convergence rate of conditional outcome model
and propensity model. Since our first assumption is one of two models are
correctly specified, without loss of generality, let us assume the conditional
outcome mean model is correctly specified while propensity model not, i.e.,
$\pi^{*}\not=\pi$ and $\mu^{*}=\mu$. Then
$\hat{\mu}-\mu=\hat{\mu}-\mu^{*}=O_{p}(r_{n})$ and $\hat{\pi}-\pi=O_{p}(1)$,
because the latter one is biased but still bounded. Thus
$(\mu-\hat{\mu})(\pi-\hat{\pi})=O_{p}(r_{n})$ represent the convergence rate
of the correctly specified model, when one of them might be biased.

In the next theorem we show that when one or both of $\hat{\pi}$ and
$\hat{\mu}$ are estimated with fast enough convergence, then the proposed
estimator is asymptotically normal after scaling.

######  Theorem 2.4 (Asymptotic Normality)

Along with the same assumptions in theorem 2.3, also assume the convergence
rate $r_{n}$ satisfies $r_{n}=o_{p}(1/\sqrt{nb})$, then we have

$\sqrt{nb}\\{\hat{\theta}_{b}(a)-\theta(a)+bias(a)\\}\overset{d}{\to}N(0,\frac{\sigma^{2}(a)\int
K^{2}(u)du}{w(a)}),$

where $bias(a)=\theta^{\prime\prime}(a)(b^{2}/2)\int u^{2}K(u)du+o(b^{2})$,
and

$\displaystyle\sigma^{2}(a)$
$\displaystyle=\mathbb{E}\\{[\xi(\bm{Z}{;}\pi^{*},\mu^{*},\kappa)-\theta(a)]^{2}\mid
A=a\\}$
$\displaystyle=\kappa^{-2}\mathbb{E}\left[\frac{var\\{Y\lvert\bm{X},A=a\\}+\\{\mu(\bm{X},a)-\mu^{*}(\bm{X},a)\\}^{2}}{\\{\pi^{*}(a\lvert\bm{X})/w^{*}(a)\\}^{2}/\\{\pi(a\lvert\bm{X})/w(a)\\}}\right]-\\{\theta(a)-\kappa^{-1}m^{*}(a)\\}^{2}.$

With all assumptions satisfied, we can see that the mean square error of
$\hat{\theta}$ is $O(1/nb+b^{4})$ and the variance-bias trade-off determines
the theoretical optimal bandwidth $b\sim n^{-1/5}$. As for data-driven
bandwidth selection, we treated $\hat{\xi}$ as known and used leave-one-out
cross-validation bandwidth selection (Härdle et al.,, 1988):

$\hat{b}_{opt}=\operatorname*{arg\,max}_{b}\sum_{i=1}^{n}\left\\{\frac{\hat{\xi}(\bm{Z}_{i}{;}\hat{\pi},\hat{\mu})-\hat{\theta}_{b}(A_{i})}{1-\hat{W}_{b}(A_{i})}\right\\}^{2},$

where
$\hat{W}_{b}(a_{i})=(1,0)\mathbb{P}_{n}\\{\bm{\gamma}_{ba_{i}}(A)K_{ba_{i}}(A)\bm{\gamma}^{T}_{ba_{i}}(A)\\}^{-1}(1,0)^{T}b^{-1}K(0)$
is the $i^{th}$ diagonal element of the hat matrix.

##  3 Global Optimal Tree Search Algorithm

We use non-greedy tree-based learning to search for a global optimal decision
tree. Following last section’s algorithm, the individual level outcome
$\theta_{i}(a)$, $i\in\\{1,2,...,N\\}$, $a\in\mathscr{A}$, is estimated with
$\hat{\theta}_{i}(a)$ as the input of supervised learning. In this section, we
will simplify $\hat{\theta}_{i}(a)$ to $\theta_{i}(a)$ for brevity, because it
is treated as the known ground truth. Then the question can be formalized to
fit an optimal and interpretable continuous treatment decision tree
$g:\bm{X}\to A$, to maximize the objective function:
$L(g)=\sum_{i=1}^{N}\theta_{i}(g(\bm{X}))$.

Regular tree regression algorithms such as CART are developed in a greedy
manner. They grow from top-down without backtracking, determining all input
variables and split rules locally to optimize the current purity function.
However, greedy algorithms have potential limitations such as failing to
update parent nodes once they have a child node and the fact that the
cumulative local optimal rules do not necessarily lead to a global optimal
rule. Greedy algorithms are likely to miss stronger splits if they are hidden
behind weaker ones, resulting in sub-optimal or overly complicated trees
(Murthy and Steven,, 1995). To overcome these limitations, we use a non-greedy
algorithm called Tree Alternating Optimization (Perpinan and Tavallali,, 2018)
to search for the global optimum with iterative optimization. Given a decision
tree, we can optimize any node and keep the rest of the tree structure
unchanged. By alternating optimization over the depth levels of the tree, the
decision tree can optimize its structure, escape from local optima, and is
more likely to converge to the global optimum than if it were to keep
optimizing based on sub-optimal prior nodes. Figure 1 shows a comparison
between greedy trees and non-greedy trees.

#### Optimization at internal node

Consider an internal node (non-leaf) and the corresponding subset of samples
$\mathscr{S}$, our target is to find the optimal partition rule $\eta$ which
divides $\mathscr{S}$ into two parts $\omega$ and
$\mathscr{S}\backslash\omega$.

It is important to note that each element in $\mathscr{S}$ will eventually be
assigned to one of the children, and both children and their descendants are
fixed decision trees denoted as $g^{left}$ and $g^{right}$. If $\bm{X}$ is
classified to the left child, it will be assigned to treatment
$g^{left}(\bm{X})$, otherwise, it will be assigned to treatment
$g^{right}(\bm{X})$. The optimization of $\eta$ can be reduced to a supervised
classification problem where we seek the optimal rule $\eta$ that maximizes
the surrogate objective function and partitions the current sample subset
$\mathscr{S}$ into two parts $\omega$ and $\mathscr{S}\backslash\omega$:

$\max_{\eta}W(\eta)=\max_{\eta}\sum_{\begin{subarray}{c}\bm{X}_{i}\in\omega\end{subarray}}\theta_{i}\\{g^{left}(\bm{X}_{i})\\}+\sum_{\begin{subarray}{c}\bm{X}_{i}\in\mathscr{S}\backslash\omega\end{subarray}}\theta_{i}\\{g^{right}(\bm{X}_{i})\\}.$

The node update rule in GoDoTree is based on the idea of simulated annealing,
which allows us to avoid getting trapped at a local optimum by introducing a
probabilistic element into the decision-making process. Specifically, let
$\eta_{j},j=1,2,\dots,p$, denote the optimal splitting rule if the $j$ th
variable is fixed as the classifier to maximize $W(\eta)$ for the current
node, then the node update rule is
$P(\eta=\eta_{j})=exp\\{\alpha_{t}W(\eta_{j})\\}/\sum_{k=1}^{p}exp\\{\alpha_{t}W(\eta_{k})\\}$,
where $\alpha_{t}$ is a increasing sequence about the iteration number $t$.

Figure 1: Comparison between greedy tree and non-greedy tree algorithms.
Greedy tree optimizes the current objective function and cannot update nodes
once they are constructed, while the proposed approach, Tree Alternating
Optimization (TAO), is a non-greedy tree search that can update any node in
the tree at any time, with other nodes fixed. TAO can alternate optimization
over depth levels of the tree and escape from local optima. The classifier
used for partitioning is denoted as $\delta$.

#### Optimization at leaf node

As for the prediction part, our target is to optimize the parameter
$\eta=a\in\mathscr{A}$ so as to maximize the surrogate objective function
$W(\eta)=\sum_{\bm{X}_{i}\in\mathscr{S}}\theta_{i}(\eta)$. Samples classified
to the same leaf node will be assigned to the same dose $\eta$.

#### TAO algorithm framework

After updating the rule for each node, global optimization of the tree can be
achieved through alternating optimization over the depth levels of the tree.
The depth levels are cycled in the order root-leaf-root until convergence,
with the criteria for convergence being either a numerically converging
objective function or a fixed tree topology for several iterations. Pruning is
not performed until convergence, and when the decision tree converges, any
leaf node with a size smaller than a pre-specified number $n_{0}$ will be
collapsed with its sibling nodes.

Algorithm 2 Tree alternating optimization

1:Initialize a decision tree with height h, set t=0.

2:while not converge do

3: for i in 1 to h do $\triangleright$ root to leave

4: for every node in height i do

5: Find the optimal splitting rule $\eta_{k}$ when k th variable is used as
predictor.

6: Find current loss function W($\eta_{k}$).

7: Update node parameter $\eta$ with probability
$P(\eta=\eta_{j})=\frac{exp(\alpha_{t}W(\eta_{j}))}{\sum_{k=1}^{p}exp(\alpha_{t}W(\eta_{k}))}$.

8: for every leaf node do

9: Find the optimal treatment.

10: for i in h to 1 do $\triangleright$ leave to root

11: Do the same as row 4-7.

12: t = t+1

##  4 Implementation: Search For Optimal Kernel

Algorithm 3 Search for optimal kernel

Step 1: Initialize $\hat{\theta}_{i}(a)=\hat{\mu}(x_{i},a)$ using BART
(Bayesian additive regression tree).

Step 2: Define distance matrix
$D(i,j)=\sup_{a}\hat{\theta}_{i}(a)+\sup_{a}\hat{\theta}_{j}(a)-\sup_{a}\\{\hat{\theta}_{i}(a)+\hat{\theta}_{j}(a)\\}$.

Step 3: Calculate the similarity matrix $S_{ij}=cor(d_{i},d_{j})$, where
$d_{i}=(d_{i1},d_{i2},...,d_{in})$ is the $i$ th row of the distance matrix
$D$.

Step 4: Calculate another similarity matrix $\tilde{S}_{ij}$, using weighted
Euclidean distance $\tilde{D}$, where
$\tilde{D}_{ij}=\sum_{k=1}^{p}w_{k}(x_{ik}-x_{jk})^{2}$, $w_{k}$ is the
variable importance estimated in Step 1.

Step 5: Define new kernels as
$K_{i}(\textbf{x}_{j})=exp\\{\min(S_{ij},\tilde{S}_{ij})/\sigma_{i}^{2}\\}$,
where $\sigma_{i}^{2}$ is selected such that
$\sum_{j=1}^{n}K_{i}(\textbf{x}_{j})\approx n_{leaf}$, the expected sample
number assigned to tree’s leaves (e.g. $n/8$).

The idea of individualized effect curve estimation is, for a specific sample,
find its neighborhood by using a kernel and giving different weights, and
study the individualized effect curve based on this weighted sub-population.
Thus the key point is to find an optimal kernel function. An ill-defined
kernel will make the estimation of individual-level effect $\theta_{i}(a)$
seriously biased from true value.

We propose a new distance measurement to evaluate the similarity of effect
curves. For continuous treatment, the similarity between two samples should be
based on their effect curves, rather than simply optimal treatment doses. Our
newly proposed distance $D$ has two useful properties. Firstly, the distance
is relevant to horizontal translation, but immutable to vertical translation.
Secondly, smooth curves are close to any other curves under our definition. A
smooth effect curve means that the outcome is not strongly associated with the
treatments, thus this sample can be assigned to any treatments and regarded as
many samples’ neighborhood. So, a variable will not influence the distance if
it only has main effect but no interaction with treatment. $\tilde{D}$ is
weighted to get rid of low information signal-noise rate. Because we only care
about the optimal dose, the variable importance $w_{k}$ here is not about the
main effect, but the interaction between variable and treatment. Finally, we
incorporate two similarities with $\min(S_{ij},\tilde{S}_{ij})$ to search for
neighborhood with both effect curve and covariate similarities. Users are also
allowed to use other similarity definitions, e.g.
$\max(S_{ij},\tilde{S}_{ij})$, which will search for neighborhood with either
effect curve or covariate similarities.

To illustrate the distance definition, we demonstrate four example effect
curves and their distances(See Figure 2 for visualization). The effect curves
of A and B have the same shape but different horizontal location/optimal dose,
thus their distance is large, which is denoted by the vertical line. A and C
have different shapes and locations, but C is a smooth curve. So we can regard
C as anyone’s neighborhood and distance between A and C is relatively small. A
and D have different vertical location but the same shape, and their distance
is zero.

Figure 2: New distance visualization with
$D(i,j)=\sup_{a}\theta_{i}(a)+\sup_{a}\theta_{j}(a)-\sup_{a}\\{\theta_{i}(a)+\theta_{j}(a)\\}$.
Four effect curves are used as examples in the left figure. Three curves in
the right figure represent the component $\theta_{i}(a)+\theta_{j}(a)$ and the
vertical lines in the right figure visualize the pairwise distance measure.

##  5 Optimal Dose Finding in Dynamic Treatment Regimes

Suppose the data $\\{(\mathbb{X}_{t},A_{t},R_{t})_{t=1}^{T}\\}$ is independent
and identically distributed with sample size $n$ and comes from either a
randomized trial or an observational study, where $t\in\\{1,2,\dots,T\\}$
denotes the $t^{th}$ stage, $\mathbb{X_{t}}$ denotes the patient
characteristics during $t^{th}$ stage, $A_{t}\in\mathscr{A}_{T}$ denotes a
bounded continuous treatment variable and $R_{t}$ denotes the reward of
current stage following $A_{t}$. Let $\mathbb{H}_{t}$ denote patient history
before treatment assignment $A_{t}$, i.e.,
$\mathbb{H}_{t}=\\{(\mathbb{X}_{v},A_{v},R_{v})_{v=1}^{t-1},\mathbb{X_{t}}\\}$.
We consider the long term outcome of interest as $Y=\psi(R_{1},\dots,R_{T})$,
where $\psi$ is a prespecified function (e.g., sum of $R_{i}$ or last value
$R_{T}$).

We denote a DTR, a sequence of individualized treatment decisions, as
$g=(g_{1},\dots,g_{T})$, where $g_{t}$ maps from patient history
$\mathbb{H}_{t}$ to a continuous treatment $A_{t}$. To define an optimal DTR,
we use the counterfactual outcome framework of causal inference in a backward
way. At the final stage $T$, let $Y^{*}(A_{1},\dots,A_{T-1},a_{T})$, or
$Y^{*}(a_{T})$ for brevity, denotes the counterfactual outcome had a patient
been treated with $A_{T}=a_{T}$ conditional on previous treatments
$(A_{1},...,A_{T-1})$, and define $Y^{*}(g_{T}):=Y^{*}(g_{T}(\mathbb{H}_{T}))$
as the counterfactual outcome under the regime $g_{T}$.

The performance of $g_{T}$ can be evaluated by the value function $V(g_{T})$
(Qian and Murphy,, 2011), which is defined as the mean counterfactual outcome
had all patients followed $g_{T}$, i.e., $V(g_{T}):=E{Y^{*}(g_{T})}$.
Therefore, the optimal rule $g_{opt}$ should satisfy $V(g_{opt})\geq V(g_{T})$
for all $g_{T}\in\mathscr{G}_{T}$, where $\mathscr{G}_{T}$ denotes the set of
all possible rules of interest. To identify the optimal DTRs, we make the
standard assumptions to link the distribution law of counterfactual data with
that of observational data (Murphy et al.,, 2001).

First, we assume the consistency, i.e., the observed outcome is the same as
the counterfactual outcome under the assigned treatment, i.e.,
$Y=Y^{*}(A_{T})$. Secondly, we assume no unmeasured confounding assumption
(NUCA), i.e., $A_{T}\perp\\!\\!\\!\perp\mathscr{Y}_{T}\mid\mathbb{H}_{T}$,
where $\mathscr{Y}_{T}=\\{Y^{*}(a_{T}):a_{T}\in\mathscr{A}_{T}\\}$ and
$\perp\\!\\!\\!\perp$ denotes statistical independence. We also assume
positivity of the probability density function $f_{A_{T}}$, i.e.,
$f_{A_{T}}(a)>\epsilon>0$, $\forall a\in\mathscr{A}_{T}$. Finally, we assume
the continuity of the counterfactual outcome mean, i.e.,
$\mathbb{E}Y^{*}(a_{T})$ is continuous about $a_{T}\in\mathscr{A}_{T}$. Under
these assumptions, the optimal rule at stage T can be written as

$g_{T}^{opt}=\operatorname*{arg\,max}_{g_{T}\in\mathscr{G}_{T}}E_{\mathbb{H}_{T}}\Big{[}E[Y\mid
A_{T}=g_{T}(\mathbb{H}_{T}),\mathbb{H}_{T}]\Big{]}.$

At an intermediate stage $t$ ($1\leq t\leq T-1$), we consider
$Y^{*}(A_{1},A_{2},\dots,g_{t},g_{t+1}^{opt},\dots,g_{T}^{opt})$, the
counterfactual outcome under optimal rules for all future stages, had a
patient following $g_{t}$ at stage $t$, given $A_{1},\dots,A_{t-1}$ (Moodie et
al.,, 2012). Similarly, under the four assumptions above, the optimal rule
$g^{opt}$ at stage t can be defined as

$\displaystyle g_{t}^{opt}=$
$\displaystyle\operatorname*{arg\,max}_{g_{t}\in\mathscr{G}_{t}}E_{\mathbb{H}_{t}}[Y^{*}(A_{1},A_{2},\dots,g_{t},g_{t+1}^{opt},\dots,g_{T}^{opt})]$
$\displaystyle=$
$\displaystyle\operatorname*{arg\,max}_{g_{t}\in\mathscr{G}_{t}}E_{\mathbb{H}_{t}}\Big{[}E[\tilde{Y}_{t}\mid
A_{t}=g_{t}(\mathbb{H}_{t}),\mathbb{H}_{t}]\Big{]},$

where $\mathscr{G}_{t}$ is the set of all potential rules at stage $t$,
$\tilde{Y}_{T}=Y$ at stage $T$, and can be defined recursively using Bellman’s
optimality at an earlier stage $t$:

$\tilde{Y}_{t}=E\Big{[}\tilde{Y}_{t+1}\mid
A_{t+1}=g_{t+1}^{opt}(\mathbb{H}_{t+1}),\mathbb{H}_{t+1}\Big{]}.$

A dynamic treatment regime $g$ consists of treatment rules at all treatment
stages $t\in\\{1,2,\dots,T\\}$. The aim is to find the optimal treatment
decision rules $g^{opt}=(g_{1}^{opt},\dots,g_{T}^{opt})$, such that when
followed by the targeted patient population, the average response outcome $Y$
conditional on individual history is optimized. In the above sections we have
already introduced treatment regime optimization with single stage. When it is
extended to multistage scenario, for each stage t, we just need to replace the
outcome $Y$ with $\tilde{PO}_{t}(a_{t}):=\hat{E}(\tilde{Y}_{t}\mid
A_{t}=a_{t},\mathbb{H}_{t})$ and replace the covariates $\mathbb{X}$ with
patient history $\mathbb{H}_{t}$.

##  6 Simulations

For scenario 1, we use a regular setting $A_{opt}=(x_{1}+x_{2})/2$ (Laber and
Zhao,, 2015), to test different models’ performance when optimal treatment is
continuous. We consider generative models in which treatments are uniformly
distributed on (0,1), covariates $X$ are uniformly distributed on the
p-dimensional unit cube $[0,1]^{p}$ and $Y=u(X)-c(X,A)+Z$,
$u(X)=k_{p}+\tau_{p}\sum_{j=1}^{p}X_{j}$, where $k_{p}$ and $\tau_{p}$ are
chosen so that $var\\{u(X)\\}=5$ and $E\\{u(X)\\}=-E\\{\sup_{a}c(X,a)\\}$, Z
is an independent standard normal variate, and

$c(x,a)\propto\frac{1}{1+10(2a-x_{1}-x_{2})^{2}}.$

For scenario 2, we consider generative model in which treatments are uniformly
distributed on (0,1), covariates $X$ are uniformly distributed on the
p-dimensional unit cube $[-1,1]^{p}$ and $Y\sim\mu(X,A)+N(0,1)$,
$\mu(X,A)=(\sum_{i=1}^{p}x_{i})/p+100(A-A_{opt})^{2}$. To illustrate the
global optimality of GoDoTree, we consider a simple scenario where greedy
algorithms will fail: $A_{opt}=0.75$ when $x_{1}x_{2}\geq 0$ and
$A_{opt}=0.25$ when $x_{1}x_{2}<0$. In this case, the interaction of $x_{1}$
and $x_{2}$ decides the pattern of optimal dose. Although $x1$ and $x_{2}$ are
informative variables, greedy algorithms will fail in this scenario since
there is no purity increment when only one of them is used as the splitting
rule.

For scenario 3 and 4, we consider a generative model in two-stage scenario.
Two-stage treatments $A_{1}$ and $A_{2}$ are uniformly distributed on (0,1),
two-stage covariates $X$ and $Z$ are both uniformly distributed on the
p-dimensional unit cube $[-1,1]^{p}$ and two-stage outcomes are generated as
$Y_{1}\sim\mu_{1}(X,A_{1})+N(0,0.1)$, $Y_{2}\sim\mu_{2}(Z,A_{2})+N(0,0.1)$,
where $\mu_{1}(X,A_{1})=(\sum_{i=1}^{p}X_{i})/p+\rho_{1}(A_{1}-A_{opt,1})^{2}$
and $\mu_{2}(Z,A_{2})=(\sum_{i=1}^{p}z_{i})/p+\rho_{2}(A_{2}-A_{opt,2})^{2}$.
For scenario 3, we consider linear optimal treatments which can be detected by
greedy algorithm:$\rho_{1}=1,\rho_{2}=2$, $g_{opt,1}(X)=0.5+(X_{1}+X_{2})/4$
and $g_{opt,2}(H_{1},Z)=0.5+(Y_{1}+Z_{1})/4$. For scenario 4, we consider a
tree type optimal treatment at stage 2 which cannot be detected by greedy
algorithm: $\rho_{1}=2,\rho_{2}=10$, $g_{opt,1}(X)=0.5+(X_{1}+X_{2})/4$,
$g_{opt,2}(H_{1},Z)=0.2$, when $Z_{1}(Y_{1}-0.1)>0$, and $g_{opt,2}=0.8$
otherwise. Here $H_{1}=(X,A_{1},Y_{1})$ denotes the history data before
$A_{1}$, the optimal dynamic treatment decision rule
$g_{opt}=(g_{opt,1},g_{opt,2})$ is to minimize the long term outcome
$Y=Y_{1}+Y_{2}$.

Three algorithms for optimal dose finding are compared with the proposed
GoDoTree algorithm: LZ (Laber and Zhao,, 2015), CART, and random dosing. LZ
uses outcome weighted learning and a greedy approach to search for the optimal
tree. CART is trained by estimates of the conditional outcome mean model to
predict the optimal dose. Random dosing is used as a benchmark. To make the
comparison fair, all algorithms use BART as the working conditional outcome
mean model, and the propensity model is assumed to be known for LZ but unknown
for GoDoTree. LZ relies on a correctly specified propensity model, while
GoDoTree has doubly robustness and only requires one of the two models to be
correctly specified.

Table 1: Comparison of the proposed method GoDoTree with other methods for determining the optimal continuous treatment $A_{opt}=(x_{1}+x_{2})/2$: n = 500, p = (10,50), $E\\{Y(g_{opt})\\}=0$; reported values are based on 100 Monte Carlo replications, using a test set of size 1000. Scenario 1 | Tree Height = 2 | Tree Height = 3
---|---|---
$E\\{Y(\hat{g}_{opt})\\}$ | RMSE of $\hat{g}_{opt}$ | $E\\{Y(\hat{g}_{opt})\\}$ | RMSE of $\hat{g}_{opt}$
p=10 | GoDoTree | 2.20 (0.19) | 0.125 (0.01) | 1.81 (0.23) | 0.112 (0.01)
LZ | 2.07 (0.27) | 0.125 (0.01) | 1.90 (0.23) | 0.118 (0.01)
CART | 2.22 (0.19) | 0.136 (0.01) | 2.07 (0.20) | 0.128 (0.01)
Random | 4.51 (0.09) | 0.353 (0.01) | 4.51 (0.09) | 0.353 (0.01)
p=50 | GoDoTree | 2.27 (0.25) | 0.138 (0.01) | 2.06 (0.29) | 0.127 (0.02)
LZ | 2.12 (0.25) | 0.128 (0.01) | 2.12 (0.24) | 0.133 (0.02)
CART | 2.32 (0.19) | 0.140 (0.01) | 2.26 (0.23) | 0.142 (0.01)
Random | 4.52 (0.10) | 0.353 (0.01) | 4.52 (0.10) | 0.353 (0.01)
Table 2: Simulation results in scenario 2, where greedy search cannot achieve global optimum: n = 500, p = (10,50), $E\\{Y(g_{opt})\\}=0$; reported values are based on 100 Monte Carlo replications, using a test set of size 1000. Scenario 2 | Tree Height = 2 | Tree Height = 3
---|---|---
$E\\{Y(\hat{g}_{opt})\\}$ | RMSE of $\hat{g}_{opt}$ | $E\\{Y(\hat{g}_{opt})\\}$ | RMSE of $\hat{g}_{opt}$
p=10 | GoDoTree | 2.51 (1.29) | 0.15 (0.04) | 2.35 (0.84) | 0.15 (0.02)
LZ | 7.91 (2.31) | 0.65 (0.04) | 7.00 (3.38) | 0.66 (0.04)
CART | 5.71 (0.71) | 0.24 (0.02) | 5.26 (0.96) | 0.23 (0.02)
Random | 14.55 (0.52) | 0.38 (0.01) | 14.55 (0.52) | 0.38 (0.01)
p=50 | GoDoTree | 4.54 (1.19) | 0.21 (0.03) | 4.34 (1.05) | 0.21 (0.03)
LZ | 10.95 (1.85) | 0.66 (0.04) | 8.43 (1.45) | 0.67 (0.05)
CART | 6.24 (0.24) | 0.25 (0.004) | 6.24 (0.27) | 0.25 (0.005)
Random | 14.53 (0.54) | 0.38 (0.01) | 14.53 (0.54) | 0.38 (0.01)

In Tables 1 and 2, we report the average performances and corresponding
standard deviations of GoDoTree and several other methods for the single-stage
scenario, using two evaluation criteria: $E\\{Y(\hat{g}_{opt})\\}$ (the
expected counterfactual outcome when patients take the decision rule
$\hat{g}_{opt}$, with an lower bound of $E\\{Y(g_{opt})\\}=0$ for both
scenarios) and root mean square error (RMSE, the square root of
$E\\{g_{opt}-\hat{g}_{opt}\\}^{2}$, measuring the difference between the
proposed decision rule $\hat{g}_{opt}$ and the optimal decision rule
$g_{opt}$). The results in Table 1 show that both LZ and GoDoTree outperform
CART, and that the decision trees are robust with additional noise variables
and varying tree height. In Table 2, we find that GoDoTree achieves the best
performance among all methods: it is the only method that recognizes the
global pattern (i.e., the interaction between $x_{1}$ and $x_{2}$), and thus
has the highest counterfactual outcome mean and the smallest RMSE for
$\hat{g}_{opt}$. In contrast, CART and LZ fail to recognize the global pattern
and get trapped in local optima due to their greedy searching approach.

Table 3 shows the average performances of the above methods in two-stage
scenarios, which are consistent with the results of the single-stage scenario.
In scenario 3, where the optimal treatment structure is simple and can be
identified by greedy algorithms, all three methods have comparable results. In
scenario 4, where the optimal treatment has an underlying tree structure that
cannot be detected by greedy algorithms, GoDoTree outperforms the other two
methods due to its global optimality.

Table 3: Simulation results of scenario 3 and scenario 4, which are two-stage DTRs: n = 500, p =10, tree height = 3, $E\\{Y(g_{opt})\\}=0$; reported values are based on 100 Monte Carlo replications, using a test set of size 1000. Two-stage DTR | $100E\\{Y(\hat{g}_{opt})\\}$ | RMSE of $\hat{g}_{opt,1}$ | RMSE of $\hat{g}_{opt,2}$
---|---|---|---
Scenario 3 | GoDoTree | 1.1 (0.6) | 0.018 (0.008) | 0.006 (0.002)
LZ | 2.5 (1.2) | 0.051 (0.015) | 0.013 (0.013)
CART | 0.9 (0.6) | 0.013 (0.003) | 0.005 (0.001)
Scenario 4 | GoDoTree | 6.2 (1.6) | 0.019 (0.008) | 0.063 (0.015)
LZ | 11.3 (3.7) | 0.078 (0.044) | 0.115 (0.040)
CART | 8.9 (0.8) | 0.016 (0.007) | 0.090 (0.005)
Figure 3: Comparison of decision trees generated by the proposed method
GoDoTree with the truth and two other methods CART & LZ. The first two levels
of GoDoTree are nearly identical to the optimal decision rules, while there
are many sub-optimal nodes in CART and LZ decision trees.

Figure 3 highlights the limitations of greedy methods. In scenario 2, when
using a single informative variable $x_{1}$ or $x_{2}$, the greedy methods are
unable to select the optimal decision rule as the first splitting node, as
this does not result in a higher objective function at the early stage.
Consequently, the greedy trees may choose a sub-optimal decision rule as the
first node, which cannot be updated later, leading to an overly complicated
tree. In the CART model fitting process, there is a local optimum that is also
the global optimum when $x_{2}\textless-0.5$, and the decision rule
$x_{1}\textless 0.1\approx 0$ is nearly optimal. However, since the first node
$x_{2}\textless-0.5$ is sub-optimal, it is even more difficult for the second
level nodes to identify the correct pattern. As a result, CART selects a sub-
optimal rule $x_{3}\textless-0.5$ again. While increasing the tree height and
size may improve the probability of detecting the correct pattern, it also
increases the risk of overfitting and reduces interpretability.

##  7 Real Application for Optimal Warfain Dose Finding

Warfarin is a commonly used anticoagulant medication that requires precise
dosing to prevent harmful blood clots. Overdosing predisposes patients to a
high risk of bleeding, while underdosing diminishes the drug’s preemptive
protection against thrombosis. The international normalized ratio (INR) is
used to measure how rapidly the blood can clot and is monitored to ensure that
the dose of Warfarin is safe and effective. For patients prescribed Warfarin,
the optimal therapeutic INR range is typically between 2 to 3 (Krittayaphong
et al.,, 2020). To convert the INR to a direct measure of reward, the
literature codes reward $R=-100\times\lvert INR-2\rvert-100\times\lvert
INR-3\rvert$, which is a concave function and reaches maximum between 2 to 3
(Chen et al., 2016b, ).

The dataset provided by Klein et al., (2009) consists of 1780 subjects,
including information on patient covariates, final therapeutic warfarin
dosages, and patient outcomes (INR). The variables include weight, height,
age, use of certain medications, gender, race, and genetic information. The
optimal dose decision tree generated by GoDoTree is shown in Figure 4, and the
results agree with well-established medical knowledge in the
literature(Tatsuno and Tatsuno,, 2014). for example, patients with VKORC1
homozygous A/A and CYP2C9 alleles 2 or 3 requires a lower dosage of warfarin.
Furthermore, our estimated results suggest to offer optimal dosage
numerically, considering the possible interaction between genes and other
biomarkers (race, gender, height, and weight). Our decision tree also suggests
that Asian people should receive a higher dose, which may appear to contradict
common sense. However, this is because our analysis only considers the
criterion of “INR located between 2 to 3” as the optimal result, while in
practice physicians may consider additional factors that could result in an
optimal INR outside of this range.

Figure 4: Optimal Warfarin dose decision tree generated by GoDoTree, with
individual and population-level estimated effect curves shown for each leaf
node. Quantile information for weight and height is used due to de-identified
data. Effect curves are normalized to have a maximum value of zero since only
the optimal dose is of interest.

##  8 Discussion

The proposed GoDoTree is a novel tree-based learning approach for continuous
dosage finding in multiple decision stages. It falls under the category of
treatment-tree algorithms, which first estimate patient-specific effect curves
and then perform a supervised tree learning. This approach has two major
advantages: global optimality and great interpretability. The individualized
counterfactual outcome estimation and non-greedy tree search allow GoDoTree to
update sub-optimal nodes iteratively and recognize underlying patterns,
resulting in better global convergence. Additionally, the implementation of
the tree-based decision rule is straightforward and easy to interpret for
users. GoDoTree also has two ideal theoretical properties: doubly robustness,
which makes the results more stable, and asymptotic normality, which makes
statistical inference possible.

Due to greedy tree’s lack of optimality, some previous efforts have been made
to search for global optimal tree search, e.g. Hu et al., (2019) optimized
sparse decision tree with penalized objective function. But most of them are
supervised learning and cannot be combined with DTR learning easily because
individual counterfactual outcome is missing and only population level causal
effect can be estimated. There are also literature using stochastic search as
the remedy, e.g., Sun and Wang, (2021) is introducing random change of tree
structure with Markov Chain Monte Carlo (MCMC) to optimize the decision rule.
However, although it has improvements over random search, it still lacks
efficiency because of the nature of unsupervised learning. An important
contribution of our work is the development of individualized counterfactual
outcome estimation, which enables the global optimal tree learning of DTRs.
The node updating in tree learning is a simple bi-partition task with
supervision, which significantly speeds up the optimization training.

GoDoTree can be easily extended to scenarios with categorical treatment and
this work mainly focus on optimal dose finding because it has not been covered
thoroughly. The same approach also works when there are multiple objective
functions to be optimized as long as with a well defined utility function. The
optimal kernel search of GoDoTree can also be used to raise the information-
noise ratio based on the estimated variable importance, making it suitable for
sparse and high-dimensional data.

GoDoTree contributes to the development of precision health by enabling
precise dose finding for each individual and tailoring interventions to every
single patient. It can help clinicians search, validate, and refine new
decision rules, especially in dose finding of new drugs and radiation
oncology. The optimal dose can be quantified, rather than roughly estimated,
or categorized by experience.

## References

  * (1) Chen, G., Zeng, D., and Kosorok, M. R. (2016a).  Personalized dose finding using outcome weighted learning.  Journal of the American Statistical Association, 111(516):1509–1521.  PMID: 28255189.
  * (2) Chen, G., Zeng, D., and Kosorok, M. R. (2016b).  Personalized dose finding using outcome weighted learning.  Journal of the American Statistical Association.
  * Chipman et al., (2010) Chipman, H. A., George, E. I., and McCulloch, R. E. (2010).  Bart: Bayesian additive regression trees.  The Annals of Applied Statistics, 4(1).
  * Hayfield and Racine, (2008) Hayfield, T. and Racine, J. S. (2008).  Nonparametric econometrics: The np package.  Journal of Statistical Software, 27(5).
  * Hu et al., (2019) Hu, X., Rudin, C., and Seltzer, M. I. (2019).  Optimal sparse decision trees.  CoRR, abs/1904.12847.
  * Härdle et al., (1988) Härdle, W., Hall, P., and Marron, J. S. (1988).  How far are automatically chosen regression smoothing parameters from their optimum?  Journal of the American Statistical Association, 83(401):86–95.
  * Kennedy et al., (2016) Kennedy, E. H., Ma, Z., McHugh, M. D., and Small, D. S. (2016).  Non-parametric methods for doubly robust estimation of continuous treatment effects.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(4):1229–1245.
  * Klein et al., (2009) Klein, I. W. P. C. T. E., Altman, R. B., Eriksson, N., Gage, B. F., Kimmel, S. E., Lee, M.-T. M., Limdi, N. A., Page, D., Roden, D. M., Wagner, M. J., Caldwell, M. D., and Johnson, J. A. (2009).  Estimation of the warfarin dose with clinical and pharmacogenetic data.  New England Journal of Medicine, 360(8):753–764.  PMID: 19228618.
  * Krittayaphong et al., (2020) Krittayaphong, R., Kunjara-Na-Ayudhya, R., Ngamjanyaporn, P., Boonyaratavej, S., Komoltri, C., Yindeengam, A., Sritara, P., and Lip, G. Y. H. (2020).  Optimal inr level in elderly and non-elderly patients with atrial fibrillation receiving warfarin: a report from the cool-af nationwide registry in thailand.  Journal of geriatric cardiology.
  * Laber and Zhao, (2015) Laber, E. B. and Zhao, Y. Q. (2015).  Tree-based methods for individualized treatment regimes.  Biometrika, 102(3):501–514.
  * Moodie et al., (2012) Moodie, E. E. M., Chakraborty, B., and Kramer, M. S. (2012).  Q-learning for estimating optimal dynamic treatment rules from observational data.  The Canadian journal of statistics, 40(4):629–645.
  * Murphy, (2005) Murphy, S. A. (2005).  A generalization error for q-learning.  Journal of Machine Learning Research, 6(37):1073–1097.
  * Murphy et al., (2001) Murphy, S. A., van der Laan, M. J., Robins, J. M., and Group, C. P. P. R. (2001).  Marginal mean models for dynamic regimes.  96(456):1410–1423.
  * Murthy and Steven, (1995) Murthy, S. K. and Steven, S. (1995).  Decision tree induction: How effective is the greedy heuristic?  In KDD, pages 222–227. AAAI Press.
  * Perpinan and Tavallali, (2018) Perpinan, M. A. C. and Tavallali, P. (2018).  Alternating optimization of decision trees, with application to learning sparse oblique trees.  In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc.
  * Qian and Murphy, (2011) Qian, M. and Murphy, S. A. (2011).  Performance guarantees for individualized treatment rules.  The Annals of Statistics, 39(2):1180 – 1210.
  * Robins et al., (2000) Robins, J. M., Hernan, M. A., and Brumback, B. (2000).  Marginal structural models and causal inference in epidemiology.  Epidemiology, 11(5).
  * Robins and Rotnitzky, (2001) Robins, J. M. and Rotnitzky, A. (2001).  Comment on the bickel and kwon article, ’inference for semiparametric models: Some questions and an answer’.  Statistica Sinica, 11(4):920–936.
  * Rubin and van der Laan, (2005) Rubin, D. B. and van der Laan, M. J. (2005).  A general imputation methodology for nonparametric regression with censored data.
  * Sun and Wang, (2021) Sun, Y. and Wang, L. (2021).  Stochastic tree search for estimating optimal dynamic treatment regimes.  Journal of the American Statistical Association, 116(533):421–432.
  * Tao et al., (2018) Tao, Y., Wang, L., and Almirall, D. (2018).  Tree-based reinforcement learning for estimating optimal dynamic treatment regimes.  The Annals of Applied Statistics, 12(3):1914 – 1938.
  * Tatsuno and Tatsuno, (2014) Tatsuno, S. Y. and Tatsuno, E. M. (2014).  Does ethnicity play a role in the dosing of warfarin in hawaii?  Hawaii journal of medicine and public health: a journal of Asia Pacific Medicine and Public Health.
  * van der Laan and Robins, (2003) van der Laan, M. J. and Robins, J. M. (2003).  Unified Methods for Censored Longitudinal Data and Causality.  Book.
  * Zhao et al., (2012) Zhao, Y., Zeng, D., Rush, A. J., and Kosorok, M. R. (2012).  Estimating individualized treatment rules using outcome weighted learning.  Journal of the American Statistical Association, 107(499):1106–1118.  PMID: 23630406.


XX XX

Department: Head Editor: Name, xxxx@email

# Enabling Cost-Benefit Analysis of Data Sync Protocols

Novak Boškov  Ari Trachtenberg and David Starobinski  Boston University

(2023)

###### Abstract

The problem of data synchronization arises in networked applications that
require some measure of consistency. Indeed data synchronization approaches
have demonstrated a significant potential for improving performance in various
applications ranging from distributed ledgers to fog-enabled storage
offloading for IoT. Although several protocols for data sets synchronization
have been proposed over the years, there is currently no widespread utility
implementing them, unlike the popular Rsync utility available for file
synchronization. To that end, we describe a new middleware called _GenSync_
that abstracts the subtleties of the state-of-the-art data synchronization
protocols, allows users to choose protocols based on a comparative evaluation
under realistic system conditions, and seamlessly integrate protocols in
existing applications through a public API. We showcase _GenSync_ through a
case study, in which we integrate it into one of the world’s largest wireless
emulators and compare the performance of its included protocols.

††journal: IEEE Computer

Introduction Replication is a common thread among disparate distributed
systems, typically arising when there is a need for fault tolerance or
availability. Multiple replicas enable an administrator to reroute requests
from failed to healthy replicas, while a repair is completed. Repaired
replicas can then be returned to a consistent state from existing healthy
replicas.

Replication can also be used as a tool to improve the performance of
distributed systems with many concurrent accesses. For instance, globally
distributed databases and content delivery networks use replication to speed
up accesses from geographically dispersed locations.

Beyond availability, fault tolerance and performance benefits, replicated
systems can also facilitate decentralization. For example, in distributed
ledger technologies (_e.g.,_ blockchains), each participant holds a replica of
an entire data set (_i.e.,_ transactions), which allows any participant to
independently verify the state of system. As long as enough participants
maintain a substantial level of consistency among their replicas, no one needs
to trust a central authority.

At the core of replication is _data synchronization (sync)_ or
_reconciliation_ , the process that brings multiple replicas into a consistent
state. Assuming that there are two parties aiming at syncing their data sets,
a naive solution would simply exchange data sets. Having both sets, the
parties can easily identify the differing elements and include them in their
own replicas. However, this is prohibitively expensive, because communicating
huge data sets can be slow for bandwidth-constrained networks or resource
constrained compute platforms (_e.g.,_ augmented reality and wearables). Worse
yet, huge sets may differ in only a few elements, making most of the
communication redundant. The waste is trebled when multiple participants
exchange their sets in a peer-to-peer fashion, as is the case for distributed
ledgers.

### Better Ways to Synchronize Data

Many protocols have been developed over time to address the weaknesses of the
naive sync approach. Perhaps the most popular is Tridgell and Mackerras’s
Rsync [1], which is included in most Linux distributions and used as a default
network-enabled synchronization utility. Originally designed for files (or bit
arrays), Rsync operates by dividing the files into chunks, exchanging the
hashes of chunks, and generating _instructions_ on how to make two files
equal. These instructions are transferred over the network and applied to
synchronize the files.

Rsync does not readily generalize from files to data sets, however, there have
been other approaches proposed for that purpose. One such approach is somewhat
similar to Rsync and uses a specialized hash known as _Bloom filter_ to
compactly represent the data set. Rather than exchanging hashes of chunks like
Rsync, this approach exchanges Bloom filters. The main drawback of this
approach, however, is non-optimal usage of bandwidth — it exchanges some
traffic even for the set elements evident on both sides of the communication
channel. A _communication-efficient_ protocol, on the other hand, should only
exchange the information related to differences between the sets [2].

### Choosing the Right Sync Protocol

Several communication-efficient protocols have been proposed in the
literature. Some rely on coding theory, while other make use of probabilistic
data structures. These subtle but significant differences in design make it
important to be able to compare them under practical conditions. Yet prior to
our work [3], there were no publicly available general tools that afforded
such a systematic comparison. Our work has shown that a protocol may dominate
under certain network conditions, but grossly underperform when the network
conditions change. Worse yet, the best protocol choice is a function not only
of network conditions, but also of the immediate compute capacities of the
nodes.

### Contributions

There are two significant obstacles for developers in choosing the right
protocol for their applications:  (1) the lack of a utility for comparative
analysis of the sync protocols under practical system conditions, and  (2) the
complexity of integrating sync protocols into existing implementations.
Therefore, the main objective of this work is to overcome these obstacles and
enable the integration of the state-of-the art sync protocols in future
applications, such as 6G-enabled Enhanced Reality (ER), Internet-of-Things
(IoT) and Internet-of-Vehicles (IoV), where sync protocol customization is
needed at the application level. We summarize our contributions as follows:

  * •

We describe our middleware _GenSync_ , to the best of our knowledge, the first
utility that enables a systematic cost-benefit analysis of utilizing different
sync protocols.

  * •

We demonstrate the use of this middleware in an independent, large-scale
wireless emulator, _Colosseum_.

  * •

We show that the choice of sync protocols can significantly impair or improve
performance.

## Applications of Data Sync

Data sync serves as a building block for a diverse collection of applications
and computing paradigms, though it is sometimes embedded deeply within the
architecture of these applications. Next, we describe example applications
that can benefit from _GenSync_ , focusing on distributed ledgers, cloud
storage services, and IoT storage offloading.

### Distributed Ledgers

Distributed ledgers record the decentralized transactions of massive amounts
of participants. The most popular among distributed ledgers, blockchains,
implement their transaction inventory as a list of transaction blocks that are
logically chained together with cryptographic hashes.

Recent advances in blockchain and distribute ledger technologies have opened a
range of new possibilities for tackling long-standing problems in related
areas such as IoT access management [4], security of federated learning in fog
computing [5], accident forensics in vehicular networks (VANET) of self-
driving cars [6], or information poisoning prevention in mission-critical
unmanned aerial vehicle networks (UAANET) [7].

Applying blockchain technology over disparate layers of Distributed Computing
Continuum Systems (DCCS) [8], which are characterized by heterogeneous compute
and network resources, pose new challenges to performance and reliability. To
cope with these challenges, several improvements to the blockchain’s
networking layer have recently been proposed, based on set reconciliation
protocols, including MempoolSync [9], Graphene [10], and Erlay [11]. The
performance of these protocols has been shown to vary significantly depending
on the network conditions and compute capabilities of nodes [3], complicating
the analysis and choice of optimum reconciliation protocol and parameters.

### Cloud Storage Services

Cloud storage services such as Apple iCloud, Dropbox, Google Cloud, and
Microsoft OneDrive have also become commonplace for modern internet users and
are known to generate tremendous amounts of data traffic for the cloud
providers [12]. To reduce the amount of data transfer, these applications
employ protocols that are commonly referred to as _delta synchronization_. The
main objective of these protocols is to determine and transmit only those
portions of data that have been updated locally. In that vein, the research
community has defined a metric called TUE (Traffic Usage Efficiency) that is
the ratio between total sync data traffic and the update size [12]. The state-
of-the art delta synchronization protocols currently in use rely on improved
versions of Rsync.

### IoT Storage Offloading

One of the significant problems in IoT applications is storage offloading. A
standard IoT setup deals with a relatively large number of IoT devices that
each produce a substantial amount of data but lack storage and compute
capacity to maintain it. Traditionally, this problem has been addressed by
moving data to the cloud for processing. However, there are two significant
drawbacks of such approaches:  (1) the data is physically transferred to
another entity, which raises various security concerns, and  (2) data access
latency may be unacceptable for real-time applications as each data access
goes through a wide-area network.  To tackle these drawbacks, Wang  _at al._
[13] have come up with a fog-based architecture, based on Rsync, that allows
for data storage within the boundaries of the entity that operates the IoT
devices. The proposed architecture has three layers. The lowest layer consists
of IoT devices that synchronize their data with the fog layer above them. The
fog layer accumulates this data and synchronizes it with the cloud in batches.

### Other Applications

There are many other distributed systems where data sync is used as a building
block. For instance, researchers proposed a dissemination protocol for
wireless sensor networks that use a variant of Bloom filter-based data sync to
reduce energy consumption and propagation delay [14]. Similarly, a
physiological value-based key agreement scheme for body area networks called
E2PKA [15] uses data sync to reduce memory footprint and energy consumption.
Data sync protocols have also been proposed as a solution to network
partitioning in information-centric networking (ICN) [16] as well as re-
establishing consistency among replicas of distributed databases [17]. As
Distributed Computing Continuum Systems [8, 18] and applications that operate
across the layers of this continuum emerge, we expect that many more
applications could take advantage of an efficient data sync middleware.

##  _GenSync_ Middleware

Figure 1: Middleware design.

We implement _GenSync_ framework as a C++ middleware library that can be used
through an API, as shown in Fig. 1. The main challenges that we tackled in
designing _GenSync_ were  (1) unifying and simplifying the parameterization of
disparate state-of-the-art data sync protocols;  (2) designing generic
protocol implementations that enable _GenSync_ utilization in various
applications; and  (3) constructing a lightweight yet versatile benchmarking
layer.  From the _GenSync_ users’ standpoint, the interaction between the
application and the middleware library is carried out through three
abstractions, named _GenSync_ , _Observation_ , and _Builder_.

GenSync is a generic representation of a data sync protocol. Its interface
consists of two main methods, _addElement_ , which adds elements to the
associated data set, and _syncBegin_ , which initiates synchronization with an
external party using the desired protocol. In a typical application, set
elements can be easily mapped to a set of unique identifiers using a hash
function and passed to _addElement_. We design _GenSync_ as an abstract
interface to allow for middleware extensions, which is particularly useful to
researchers that design novel set reconciliation protocols and want to
benchmark against the state-of-the art, and the practitioners that design
platform-specific implementations of existing protocols. A custom sync
implementation needs only to implement the two _GenSync_ methods. By
implementing _addElement_ , users can control how data set is being
maintained, while _syncBegin_ can be used to provide the core implementation
of the custom protocol.

Observation is the summarized result of the sync and will be generated by
_syncBegin_. It represents a collection of execution statistics including the
measure of success, exact protocol parameters that have been used and
monitoring information, such as the communication cost and time expended on
data transfer and computation.

Builder is an auxiliary abstraction that facilitates the creation and
connection of _GenSyncs_ to each other. For example, an application may want
to instantiate multiple _GenSync_ objects (say, one for each of several
neighbors in a peer-to-peer system). The builder allows the application to
attach each _GenSync_ object to remote peers through _Communicants_ , which
abstract out a communication channel (_e.g.,_ a TCP, UDP, or local Unix
socket). The _Builder_ abstraction is particularly useful in heterogeneous
environments, where peers may connect using different underlying transport or
even physical-layer protocols, or may want to use individually optimized sync
protocols for different neighbors. Code listing 2 shows how _Builder_ and
_GenSync_ can be chained together to conduct the sync and produce an
_Observation_.

#include <GenSync.h>// Point to a remote and pick a protocolauto builder =
GenSync::Builder();builder.setProtocol(GenSync::CPI);builder.setCommunicant(GenSync::socket);builder.setHost("the.peer.remote.addr");GenSync
gs = builder.build();// Add datafor (auto data_point : data_set)
gs.addElement(hash(data_point));// Perform syncif ( gs.syncBegin() ) { // Get
execution statistics Observation ob = gs.getObservation();
ob.communicationTime; ob.computationTime; ob.bytesTransmitted;} else { // Sync
failed}

Figure 2: Illustrative usage of _GenSync_ from an application.

### Benchmarking Layer

To allow for performance evaluation under realistic system conditions,
_GenSync_ is equipped with a _benchmarking_ layer. The benchmarking layer
creates an execution environment, based on the cgroups feature of the Linux
kernel, that simulates the target system with a given set of system
parameters. Developers can readily extract the _Observations_ from this
simulated environment for further analysis, as we expose _GenSync_ ’s
benchmarking layer through a script wherein developers can configure the
desired system conditions and the sync protocol to evaluate (see listing 3).

# Protocol identifierprotocol=CPI# Latency in millisecondslatency=20#
Bandwidth in Mbps (in two directions)bandwidth="10/25"# Packet loss
(percentage)packet_loss=0.01# Percentage of CPU cycles used for
synccpu_server=100cpu_client=20# Repeat each experimentrepeat=100

Figure 3: _GenSync_ ’s benchmarking layer configuration script.

## Included Protocols

_GenSync_ includes a number of sync protocols that are based on compact
auxiliary data structures (_sketches_) through a similar high-level structure.
Roughly speaking, the structure of these protocols consists of four phases: 1
compute sketches of local data, 2 exchange the sketches between syncing hosts,
3 compute local differences, and 4 exchange differences. The protocols
themselves are distinguished through their choice of sketches and how the
sketches are utilized to infer the differences between the sets. Specifically,
users can currently select from the following protocols.

CPI (Characteristic Polynomial Interpolation) sync is based on a representing
data sets as characteristic polynomials. In phase 1, both parties encode their
elements as zeros of a characteristic polynomial; they exchange evaluations of
these polynomials in 2. In 3, one party extrapolates the rational function
resulting from dividing these polynomials, and extracts the roots of this
function to determine the set differences, which are then exchanged in 4.

Cuckoo sync uses Cuckoo filters as sketches. In 1, parties insert all their
elements into a Cuckoo filter, and exchange them in stage 2. In 3, each party
queries its own data set elements against the Cuckoo filter of the other
party. Any element for which the Cuckoo filter returns a negative answer is
certainly only locally available and should thus be sent over in 4.

IBLT (Invertible Bloom Lookup Table) sync uses IBLTs as sketches, which makes
it somewhat similar to Cuckoo. In 1, each party constructs their own IBLT and
exchange them in 2. In 3, one party can subtract the other party’s IBLT from
their own to learn the elements that it needs to exchange in 4.

Figure 4: Logarithm of sync time as a function of set cardinality (size) and
the number of differences. Figure 5: Protocol transfer size as a function of
the number of mutual differences. Data set size is constant at 10 000. Figure
6: Protocol transfer size as a function of the data set size. The number of
mutual differences is constant at 100.

### Navigating Trade-offs

Typical performance metrics for data sync protocols are _transferred data
size_ (sum of all traffic until the sets are in sync) and _total sync time_.
The transferred data size depends on _data set parameters_ (_i.e.,_ size of
the sets and the number of their mutual differences) and the protocols’
theoretical bounds on communication complexity. The total sync time, however,
depends on _system parameters_ , which includes network bandwidth and latency
(jointly referred to as _network conditions_), and the compute capabilities of
the nodes (_compute conditions_). Given the systems’ complexity, the total
sync time cannot easily be estimated using only the theoretical bounds. Worse
yet, a bad protocol choice for the given system parameters can cause a 5x loss
in the total sync time performance [3].

Using the _GenSync_ ’s benchmarking layer, we modeled a bandwidth-constrained
system to explore the effects of data set parameters on total sync time. We
varied set size from ten thousand to one hundred thousand, and the number of
differences between zero and 300\. We ran _GenSync_ ’s benchmarking layer on
an Intel Core i7-7700 experimental server with 5.18.10 version of the Linux
kernel. The benchmarking layer parameters appear in listing 3. The resulting
surfaces are plotted in Fig. 4, where we can observe the following trends.

The sync time of Cuckoo is largely _invariant_ to the number of differences.
Trend 1

In other words, Cuckoo sync performs relatively well for very small sets, but
worsens as the set size increases (regardless of the differences count). This
can be explained with two observations:  (1) we are dealing with a bandwidth-
constrained network, and  (2) the transfer size for Cuckoo increases (in
steps) with the set size (Fig. 6) but stays almost constant as a function of
differences count (Fig. 5).  The slight increase in transfer size that we
observe in Fig. 5 is mostly due to the final transfer of the differences
themselves.

The sync times of IBLT and CPI are largely _invariant_ to the size of the data
sets being synced. Trend 2

Figure 7: Bandwidth (up) and latency (down) traces from Colosseum [19], one of
the world’s largest wireless network emulators.

That is, IBLT and CPI perform well relative to Cuckoo for very similar sets
(_i.e.,_ differing in only a few elements). IBLT and CPI generally do _not_
transfer data for the elements that are common. Moreover, CPI transfers a
nearly optimal amount of data per difference [2], whereas IBLT transfers more
(Fig. 6). Whether this discrepancy in transferred data size will result in
CPI’s dominance (with respect to total sync time) is the matter of system
parameters, and can be evaluated through _GenSync_ ’s benchmarking layer.

## Sync on the Edge

To showcase the _GenSync_ middleware’s power in estimating the actual sync
performance in practical systems, we apply it to Colosseum, one of the world’s
largest wireless network emulators [19], capable of emulating various real
world radio frequency _scenarios_ using software-defined radio (SDR)
technology. We use Colosseum’s “Boston” cellular network scenario to emulate
the cellular network in the vicinity of Boston Common in Boston,
Massachusetts. The scenario has _stationary_ and _pedestrian_ regimes, where
the latter captures moderate user movement relative to the base station (see
Table 1).

Scenario Regime |  User Speed |  Base Station Distance |  Scenario Duration
---|---|---|---
Stationary |  0 m/s |  20 m |  600 s
Pedestrian |  5 m/s |  20 m |  600 s
Table 1: Emulation parameters for the Colosseum’s Boston scenario [19].

The resulting network traces (the available bandwidth and latency) are plotted
in Fig. 7, with their extremes annotated. These traces show that user movement
results in wider oscillations of available bandwidth and latency. Using the
average values of bandwidth and latency during the extreme periods, we define
two sets of network conditions against which to evaluate our sync protocols:

  1. (1)

_bad_ (bandwidth 1 Mbps, latency 50 ms), and

  2. (2)

_good_ (bandwidth 7 Mbps, latency 30 ms).

The resulting total sync time for the two sets of network conditions and the
three _GenSync_ protocols is plotted in Fig. 8. IBLT performs the best in both
bad and good network conditions. The reason for this effect lies in IBLT’s
balance of low computational complexity (linear in the number of differences)
and communication cost that, although higher than CPI[2], beats Cuckoo for
large data sets. As the average bandwidth in both good and bad network
conditions does not drop below a critical point, IBLT also beats CPI.

Figure 8: Sync time for the three main sync protocols in _good_ (solid line)
and _bad_ (dashed line) network conditions. Data cardinality is $10^{8}$.
Confidence intervals are shaded. Figure 9: Bandwidth consumed for the two most
bandwidth-efficient sync protocols that still complete under 2 seconds in
_bad_ network conditions.

However, an application-specific protocol that is being constructed using
_GenSync_ may have additional performance objectives. For instance, it may
want to be conservative about the amount of _consumed bandwidth_ , while still
keeping a reasonable sync time (even under bad network conditions). Since the
synchronizing device may run several applications concurrently, this kind of
bandwidth budgeting could be an important dimension to consider when designing
an application-specific sync protocol. Bandwidth savings in the sync protocol
could generate substantial performance gains for Quality-of-Experience-
critical processes, such as streaming point clouds in augmented reality (AR)
applications [20]. In Fig. 9 we plot the amounts of bandwidth that IBLT and
CPI consume when sync time is constrained to two seconds. In this case, CPI
achieves almost nine times better performance across all difference counts.
The reason for CPI’s dominance over IBLT in this scenario is CPI’s nearly
optional communication cost, while IBLT adds a multiplicative constant to that
cost.

## Conclusions and Future Goals

The _GenSync_ middleware library is the first open and general framework that
(1) enables comparative evaluation of state-of-the-art data synchronization
protocols in practical environments through a versatile benchmarking layer,
and  (2) allows developers to seamlessly integrate the protocol of their
choice into their applications.

A current limitation of _GenSync_ is that is non-adaptive, in the sense that
it cannot intelligently detect the sways in system conditions (_e.g.,_
available bandwidth) and automatically replace the current protocol with
better one for the new conditions. We leave this promising direction for
future work.

We also suggest exploring the extension of _GenSync_ ’s protocols to file (bit
array) sync, with the hope of improving the venerable Rsync-based techniques.

## Acknowledgments

The authors would like to thank Sevval Simsek for discussions regarding the
Colosseum framework. The authors would also like to thank Red Hat, the Boston
University Red Hat Collaboratory (award # 2022-01-RH03), and the US National
Science Foundation (award # CNS-2210029) for their support.

## References

  * [1] A. Tridgell, “Efficient Algorithms for Sorting and Synchronization,” Ph.D. dissertation, The Australian National University, 1999.
  * [2] Y. Minsky, A. Trachtenberg, and R. Zippel, “Set Reconciliation with Nearly Optimal Communication Complexity,” _IEEE Transactions on Information Theory_ , vol. 49, no. 9, pp. 2213–2218, 2003.
  * [3] N. Boškov, A. Trachtenberg, and D. Starobinski, “GenSync: A New Framework for Benchmarking and Optimizing Reconciliation of Data,” _IEEE Transactions on Network and Service Management_ , pp. 1–1, 2022.
  * [4] O. Novo, “Blockchain Meets IoT: An Architecture for Scalable Access Management in IoT,” _IEEE Internet of Things Journal_ , vol. 5, no. 2, pp. 1184–1195, 2018.
  * [5] W. Issa, N. Moustafa, B. Turnbull, N. Sohrabi, and Z. Tari, “Blockchain-Based Federated Learning for Securing Internet of Things: A Comprehensive Survey,” _ACM Computing Surveys_ , August 2022.
  * [6] M. Cebe, E. Erdin, K. Akkaya, H. Aksu, and S. Uluagac, “Block4Forensic: An Integrated Lightweight Blockchain Framework for Forensics Applications of Connected Vehicles,” _IEEE Communications Magazine_ , vol. 56, no. 10, pp. 50–57, 2018.
  * [7] K. Lei, Q. Zhang, J. Lou, B. Bai, and K. Xu, “Securing ICN-Based UAV Ad Hoc Networks with Blockchain,” _IEEE Communications Magazine_ , vol. 57, no. 6, pp. 26–32, 2019.
  * [8] S. Dustdar, V. Casamajor Pujol, and P. K. Donta, “On distributed computing continuum systems,” _IEEE Transactions on Knowledge and Data Engineering_ , pp. 1–1, 2022.
  * [9] M. A. Imtiaz, D. Starobinski, A. Trachtenberg, and N. Younis, “Churn in the Bitcoin Network,” _IEEE Transactions on Network and Service Management_ , pp. 1–1, 2021.
  * [10] A. P. Ozisik, G. Andresen, B. N. Levine, D. Tapp, G. Bissias, and S. Katkuri, “Graphene: efficient interactive set reconciliation applied to blockchain propagation,” in _Proceedings of the ACM Special Interest Group on Data Communication_ , 2019, pp. 303–317.
  * [11] G. Naumenko, G. Maxwell, P. Wuille, A. Fedorova, and I. Beschastnikh, “Erlay: Efficient transaction relay for bitcoin,” in _Proceedings of the 2019 ACM SIGSAC Conference on Computer and Communications Security_ , 2019, pp. 817–831.
  * [12] Z. Li, C. Jin, T. Xu, C. Wilson, Y. Liu, L. Cheng, Y. Liu, Y. Dai, and Z.-L. Zhang, “Towards Network-Level Efficiency for Cloud Storage Services,” in _Proceedings of the Internet Measurement Conference (IMC)_. New York, NY, USA: Association for Computing Machinery, 2014, p. 115–128.
  * [13] T. Wang, J. Zhou, A. Liu, M. Z. A. Bhuiyan, G. Wang, and W. Jia, “Fog-based computing and storage offloading for data synchronization in IoT,” _IEEE Internet of Things Journal_ , vol. 6, no. 3, pp. 4272–4282, 2018.
  * [14] T. Chen, D. Guo, X. Liu, H. Chen, X. Luo, and J. Liu, “Bdp: A bloom filters based dissemination protocol in wireless sensor networks,” in _2009 IEEE 6th International Conference on Mobile Adhoc and Sensor Systems_ , 2009, pp. 593–602.
  * [15] W. Choi, I. S. Kim, and D. H. Lee, “E2PKA: An Energy-Efficient and PV-Based Key Agreement Scheme for Body Area Networks,” _Wireless Personal Communications_ , vol. 97, no. 1, pp. 977–998, Nov 2017.
  * [16] Z. Zhu and A. Afanasyev, “Let’s ChronoSync: Decentralized dataset state synchronization in Named Data Networking,” in _2013 21st IEEE International Conference on Network Protocols (ICNP)_ , 2013, pp. 1–10.
  * [17] D. Chen, C. Konrad, K. Yi, W. Yu, and Q. Zhang, “Robust Set Reconciliation,” in _Proceedings of the ACM International Conference on Management of Data (SIGMOD)_ , 2014, pp. 135–146.
  * [18] P. K. Donta and S. Dustdar, “The Promising Role of Representation Learning for Distributed Computing Continuum Systems,” in _2022 IEEE International Conference on Service-Oriented System Engineering (SOSE)_ , 2022, pp. 126–132.
  * [19] L. Bonati, P. Johari, M. Polese, S. D’Oro, S. Mohanti, M. Tehrani-Moayyed, D. Villa, S. Shrivastava, C. Tassie, K. Yoder, A. Bagga, P. Patel, V. Petkov, M. Seltser, F. Restuccia, A. Gosain, K. R. Chowdhury, S. Basagni, and T. Melodia, “Colosseum: Large-Scale Wireless Experimentation Through Hardware-in-the-Loop Network Emulation,” in _Proc. of IEEE International Symposium on Dynamic Spectrum Access Networks (DySPAN)_ , December 2021.
  * [20] L. Wang, C. Li, W. Dai, J. Zou, and H. Xiong, “QoE-Driven and Tile-Based Adaptive Streaming for Point Clouds,” in _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2021, pp. 1930–1934.

Novak Boškov is a Ph.D. candidate at Boston University’s Department of
Electrical and Computer Engineering. His research interests include
decentralized networks, cloud computing, and cybersecurity. Contact him at
<EMAIL_ADDRESS>
Ari Trachtenberg is a professor at Boston University, with appointments in
Electrical and Computer Engineering, Computer Science, and Systems
Engineering. His research interests include cybersecurity, distributed
systems, and information theory. Contact him at<EMAIL_ADDRESS>
David Starobinski is a professor at Boston University, with appointments in
Electrical and Computer Engineering, Computer Science, and Systems
Engineering. His research interests include cybersecurity, wireless
networking, and network economics. Contact him at<EMAIL_ADDRESS>