# Limits of Approximating the Median Treatment Effect111The author ordering is
alphabetical, as is customary in the theoretical computer science community.

Raghavendra Addanki
Adobe Research
<EMAIL_ADDRESS>Siddharth Bhandari
Toyota Technological Institute at Chicago
<EMAIL_ADDRESS>
###### Abstract

Average Treatment Effect (ATE) estimation is a well-studied problem in causal
inference. However, it does not necessarily capture the heterogeneity in the
data, and several approaches have been proposed to tackle the issue, including
estimating the Quantile Treatment Effects. In the finite population setting
containing $n$ individuals, with treatment and control values denoted by the
potential outcome vectors $\mathbf{a},\mathbf{b}$, much of the prior work
focused on estimating median$(\mathbf{a})-$ median$(\mathbf{b})$, where
median($\mathbf{x}$) denotes the median value in the sorted ordering of all
the values in vector $\mathbf{x}$. It is known that estimating the difference
of medians is easier than the desired estimand of median$(\mathbf{a-b})$,
called the Median Treatment Effect (MTE). The fundamental problem of causal
inference – for every individual $i$, we can only observe one of the potential
outcome values, i.e., either the value $a_{i}$ or $b_{i}$, but not both, makes
estimating MTE particularly challenging. In this work, we argue that MTE is
not estimable and detail a novel notion of approximation that relies on the
sorted order of the values in $\mathbf{a-b}$. Next, we identify a quantity
called _variability_ that exactly captures the complexity of MTE estimation.
By drawing connections to the notions of instance-optimality studied in
theoretical computer science, we show that _every_ algorithm for estimating
the MTE obtains an approximation error that is no better than the error of an
algorithm that computes variability. Finally, we provide a simple linear time
algorithm for computing the variability exactly. Unlike much prior work, a
particular highlight of our work is that we make no assumptions about how the
potential outcome vectors are generated or how they are correlated, except
that the potential outcome values are $k$-ary, i.e., take one of $k$ discrete
values.

##  1 Introduction

In the realm of experimentation, much attention is dedicated to constructing
efficient estimators for the average treatment effect (ATE), where machine
learning-based double robust estimation techniques are considered state-of-
the-art [Chernozhukov et al., 2018; Kennedy et al., 2017]. However, making
policy or downstream task decisions only based on ATE, might be sub-optimal,
as there is a loss of valuable information from the heterogenous regions in
the data. There have been many methods proposed in the literature for tackling
this issue [Imai and Ratkovic, 2013; Athey and Imbens, 2016; Ding et al.,
2019], and we focus on investigating variations across different outcome
distribution quantiles. Unlike ATE, which averages out the responses, Quantile
Treatment Effects (QTE) [Chernozhukov and Hansen, 2005] provide a more nuanced
perspective by unveiling whether effects differ among individuals or outcomes,
allowing for a more detailed understanding. This approach proves invaluable in
discerning if a treatment is more impactful at specific points in the
distribution, facilitating the study of policies affecting diverse population
segments, and groups unfairly affected. Several cautionary tales exist in the
literature [Heckman et al., 1997], with recent examples related to Netflix’s
A/B testing oversight of treatment effects disparities among user engagement
levels [Bojinov et al., 2020].

We study the finite populating setting on $n$ individuals, where an
experimenter assigns the individuals to treatment and control groups. Suppose
that the treatment and control values of all the individuals be denoted by
$\mathbf{a}$ and $\mathbf{b}$ respectively. Much of the work on QTE estimation
has focused on estimating the difference of quantiles of treatment and control
[Howard and Ramdas, 2022; Martinez-Taboada and Kennedy, 2023]. E.g., for
median estimation, the estimand studied is
$\textit{median}(\mathbf{a})-\textit{median}(\mathbf{b})$. Although efficient
estimators are known, such an estimand has a limitation, as it is not based on
the causal effect vector $\mathbf{a}-\mathbf{b}$, which of course is not
observable. It might not necessarily capture the causal effect entirely. In
contrast, some prior work considered the estimand
$\textit{median}(\mathbf{a}-\mathbf{b})$, broadly speaking [Heckman et al.,
1997; Leqi and Kennedy, 2021; Kallus, 2022; Kennedy, 2023]. Recently, Kallus
[2022] studied a related problem of identifying the fraction of individuals
that are negatively affected (under binary potential outcomes). They provide
efficient estimators, but assume the availability of covariates and make
distributional assumptions.

In this work, we study the estimation of
$\textit{median}(\mathbf{a}-\mathbf{b})$, henceforth called the Median
Treatment Effect (MTE) estimation, without any assumptions, except that the
potential outcomes, i.e., the vectors $\mathbf{a}$ and $\mathbf{b}$ are binary
(or $k$-ary). To the best of our knowledge, ours is the first work that
tackles this problem in a general setting with minimal or no assumptions. We
want to highlight that we do not make any distributional assumptions, and the
vectors $\mathbf{a}$ and $\mathbf{b}$ can be arbitrarily correlated. Our
results contribute to the growing body of work that identifies connections
between theoretical computer science and causal inference [Pouget-Abadie et
al., 2019; Addanki et al., 2022; Harshaw et al., 2023]. In our work, we answer
the following questions.

  1. 1.

Is Median Estimable? We show that MTE estimation is infeasible, in stark
contrast to ATE estimation. In particular, we show that _any_ estimator
(algorithm) is unable to distinguish between two distributions from which the
treatment and control values are drawn if they have the same marginal
distributions. Further, the two distributions can have a large difference in
their median values.

  2. 2.

Can we approximate MTE? To overcome the difficulty, we define a notion of
approximation that does not rely on the estimate being additively close to the
median but instead approximates the median by being close to in terms of
quantiles, i.e., the rank in the sorted order of the values in $\mathbf{a-b}$.
Further, we show that we can approximately estimate MTE under this notion.

  3. 3.

Limits of Approximating MTE. We identify a quantity called _variability_ ,
that essentially captures the complexity of MTE estimation. We show that any
algorithm that computes variability obtains the minimum approximation error,
and also provides an efficient greedy algorithm for computing the variability.

###  1.1 Technical Overview

Throughout the overview, we focus on binary outcomes (i.e., $k=2$), where
$\mathbf{a},\mathbf{b}$ are binary. All our results extend to the situation
where we have $k$-ary outcomes. When $k=2$, the MTE which is the median of the
difference vector $\mathbf{a-b}$ can be only one of the following values:
$\\{-1,0,1\\}$. An experimenter provides us with an assignment while ensuring
that every individual $i$ is assigned to either the treatment or control group
but not both. Hence, we observe either $a_{i}$ or $b_{i}$ but not both.
Further, the design may be adaptive or non-adaptive, i.e., the experimenter
can choose to assign an individual based on all the previously assigned
individuals and their responses. We also use $[n]$ to denote the set
$\\{1,2,\ldots,n\\}$ and for a distribution $\eta$ we use $\eta^{\otimes n}$
to denote the $n$-fold product distribution.

####  1.1.1 Is Median Estimable?

Consider two distributions $\mu_{1}$: $\begin{bmatrix}1/3&0\\\
1/3&1/3\end{bmatrix}\mbox{ and }\mu_{2}:\begin{bmatrix}0&1/3\\\
2/3&0\end{bmatrix}$, that encode the joint distributions of the pair
$(\mathbf{a},\mathbf{b})$. So, essentially, we draw $n$ samples from $\mu_{1}$
or $\mu_{2}$ to construct $(\mathbf{a},\mathbf{b})$. Here, the first row
encodes the probabilities: For an $i$th individual drawn from $\mu_{1}$:
$\Pr_{\mu_{1}}[(a_{i},b_{i})=(0,0)]=\mu_{1}[0,0]=1/3$,
$\Pr_{\mu_{1}}[(a_{i},b_{i})=(0,1)]=\mu_{1}[0,1]=0$, and so on.

Observe that for both $\mu_{1}$ and $\mu_{2}$, the marginal distributions are
the same, i.e., $\Pr_{\mu_{1}}[a_{i}=0]=\Pr_{\mu_{2}}[a_{i}=0]=1/3$ and
$\Pr_{\mu_{1}}[b_{i}=0]=\Pr_{\mu_{2}}[b_{i}=0]=2/3$. Therefore, the ATE of
$\mathbf{a-b}$ for both the distributions is the same and equals $1/3$.
However, the MTE for samples drawn from $\mu_{1}$ is $0$, and for $\mu_{2}$ is
$1$. This is because the total fraction of $1s$ in a typical sample from
$\mu_{2}$ is $2/3$ and hence contains the median, which is the middle element
in the sorted order. However, in $\mu_{1}$, the fraction of $-1$s is $0$ and
the fraction of $0$s is $\mu_{1}[0,0]+\mu_{1}[1,1]=2/3$, which exceeds the
$1/2$ threshold, therefore, the MTE for $\mu_{1}$ is $0$. This example
illustrates the inherent difficulty in estimating the median from the observed
responses in the population. We formalize the notion and provide a brief
intuition.

We can think of any algorithm for estimating the MTE as a rooted and labeled
decision tree of depth $n$. Each non-leaf node $u$, is labeled with a
distribution $\gamma_{u}$ defined over $\\{a,b\\}\times[n]$, where $a$
indicates treatment and $b$ indicates control. For each element say $(c,i)$ in
the support there are two edges out of $u$ labelled $(c,i,0)$ and $(c,i,1)$
respectively. We start from the root, and whenever we are at a non-leaf node
$u$, our algorithm samples $(c,i)$ according to $\gamma_{u}$ and assigns the
$i$th individual to treatment if $c=a$ and control if $c=b$. Suppose we see a
response value of $0$ in individual $i$, then, we follow the edge $(c,i,0)$;
if $i$th we observe a response value $1$, we follow the edge $(c,i,1)$.
Proceeding in this way we trace a path, say $P$, from the root to a leaf of
the decision tree. The feasibility constraint is that along any path $P$ that
has a non-zero probability of being followed, any individual $i\in[n]$ is
assigned exactly once. We define such decision trees as _Randomized Feasible
Decision Tree_ (RFDTs) (see Definitions 2.11 and 2.13 for a precise
definition). At the end of following a path $P$, we output an estimate for the
MTE say $\widehat{m}\in\\{-1,0,1\\}$ (we can think of the leaf node of $P$ as
being labeled with the output). Notice that our definition of RFDTs is general
enough to capture non-adaptive designs. In fact, non-adaptive designs are
those where the distribution $\gamma_{u}$ depends only on the depth of $u$ in
the decision tree.

A crucial observation is that any RFDT, $R_{n}$ is incapable of distinguishing
between two joint distributions over $\\{0,1\\}\times\\{0,1\\}$, say $\eta$
and $\tau$, which have the same marginals, i.e., if we sample $(a,b)\sim\eta$
and call the marginal of $a$ as $\eta_{a}$ and of $b$ as $\eta_{b}$ (similarly
for $\tau_{a}$ and $\tau_{b}$), then, $\eta_{a}=\tau_{a}$ and
$\eta_{b}=\tau_{b}$. Specifically, any RFDT $R_{n}$ is unable to distinguish
whether the input is drawn from the distribution $\eta^{\otimes n}$ or
$\tau^{\otimes n}$. The reason is that for any path $P$ of $R_{n}$ the
probability that $P$ is followed on the input $(\mathbf{a,b})$ drawn from
$\eta^{\otimes n}$ is a function of the distributions associated with the
nodes of $P$, and the marginal distributions $\eta_{a}\text{ and }\eta_{b}$.
Since, this probability only depends on the marginals $\eta_{a}=\tau_{a}$ and
$\eta_{b}=\tau_{b}$, hence the probability of following $P$ on $\tau^{\otimes
n}$ remains the same (see Lemma 3.1 for a detailed proof)

####  1.1.2 Approximating the Median

As we saw previously, for two distributions $\eta$ and $\tau$, any algorithm
(RFDT) is incapable of differentiating between $\eta^{\otimes n}$ and
$\tau^{\otimes n}$. We can only (approximately) estimate the marginals from
the responses, and thus any algorithm that estimates MTE, must in principle
use only the marginals. However, the MTEs of treatment and control vectors
drawn from $\eta^{\otimes n}$ and $\tau^{\otimes n}$ can be quite different,
as shown in the example taken previously. Therefore, we realize that an
unbiased estimation of the MTE is impossible. Notice the contrast with ATE
estimation, as the ATE value is the same whether the data is generated
according to $\eta^{\otimes n}$ and $\tau^{\otimes n}$. This is because the
ATE is linear in the outcomes $\mathbf{a},\mathbf{b}$ but the MTE is not.

Faced with the above roadblock in estimating the MTE, a natural question to
ask if we can obtain an additive approximation to the MTE, i.e., output an
$m\in\mathbb{R}$ such that the MTE is within $\pm\epsilon$ of the MTE for an
approximation parameter $\epsilon>0$. However, for the situation of $k=2$ (or
other values of $k$) this is not a desirable notion as we already know that
the MTE is one of $\\{-1,0,1\\}$, so either we need to set $\epsilon=1$, which
is too large, or know the MTE exactly.

To overcome this, we introduce a notion of approximation that approximates the
median by being close to in terms of quantiles, i.e., the rank in the sorted
order. More precisely, we define $\widehat{m}$ to be an $\epsilon$-_quantile-
approximation_ of the median if it lies in the $0.5\pm\epsilon$-quantile band
around the median, denoted by $(\mathbf{a-b})^{0.5+\epsilon}_{0.5-\epsilon}$.
An $\epsilon$-quantile band around the median is the set of $\epsilon\cdot n$
values of the vector $\mathbf{a-b}$ in its sorted order, above and below the
median.

####  1.1.3 Optimality of Median Approximation

To accommodate the above notion of approximation we add a parameter to the
output of an algorithm (RFDT) for estimating the MTE. In addition to an
estimate $\widehat{m}$ of the MTE, we also output the _width_ parameter
$\epsilon$, with our guess being that
$\widehat{m}\in(\mathbf{a-b})^{0.5+\epsilon}_{0.5-\epsilon}$.

How should we measure the performance of any algorithm (RFDT), say $R_{n}$
($n$ denotes the number of individuals), that approximately estimates MTE? We
measure it via two important parameters: how often is $R_{n}$ correct, i.e.,
$m$ is indeed an $\epsilon$-quantile approximation to the MTE, and what is the
expected width $R_{n}$ outputs. More precisely, for some small $\delta>0$ we
desire that for _all inputs_ $(\mathbf{a,b})$ (there is no randomness in the
inputs) with probability at least $1-\delta$ over the internal randomness of
$R_{n}$ we have $\Pr[m\in(\mathbf{a-b})^{0.5+\epsilon}_{0.5-\epsilon}]$ where
$(m,\epsilon)$ denotes the random output of $R_{n}$ on the input
$(\mathbf{a,b})$. The width of $R_{n}$ on the input $(\mathbf{a,b})$, denoted
by $\epsilon(\mathbf{a,b})$ is $\operatorname*{{\mathbb{E}}}[\epsilon]$, where
the expectation is over the randomness of the algorithm.

Naturally, we desire to construct an RFDT $R_{n}$ such that it has a small
error ($\delta$) and a small width for any input $(\mathbf{a,b})$. But how
small can we make the width? Recall the input distributions $\mu_{1},\mu_{2}$
discussed in Section 1.1.1. We know that $R_{n}$ can’t distinguish between
$\mu_{1}^{\otimes n}$ and $\mu_{2}^{\otimes n}$. Therefore, if upon observing
the responses, $R_{n}$ decides to output $\widehat{m}=0$, then to be correct
with high probability it must output the width $\epsilon$ as at least $1/6$:
this is because for a typical sample $\mathbf{a,b}$ from $\mu_{2}^{\otimes n}$
we have that the rank of $0$ in the sorted version of $\mathbf{a-b}$ is $n/3$.
Similarly, we observe that the minimum widths required to be output by $R_{n}$
for each of the estimates $-1,0,\text{ and }1$ are $1/2,1/6,\text{ and }1/6$
respectively. Hence, $(1/6)^{th}$ is in principle a lower bound on the width
of $R_{n}$ for inputs generated from $\mu_{1}^{\otimes n}$, $\mu_{2}^{\otimes
n}$, and in fact from any $\mu^{\otimes n}$ that shares the same marginals.

We extend and formalize this notion of the minimum width required for any pair
of marginal distributions $\eta_{a}$ over treatment values and $\eta_{b}$ over
control values. For instance, for the distribution $\mu_{1}$, $\eta_{a}$ was
$(1/3,2/3)$ (denoting the probabilities of being $0$ and $1$) and $\eta_{b}$
was $(2/3,1/3)$. For a particular estimate $r\in\\{-1,0,1\\}$ we analyze the
maximum deviation of its rank from $n/2$ across all joint distributions with
the same marginals $\eta_{a},\eta_{b}$. We call this the _width-of- $r$_ for
$\eta_{a},\eta_{b}$, and the minimum width-of-$r$ as the _minimum median
width_ of $\eta_{a},\eta_{b}$. See Definition 2.9 for precise definitions. We
capture the above via our notion of _variability_ (see Definition 2.8) which
is a function of only the marginals and the estimate. Informally, it measures
how much the rank of a particular estimate $r$ can vary across joint
distributions consistent with the given marginals. Moreover, we show that the
minimum median width is always upper bounded by
$\frac{1}{2}\cdot\frac{(2k-3)}{(2k-1)}$ when $\mathbf{a},\mathbf{b}$ are
$k$-ary (see Lemma 3.3). For $k=2$, this value turns out to be $1/6$ and the
distributions described above are thus also the hard instances. Hard instances
that obtain the upper bound on the minimum width for every $k$ are also
detailed in Section 3.

So far, we have shown that variability (or minimum median width) essentially
serves as a limit to MTE estimation. We argue, perhaps surprisingly that,
variability also lends a tight computational handle on the MTE in the
following sense. First, in Section 4, we provide an algorithm that computes
variability exactly in time $O(k)$. We employ a greedy strategy in trying to
identify an appropriate joint distribution that explains the marginals and
also allows us to compute the variability. Second, we show the algorithm that
outputs the minimum median width of the marginals estimated from the observed
responses is instance optimal. The algorithm is called Median-Estimator
(Algorithm 1) and the precise statement of optimality is given in Theorem 3.2.

###  1.2 Our Contributions

In this section, we highlight our key results. For the sake of brevity, we
assume the following notation for all the below-stated results. Let
$\textsc{PO}_{k}$ denote the set of $k$-ary potential outcomes, $\eta$ be any
joint distribution on $\textsc{PO}_{k}\times\textsc{PO}_{k}$ with marginals
$\eta_{a}$ and $\eta_{b}$ for the treatment and control respectively, and
$\epsilon^{*}(\eta_{a},\eta_{b})$ denote the minimum median width (see
Definition 2.9).

In the _instance-optimality_ result described below, we show that any
algorithm, encoded using an RFDT $R_{n}$, has an expected error that is at
least the minimum median width $\epsilon^{*}(\eta_{a},\eta_{b})$. In other
words, $\epsilon^{*}(\eta_{a},\eta_{b})$ essentially captures the minimum
possible approximation error of any algorithm.

######  Theorem 1.1 (Lower Bound on Width of any RFDT (Informal)).

Let $R_{n}$ be an RFDT with width $\epsilon^{R_{n}}(\mathbf{a},\mathbf{b})$
where $(\mathbf{a,b})$ are the treatment and control vectors. If for all
$\mathbf{(a,b)}\in\textsc{PO}_{k}^{n}\times\textsc{PO}_{k}^{n}$, with high
probability $R_{n}$ outputs an median estimate $\widehat{m}$ and width
$\epsilon$ such that
$\widehat{m}\in(\mathbf{a-b})^{0.5+\epsilon}_{0.5-\epsilon}$, then, we have:

$\operatorname*{{\mathbb{E}}}\limits_{(\mathbf{a},\mathbf{b})\leftarrow\eta^{\otimes
n}}[\epsilon^{R_{n}}(\mathbf{a},\mathbf{b})]\gtrsim\epsilon^{*}(\eta_{a},\eta_{b})$

Let $\psi_{k}=\frac{1}{2}\cdot\frac{2k-3}{2k-1}$. The theorem below states
that the minimum median width $\epsilon^{*}(\eta_{a},\eta_{b})$ is always
upper bounded by $\psi_{k}$ and also shows a hardness result by identifying
marginals $\bar{\eta}_{a},\bar{\eta}_{b}$ that satisfy the upper bound.

######  Theorem 1.2 (Tight Bounds for Minimum Median Width (Informal)).

For any $\eta_{a},\eta_{b}$, we have: $\epsilon^{*}(\eta_{a},\eta_{b})$
$\leq\psi_{k}$. And $\exists\bar{\eta}_{a},\bar{\eta}_{b}$, such that
$\epsilon^{*}(\bar{\eta}_{a},\bar{\eta}_{b})\geq\psi_{k}$, which implies
$\epsilon^{*}(\bar{\eta}_{a},\bar{\eta}_{b})=\psi_{k}$. Hence, when $k=2$, by
Theorem 1.1 and the above, $(1/6)^{th}$ is the fundamental limit on
approximating the MTE in the quantile sense.

Finally, given observed potential outcome vectors (obtained using a Bernoulli
design) $\widehat{\mathbf{a}},\widehat{\mathbf{b}}$, we provide a highly
efficient greedy algorithm with linear running time, that returns a median
estimate with a width roughly equal to the minimum median width of the
empirical frequency distribution. In light of Theorem 1.1, we know that
essentially this is optimal.

######  Theorem 1.3 (Median Estimator Algorithm (Informal)).

Consider the _Bernoulli Design_ for assignment, i.e., each individual is
assigned to control/treatment groups independently with probability $1/2$. Fix
arbitrary treatment and control vectors $\mathbf{a,b}$ and let the observed
response vectors be $\widehat{\mathbf{a}},\widehat{\mathbf{b}}$. Then, there
is an algorithm that returns a median estimator $\widehat{m}$, with width
$\epsilon\approx\epsilon^{*}(\eta_{a},\eta_{b})$ and runs in $O(n+k^{2})$
time. Further, with high probability
$\widehat{m}\in(\mathbf{a-b})^{0.5+\epsilon}_{0.5-\epsilon}$.

##  2 Preliminaries and Problem Formulation

We use $a$ (and $b$) to denote the treatment (and control) value for a
particular individual and $\mathbf{a}\ (\text{and }\mathbf{b})$ to denote the
treatment (and control) values of all the $n$ individuals. As we do not
observe all the treatment and control response values for every individual, we
denote the partially observed responses using $\widehat{\mathbf{a}}$ and
$\widehat{\mathbf{b}}$. Let $[n]=\\{1,2,\ldots,n\\}$ and let $\forall
x\in\mathbf{R}$, $(x)^{+}=\max\\{x,0\\}$. For a distribution $\eta$ we use
$\eta^{\otimes n}$ to denote the $n$-fold product distribution.

######  Definition 2.1 ($k$-ary outcomes).

Let
$\textsc{PO}_{k}\stackrel{{\scriptstyle{def}}}{{=}}\text{\\{}0,1,\cdots,k-1\\}$
be the set of all $k$-ary outcomes for both treatment and control groups,
i.e., $\mathbf{a}$ and $\mathbf{b}$. Each potential outcome value for an
individual $i\in[n]$ lies in the set $\textsc{PO}_{k}$, i.e.,
$a_{i},b_{i}\in\textsc{PO}_{k}\quad\forall i\in[n]$. Therefore,
$\textsc{ITE}_{i}\in\\{-(k-1),\cdots,k-1\\}\ \forall i\in[n]$.

We make the following assumptions that are common in causality literature.

######  Assumption 2.2 (Independence).

The experimental design is independent of the potential outcomes
${\mathbf{a}},{\mathbf{b}}$.

######  Assumption 2.3 (SUTVA).

The treatment assignment of a single individual does not affect the outcomes
of other units.

Now, we begin with the definition of a quantile of any vector.

######  Definition 2.4 (Quantiles of a vector).

Let $\mathbf{v}\in\mathbb{R}^{n}$ be a vector, and let $\ell,u\in[0,1]$ be
numbers such that $\ell\leq u$. For any $r\in\mathbb{R}$ define
$q_{\ell}(r,\mathbf{v})\coloneqq\frac{\lvert\\{i\in[n]\mid\mathbf{v}_{i}<r\\}\rvert}{n}$
and $q_{u}(r,\mathbf{v})\coloneqq\frac{\lvert\\{i\in[n]\mid\mathbf{v}_{i}\leq
r\\}\rvert}{n}$. Finally, define the set of _$\ell$ -to-$u$ quantiles of
$\mathbf{v}$_ as:

$\big{(}\mathbf{v}\big{)}_{\ell}^{u}\coloneqq\\{r\in\mathbb{R}\mid[\ell,u]\cap[q_{\ell}(r,\mathbf{v}),q_{u}(r,\mathbf{v})]\neq\emptyset\\}.$

Next, we describe the concept of variability. Before that, we introduce some
helpful related notions. Let $\eta$ be a joint distribution over
$\textsc{PO}_{k}\times\textsc{PO}_{k}$. We interpret it as a joint
distribution on $a,b$, i.e.,

$\eta_{xy}\stackrel{{\scriptstyle{def}}}{{=}}\text{P}r[a=x,b=y].$

The marginal distribution denoted by $\eta_{a}$ is defined as: $\forall
x\in\textsc{PO}_{k}:\eta_{a}(x)\stackrel{{\scriptstyle{def}}}{{=}}\text{}\sum\limits_{y\in\textsc{PO}_{k}}\eta_{xy}$.
$\eta_{b}$ is defined similarly. Let $\Delta_{k}$ be the simplex of all
probability distributions over $\textsc{PO}_{k}$:

$\Delta_{k}\stackrel{{\scriptstyle{def}}}{{=}}\text{}\left\\{(g_{0},g_{1},\ldots,g_{k-1})\mid
g_{i}\geq 0\text{ for all }i,\text{ and }\sum_{i=0}^{k-1}g_{i}=1\right\\}.$

Given marginal distributions $\eta_{a}$ and $\eta_{b}$, we define
$\mathcal{J}(\eta_{a},\eta_{b})$ as the set of joint distributions $\eta$ with
marginals $\eta_{a}$ and $\eta_{b}$ for the potential outcomes.

######  Definition 2.5 (Lower Quantile).

Given a joint distribution $\eta$, the lower quantile of a particular value
$r\in\\{-(k-1),\cdots,k-1\\}$ in the range of ITE, is defined as the total
probability mass of all possible
$(x,y)\in\textsc{PO}_{k}\times\textsc{PO}_{k}$ that have an ITE value strictly
less than $r$. Formally, we have:

$q_{\ell}(r,\eta)\stackrel{{\scriptstyle{def}}}{{=}}\text{}\sum_{x\in\textsc{PO}_{k},y\in\textsc{PO}_{k}}\mathbf{1}\\{x-y<r\\}\eta_{xy}.$

######  Definition 2.6 (Upper Quantile).

Given a joint distribution $\eta$, the upper quantile of a particular value
$r\in\\{-(k-1),\cdots,k-1\\}$ in the range of ITE, is defined as the total
probability mass of all possible
$(x,y)\in\textsc{PO}_{k}\times\textsc{PO}_{k}$ that have an ITE value of at
most $r$. Formally, we have:

$q_{u}(r,\eta)\stackrel{{\scriptstyle{def}}}{{=}}\text{}\sum_{x\in\textsc{PO}_{k},y\in\textsc{PO}_{k}}\mathbf{1}\\{x-y\leq
r\\}\eta_{xy}.$

Next, we prove a proposition that is useful when we want to translate between
the lower quantile of a vector and a lower quantile of a distribution. It
roughly states that for a typical sample $(\mathbf{a,b})$ from $\eta^{\otimes
n}$ (where $\eta$ is a joint distribution over
$\textsc{PO}_{k}\times\textsc{PO}_{k}$) we have
$q_{\ell}(r,\eta)=q_{\ell}(r,\mathbf{a}-\mathbf{b})$ and
$q_{u}(r,\eta)=q_{u}(r,\mathbf{a}-\mathbf{b})$. Further,
$(\mathbf{a,b})\in\textsc{PO}_{k}^{n}\times\textsc{PO}_{k}^{n}$ there is a
joint distribution $\eta$ such that $(\mathbf{a,b})$ is typical for
$\eta^{\otimes n}$.

######  Proposition 2.7.

Let $\eta$ be a joint distribution over
$\textsc{PO}_{k}\times\textsc{PO}_{k}$. A sample $(\mathbf{a,b})$ from
$\eta^{\otimes n}$ is said to be typical if for all
$(x,y)\in\textsc{PO}_{k}\times\textsc{PO}_{k}$ we have that the empirical
frequency of $\rho(x,y)=\frac{|i\in[n]:(a_{i},b_{i})=(x,y)|}{n}$ is
$\eta_{xy}$. Suppose $\mathbf{a,b}$ is a typical sample from $\eta^{\otimes
n}$, then, for all $r\in\mathbb{R}$ we have
$q_{\ell}(r,\eta)=q_{\ell}(r,\mathbf{a}-\mathbf{b})$ and
$q_{u}(r,\eta)=q_{u}(r,\mathbf{a}-\mathbf{b})$. Further, for any
$(\mathbf{a,b})\in\textsc{PO}_{k}^{n}\times\textsc{PO}_{k}^{n}$ there is a
joint distribution $\eta$ such that $(\mathbf{a,b})$ is typical for
$\eta^{\otimes n}$.

###### Proof.

Suppose $\mathbf{a,b}$ is a typical sample from $\eta^{\otimes n}$. Let
$\mathbf{v=a-b}$. Then,

$\displaystyle
q_{\ell}(r,\mathbf{v})=\frac{\lvert\\{i\in[n]\mid\mathbf{v}_{i}<r\\}\rvert}{n}=\sum_{x\in\textsc{PO}_{k},y\in\textsc{PO}_{k}}\mathbf{1}\\{x-y<r\\}\rho_{xy}$
$\displaystyle=\sum_{x\in\textsc{PO}_{k},y\in\textsc{PO}_{k}}\mathbf{1}\\{x-y<r\\}\eta_{xy}$
$\displaystyle=q_{\ell}(r,\eta).$

A similar argument holds for $q_{u}(r,\eta)=q_{u}(r,\mathbf{a}-\mathbf{b})$.
Next, given any
$(\mathbf{a,b})\in\textsc{PO}_{k}^{n}\times\textsc{PO}_{k}^{n}$ we associate a
joint distribution $\eta$ over $\textsc{PO}_{k}\times\textsc{PO}_{k}$ as
follows: sample $i\in[n]$ uniformly at random and output $(a_{i},b_{i})$. It
is easy to see that $(\mathbf{a,b})$ is typical for $\eta^{\otimes n}$. ∎

Next, we define the concept of variability which captures how much the
quantiles of a given estimate $r$ can vary across all joint distributions that
explain the given marginals.

######  Definition 2.8 (Variability).

Variability is a function
$\nu:\mathbb{R}\times\Delta_{k}\times\Delta_{k}\rightarrow[0,1]\times[0,1]$,
that takes as input an estimate $r$ and marginal distributions for the
potential outcomes $\eta_{a},\eta_{b}$, and returns a tuple as follows:

$\nu(r,\eta_{a},\eta_{b})\stackrel{{\scriptstyle{def}}}{{=}}\text{}(\nu_{l}(r,\eta_{a},\eta_{b}),\nu_{u}(r,\eta_{a},\eta_{b})),\text{
where:}$ $\displaystyle\nu_{l}(r,\eta_{a},\eta_{b})$
$\displaystyle=\max_{\eta\in\mathcal{J}(\eta_{a},\eta_{b})}q_{l}(r,\eta)$
$\displaystyle\nu_{u}(r,\eta_{a},\eta_{b})$
$\displaystyle=\min_{\eta\in\mathcal{J}(\eta_{a},\eta_{b})}q_{u}(r,\eta).$

Next, we understand that given the variability of an estimate $r$, how far can
it be from the median in a quantile sense.

######  Definition 2.9 (Minimum Median Width of Marginals).

Let $\eta_{a}$ and $\eta_{b}$ be distributions over $\textsc{PO}_{k}$.
Further, for any $r\in\mathbb{R}$ the _width-of- $r$_ for the pair
$(\eta_{a},\eta_{b})$ is defined as

$\varepsilon(r,\eta_{a},\eta_{b}):=\max\\{\left(\nu_{l}(r,\eta_{a},\eta_{b})-1/2\right)^{+},(1/2-\nu_{u}(r,\eta_{a},\eta_{b}))^{+}\\}.$

Then, the _minimum median width_ of the pair $(\eta_{a},\eta_{b})$, denoted by
$\epsilon^{*}(\eta_{a},\eta_{b})$, is defined as

$\min_{r\in\\{-(k-1),\cdots,(k-1)\\}}\varepsilon(r,\eta_{a},\eta_{b}).$

######  Remark 2.10.

Notice that for any sample $(\mathbf{a},\mathbf{b})$ from $\eta^{\otimes n}$
such that the marginals of $\eta$ are $\eta_{a}$ and $\eta_{b}$ we have that
$r\in(\mathbf{a-b})^{0.5+\epsilon}_{0.5-\epsilon}$ where
$\epsilon=\varepsilon(r,\eta_{a},\eta_{b})$.

###  2.1 Formalizing the Algorithm

Next, we formalize the meaning of an algorithm that respects the inference
constraints of observing either $a_{i}$ or $b_{i}$ but not both. At the end of
observing the vectors $a$ and $b$ in a random manner (even in an adaptive
manner), it outputs a member $m$ and a parameter $\epsilon\in[0,0.5]$ and the
output is correct if
$m\in\big{(}\mathbf{a-b}\big{)}_{0.5-\epsilon}^{0.5+\epsilon}$.

We first formalize the case of a deterministic algorithm, viewed as a decision
tree.

######  Definition 2.11 (Feasible Decision Tree (FDT)).

A _feasible decision tree_ (FDT) for block-length $n$ is a rooted and labeled
$k$-ary tree of depth $n$. Each non-leaf node is labeled with some element
from the set $\\{a,b\\}\times[n]$, and each leaf node is labeled with an
element from $\mathbb{R}\times[0,0.5]$. The edges of the tree are labeled with
elements from $\textsc{PO}_{k}$ such that from any non-leaf node the set of
labels on the edges to its children is $\textsc{PO}_{k}$. Finally, the tree
should also satisfy the following _feasibility constraint_. For any path $P$
from the root to some leaf and for all $i\in[n]$:

$\lvert\\{\operatorname{lab}(u)\mid u\in P\text{ and is a non-leaf
node}\\}\cap(\\{a,b\\}\times\\{i\\})\rvert=1.$

Here, $\operatorname{lab}(u)$ refers to the label of a vertex $u$.

######  Definition 2.12 (Output of an FDT).

Let $A_{n}$ be an FDT for block-length $n$. On input
$(\mathbf{a,b})\in\textsc{PO}_{k}^{n}\times\textsc{PO}_{k}^{n}$ the output of
$A_{n}$, denoted as $A_{n}(\mathbf{a,b})$, is the label of the unique leaf
obtained by following the path from the root to the leaf, where at each non-
leaf node $u$, labeled $(c,i)$, we follow the edge from $u$ labeled with the
same value as $c_{i}$: for instance, if $c=a$ and $a_{i}=0$, then we follow
the edge from $u$ which is labeled $0$. So, the output is a tuple
$(m,\epsilon)\in\mathbb{R}\times[0,0.5]$.

Next, we define the case of a randomized algorithm.

######  Definition 2.13 (Randomized FDT (RFDT)).

A randomized feasible decision tree (RFDT) for block-length $n$ is a
distribution over the set of all FDTs for block-length $n$.

######  Remark 2.14.

Definition 2.13 is equivalent to the description of RFDTs given in section
1.1.1. To see this, we can imagine pre-sampling all the internal randomness as
required in the description of an RFDT in Section 1.1.1. In particular, for
all non-leaf nodes $u$ of the RFDT, we fix a sample according to $\gamma_{u}$
before proceeding along the RFDT. Given the samples from $\gamma_{u}$ for each
$u$, the RFDT is now actually an FDT. Hence, overall the RFDT is a
distribution over FDTs.

######  Definition 2.15 (Output of an RFDT).

Let $R_{n}$ be an RFDT for block-length $n$ and let $\mu_{n}$ be its
distribution over the set of FDTs. On input
$(\mathbf{a,b})\in\textsc{PO}_{k}^{n}\times\textsc{PO}_{k}^{n}$ the output of
$R_{n}$, denoted as $R_{n}(\mathbf{a,b})$, is the random variable obtained in
the following manner. Sample an FDT $A_{n}$ according to $\mu_{n}$ and output
$A_{n}(\mathbf{a,b})$.

Next, we define the notions of width and error of an RFDT.

######  Definition 2.16 (Error and Width of an RFDT).

Fix $\delta\in[0,1]$. An RFDT,$R_{n}$, for block-length $n$, is said to be a
median approximator with error $\delta$ or a $\delta$-error median
approximator, if for all
$(\mathbf{a},\mathbf{b})\in\textsc{PO}_{k}^{n}\times\textsc{PO}_{k}^{n}$:

$\Pr\limits_{(m,\epsilon)\leftarrow
R_{n}(\mathbf{a,b})}\left[m\in\big{(}a-b\big{)}_{0.5-\epsilon}^{0.5+\epsilon}\right]\geq
1-\delta.$

Further, the width of $R_{n}$ on input $(\mathbf{a},\mathbf{b})$, denoted as
$\epsilon^{R_{n}}(\mathbf{a},\mathbf{b})$, is defined as
$\operatorname*{{\mathbb{E}}}\limits_{(m,\epsilon)\leftarrow
R_{n}(\mathbf{a,b})}[\epsilon]$.

##  3 Variability & Minimum Median Width

We first prove Lemma 3.1 which roughly states that any RFDT is incapable of
distinguishing between two joint distributions over
$\textsc{PO}_{k}\times\textsc{PO}_{k}$, say $\eta$ and $\tau$, which have the
same marginals, i.e., $\eta_{a}=\tau_{a}$ and $\eta_{b}=\tau_{b}$.
Specifically, any RFDT $R_{n}$ is unable to distinguish whether the input is
drawn from the distribution $\eta^{\otimes n}$ or $\tau^{\otimes n}$.

The proof proceeds by noting that any RFDT $R_{n}$ is a distribution over
FDTs. Further, for any FDT $A_{n}$ and for any path $P$ in $A_{n}$ the set of
inputs $(\mathbf{a},\mathbf{b})$ where $A_{n}$ follows the path $P$, denoted
by $F(A_{n},P)$, is described by constraints on exactly one of $a_{i}$ or
$b_{i}$ for $i\in[n]$. Hence,
$\Pr_{(\mathbf{a},\mathbf{b})\leftarrow\eta^{\otimes
n}}[(\mathbf{a},\mathbf{b})\in F(A_{n},P)]$ is only a function of the
marginals $\eta_{a},\eta_{b}$. Thus, the probability of $F(A_{n},P)$ remains
the same under $\tau$.

######  Lemma 3.1.

Let $\eta,\tau$ be joint distributions over
$\textsc{PO}_{k}\times\textsc{PO}_{k}$ with same marginals for $a,b$, i.e.,
$\eta_{a}=\tau_{a}$ and $\eta_{b}=\tau_{b}$, and let $R_{n}$ be an RFDT for
block-length $n$. Further, let $R_{n}(\eta),R_{n}(\tau)$ be the distributions
of the output of $R_{n}$ when the input $(\mathbf{a,b})$ is sampled according
to $\eta^{\otimes n}$ and $\tau^{\otimes n}$ respectively. Then,
$R_{n}(\eta)=R_{n}(\tau)$. More precisely, for every
$(m^{\prime},\epsilon^{\prime})\in\mathbb{R}\times[0,0.5]$ we have

$\Pr\limits_{\begin{subarray}{c}(\mathbf{a},\mathbf{b})\leftarrow\eta^{\otimes
n}\\\ (m,\epsilon)\leftarrow
R_{n}(a,b)\end{subarray}}[(m,\epsilon)=(m^{\prime},\epsilon^{\prime})]=\Pr\limits_{\begin{subarray}{c}(\mathbf{a},\mathbf{b})\leftarrow\tau^{\otimes
n}\\\ (m,\epsilon)\leftarrow
R_{n}(a,b)\end{subarray}}[(m,\epsilon)=(m^{\prime},\epsilon^{\prime})].$

###### Proof.

Fix some $(m^{\prime},\epsilon^{\prime})$. Let $\operatorname{supp}(R_{n})$
denote the support of the distribution of $R_{n}$, i.e., the support of the
distribution over FDTs of block-length $n$ which represents $R_{n}$. Fix any
FDT $A_{n}\in\operatorname{supp}(R_{n})$ and any path from the root to leaf in
$A_{n}$, say $P$, such that the label of the leaf node is
$(m^{\prime},\epsilon^{\prime})$. Further, let $F(A_{n},P)$ denote those
inputs,
$(\mathbf{a},\mathbf{b})\in\textsc{PO}_{k}^{n}\times\textsc{PO}_{k}^{n}$ such
that $A_{n}$ follows the path $P$ on the input $(\mathbf{a},\mathbf{b})$. We
will show that

$\Pr\limits_{\mathbf{a},\mathbf{b}\leftarrow\eta^{\otimes
n}}[F(A_{n},P)]=\Pr\limits_{\mathbf{a},\mathbf{b}\leftarrow\tau^{\otimes
n}}[F(A_{n},P)].$

To see this, note that there is a function $f:[n]\to\textsc{PO}_{k}$ and a
subset $S\subseteq[n]$ such that $F(A_{n},P)=\\{(\mathbf{a},\mathbf{b})\mid
a_{i}=f(i)\text{ for }i\in S,b_{i}=f(i)\text{ for }i\notin S\\}$. This follows
from the feasibility constraint on $A_{n}$ as mentioned in Definition 2.11.
Hence, we have

$\displaystyle\Pr\limits_{\mathbf{a},\mathbf{b}\leftarrow\eta^{\otimes
n}}[F(A_{n},P)]$ $\displaystyle=\prod_{i\in
S}\eta_{a}(f(i))\times\prod_{i\notin S}\eta_{b}(f(i))$
$\displaystyle=\prod_{i\in S}\tau_{a}(f(i))\times\prod_{i\notin
S}\tau_{b}(f(i))=\Pr\limits_{\mathbf{a},\mathbf{b}\leftarrow\tau^{\otimes
n}}[F(A_{n},P)].$

Taking a union over all paths $P$ in $A_{n}$ whose leaves are labelled with
$(m^{\prime},\epsilon^{\prime})$ it follows that

$\Pr\limits_{\mathbf{a},\mathbf{b}\leftarrow\eta^{\otimes
n}}[A_{n}((\mathbf{a},\mathbf{b}))=(m^{\prime},\epsilon^{\prime})]=\Pr\limits_{\mathbf{a},\mathbf{b}\leftarrow\tau^{\otimes
n}}[A_{n}((\mathbf{a},\mathbf{b}))=(m^{\prime},\epsilon^{\prime})].$

Since the above statement is true for every fixed $A_{n}$ in the support of
$R_{n}$ the statement of the lemma follows. ∎

Using Lemma 3.1 we prove the theorem below which serves to capture the
hardness of approximation the MTE in terms of the variability of the marginal
distributions of the treatment and control vector. As a consequence, if we can
devise an algorithm for computing $\epsilon^{*}(\eta_{a},\eta_{b})$ exactly
(section Section 4), then, then our algorithm is essentially optimal.

The proof idea is as follows. Fix a $\delta$-error median approximator
$R_{n}$. Given a joint distribution $\eta$ which has marginals $\eta_{a}$ and
$\eta_{b}$, for an estimate $r\in\\{-(k-1),\ldots,k-1\\}$, let $\gamma_{r}$
denote the probability that $R_{n}$ outputs $m=r$ and
$\epsilon<\epsilon(r,\eta_{a},\eta_{b})$. Now, consider the joint distribution
$\eta^{(r)}$ which witnesses the width-of-$r$ value, i.e.,such that at least
one of $\epsilon(r)=\left(q_{\ell}(r,\eta^{(r)})-1/2\right)^{+}$ or
$\epsilon(r)=(1/2-q_{u}(r,\eta^{(r)}))^{+}$ is true. For $\eta^{(r)}$ the
probability that $R_{n}$ outputs $m=r$ and
$\epsilon<\epsilon(r,\eta_{a},\eta_{b})$ is still $\gamma_{r}$. However if
$\epsilon<\epsilon(r,\eta_{a},\eta_{b})$, then, for a typical sample from
$\left(\eta^{(r)}\right)^{\otimes n}$ it will be the case that
$r\notin(\mathbf{a-b})^{0.5+\epsilon}_{0.5-\epsilon}$, and this counts as an
error. Hence, $\gamma_{r}$ is not too much larger than $\delta$. Overall,
taking a union over all values for $r$, we get that except for around a
$2k\delta$ probability the output of $R_{n}$ has to be at least
$\epsilon^{*}(\eta_{a},\eta_{b})$.

######  Theorem 3.2 (Lower Bound on Width of any RFDT).

Let $\eta$ be any joint distribution on $\textsc{PO}_{k}\times\textsc{PO}_{k}$
with marginals $\eta_{a}$ and $\eta_{b}$ for the treatment and control
respectively. Further, let $R_{n}$ be a $\delta$-error median approximator for
block-length $n$. Then, for any $\beta>0$ we have

$\operatorname*{{\mathbb{E}}}\limits_{(\mathbf{a},\mathbf{b})\leftarrow\eta^{\otimes
n}}[\epsilon^{R_{n}}(\mathbf{a},\mathbf{b})]\geq(\epsilon^{*}(\eta_{a},\eta_{b})-\beta)\cdot(1-2k\delta)+2k\exp(-2\beta^{2}\cdot
n).$

###### Proof.

For $r\in\\{-(k-1),\ldots,k-1\\}$ recall that (see width-of-$r$ in Definition
2.9)

$\varepsilon(r,\eta_{a},\eta_{b})=\max\\{\left(\nu_{\ell}(r,\eta_{a},\eta_{b})-1/2\right)^{+},(1/2-\nu_{u}(r,\eta_{a},\eta_{b}))^{+}\\},$

where $\nu_{\ell}$ and $\nu_{u}$ are as defined in Definition 2.8. For the
remainder of the proof we will use $\epsilon(r)$ and $\epsilon^{*}$ to denote
$\epsilon(r,\eta_{a},\eta_{b})$ and $\epsilon^{*}(\eta_{a},\eta_{b})$
respectively. Note that $\epsilon^{*}=\min_{r}\epsilon(r)$.

Further, let $\eta^{(r)}$ denote a joint distribution over
$\textsc{PO}_{k}\times\textsc{PO}_{k}$ with marginals $\eta_{a},\eta_{b}$ such
that at least one of $\epsilon(r)=\left(q_{\ell}(r,\eta^{(r)})-1/2\right)^{+}$
or $\epsilon(r)=(1/2-q_{u}(r,\eta^{(r)}))^{+}$ is true. Now, let

$\gamma_{r}=\Pr\limits_{\begin{subarray}{c}(\mathbf{a},\mathbf{b})\leftarrow\eta^{\otimes
n}\\\ (m,\epsilon)\leftarrow
R_{n}(\mathbf{a},\mathbf{b})\end{subarray}}[m=r,\epsilon<\epsilon(r)-\beta].$

Notice that by Lemma 3.1 we have

$\gamma_{r}=\Pr\limits_{\begin{subarray}{c}(\mathbf{a},\mathbf{b})\leftarrow\left(\eta^{(r)}\right)^{\otimes
n}\\\ (m,\epsilon)\leftarrow
R_{n}(\mathbf{a},\mathbf{b})\end{subarray}}[m=r,\epsilon<\epsilon(r)-\beta],$

i.e., the value of $\gamma_{r}$ remains unchanged when inputs are sampled from
$\left(\eta^{(r)}\right)^{\otimes n}$ instead of $\eta^{\otimes n}$.

Now, for an input $(\mathbf{a},\mathbf{b})$ sampled according to
$\left(\eta^{(r)}\right)^{\otimes n}$ we have with at least
$1-\exp(-2\beta^{2}\cdot n)$ probability that
$r\notin(\mathbf{a}-\mathbf{b})^{0.5+\epsilon(r)-\beta}_{0.5-\epsilon(r)+\beta}$.
This follows from the Chernoff bound. Suppose
$\epsilon(r)=\left(q_{\ell}(r,\eta^{(r)})-1/2\right)$, and note that with
probability at least $1-\exp(-2\beta^{2}\cdot n)$ the empirical frequency of
ITE’s less than $r$ will be at least $q_{\ell}(r,\eta^{(r)})-\beta$, and hence
$q_{\ell}(r,\mathbf{a-b})\geq q_{\ell}(r,\eta^{(r)})-\beta$ with probability
at least $1-\exp(-2\beta^{2}\cdot n)$.

Notice that this yields $\gamma_{r}\leq\delta+\exp(-2\beta^{2}\cdot n)$ as

$\displaystyle\gamma_{r}$
$\displaystyle=\Pr\limits_{\begin{subarray}{c}(\mathbf{a},\mathbf{b})\leftarrow\left(\eta^{(r)}\right)^{\otimes
n}\\\ (m,\epsilon)\leftarrow
R_{n}(\mathbf{a},\mathbf{b})\end{subarray}}[m=r,\epsilon<\epsilon(r)-\beta]$
$\displaystyle\leq\Pr\limits_{\begin{subarray}{c}(\mathbf{a},\mathbf{b})\leftarrow\left(\eta^{(r)}\right)^{\otimes
n}\\\ (m,\epsilon)\leftarrow
R_{n}(\mathbf{a},\mathbf{b})\end{subarray}}[m=r,\epsilon<\epsilon(r)-\beta,r\notin(\mathbf{a}-\mathbf{b})^{0.5+\epsilon(r)-\beta}_{0.5-\epsilon(r)+\beta}]+\exp(-2\beta^{2}\cdot
n)$
$\displaystyle\leq\Pr\limits_{\begin{subarray}{c}(\mathbf{a},\mathbf{b})\leftarrow\left(\eta^{(r)}\right)^{\otimes
n}\\\ (m,\epsilon)\leftarrow
R_{n}(\mathbf{a},\mathbf{b})\end{subarray}}[m\notin(\mathbf{a}-\mathbf{b})^{0.5+\epsilon}_{0.5-\epsilon}]+\exp(-2\beta^{2}\cdot
n)\leq\delta+\exp(-2\beta^{2}\cdot n).$

Thus,

$\displaystyle\operatorname*{{\mathbb{E}}}\limits_{(\mathbf{a},\mathbf{b})\leftarrow\eta^{\otimes
n}}[\epsilon^{R_{n}}(\mathbf{a},\mathbf{b})]$
$\displaystyle\geq(\epsilon^{*}-\beta)\times\Pr\limits_{\begin{subarray}{c}(\mathbf{a},\mathbf{b})\leftarrow\eta^{\otimes
n}\\\ (m,\epsilon)\leftarrow
R_{n}(\mathbf{a},\mathbf{b})\end{subarray}}[\epsilon\geq\epsilon^{*}-\beta],$
which by a union bound over all $r\in\textsc{ITE}$ gives
$\displaystyle\geq(\epsilon^{*}-\beta)\times\left(1-\sum_{r\in\\{-(k-1),\ldots,k-1\\}}\gamma_{r}\right)$
$\displaystyle\geq(\epsilon^{*}(\eta_{a},\eta_{b})-\beta)\cdot(1-2k\delta)+2k\exp(-2\beta^{2}\cdot
n).$

∎

###  3.1 Bounds on Variability

In this section, we show that the value of $\epsilon^{*}(\eta_{a},\eta_{b})$
is upper bounded by $\frac{1}{2}\cdot\frac{(2k-3)}{(2k-1)}$. In doing so, we
have essentially shown that any algorithm that computes the variability
exactly obtains an error that is upper bounded by
$\frac{1}{2}\cdot\frac{(2k-3)}{(2k-1)}$.

######  Lemma 3.3 (Upper Bound on Minimum Median Width).

Given two distributions $\eta_{a}$ and $\eta_{b}$ over $\textsc{PO}_{k}$, we
have
$\epsilon^{*}(\eta_{a},\eta_{b})\leq\frac{1}{2}\cdot\frac{(2k-3)}{(2k-1)}$.

###### Proof.

For contradiction, let’s assume that $\exists\eta_{a},\eta_{b}\text{ such that
}\epsilon^{*}(\eta_{a},\eta_{b})>\frac{1}{2}\cdot\frac{(2k-3)}{(2k-1)}$. For
brevity let $\psi_{k}=\frac{1}{(2k-1)}$ and let
$\epsilon(r,\eta_{a},\eta_{b})$, $\nu_{\ell}(r,\eta_{a},\eta_{b})$ and
$\nu_{u}(r,\eta_{a},\eta_{b})$ be denoted as $\epsilon(r),\nu_{\ell}(r)$ and
$\nu_{u}(r)$ respectively. Therefore, we have:

$\varepsilon(r)>1/2-\psi_{k}\quad\forall r\in\\{-(k-1),\cdots,(k-1)\\},$

where from Definition 2.9

$\varepsilon(r)=\max\\{\left(\nu_{l}(r)-1/2\right)^{+},(1/2-\nu_{u}(r))^{+}\\}.$

Notice that $\nu_{\ell}(-(k-1))=0$, thus, $\nu_{u}(-(k-1))<\psi_{k}$. However,
$\nu_{u}(k-1)=1$, thus $\nu_{\ell}((k-1))>1-\psi_{k}$. Now, let
$r\in\textsc{ITE}$ be the smallest value such that $\nu_{\ell}(r)>1-\psi_{k}$,
therefore, $\nu_{\ell}(r-1)\leq 1-\psi_{k}$ and hence $\nu_{u}(r-1)<\psi_{k}$.
We show that the existence of such an $r$ is impossible by analyzing
$\operatorname*{{\mathbb{E}}}_{X\leftarrow\eta_{a}}[X]-\operatorname*{{\mathbb{E}}}_{Y\leftarrow\eta_{b}}[Y]$.

Notice that for any joint distribution $\eta$ such that its marginals are
$\eta_{a}$ and $\eta_{b}$ and $q_{u}(r-1,\eta)<\psi_{k}$ we have

$\displaystyle\operatorname*{{\mathbb{E}}}_{X\leftarrow\eta_{a}}[X]-\operatorname*{{\mathbb{E}}}_{Y\leftarrow\eta_{b}}[Y]$
$\displaystyle=\operatorname*{{\mathbb{E}}}_{X,Y\leftarrow\eta}[X-Y]>-(k-1)\cdot\psi_{k}+(r)\cdot(1-\psi_{k}).$

Similarly for any joint distribution $\eta^{\prime}$ such that its marginals
are $\eta_{a}$ and $\eta_{b}$ and $q_{\ell}(r,\eta)>1-\psi_{k}$ we have

$\displaystyle\operatorname*{{\mathbb{E}}}_{X\leftarrow\eta_{a}}[X]-\operatorname*{{\mathbb{E}}}_{Y\leftarrow\eta_{b}}[Y]$
$\displaystyle=\operatorname*{{\mathbb{E}}}_{X,Y\leftarrow\eta^{\prime}}[X-Y]<(r-1)\cdot(1-\psi_{k})+(k-1)\cdot\psi_{k}.$

Hence, we obtain

$\displaystyle-(k-1)\cdot\psi_{k}+(r)\cdot(1-\psi_{k})$
$\displaystyle<(r-1)\cdot(1-\psi_{k})+(k-1)\cdot\psi_{k}$
$\displaystyle\implies 1-\psi_{k}<2(k-1)\psi_{k}$
$\displaystyle\implies\frac{2k-2}{2k-1}<\frac{2k-2}{2k-1},$

which is a contradiction. Therefore, the posited $r$ can’t exist and this
completes the proof. ∎

We show a tightness result by arguing that there exists $\eta_{a},\eta_{b}$
that satisfies the stated upper bound.

######  Theorem 3.4 (Extremal distributions for Minimum Median Width).

For any $k\in\mathcal{N}$, $r\in\textsc{PO}_{k}$, let $\eta_{a},\eta_{b}$ be
distributions over $\textsc{PO}_{k}$ defined as:

$\eta_{a}(r)=\frac{1}{2k-1}\cdot(1+\mathbf{1}\\{r>0\\})\mbox{ and
}\eta_{b}(r)=\frac{1}{2k-1}\cdot(1+\mathbf{1}\\{r<k-1\\}).$

Then, $\epsilon^{*}(\eta_{a},\eta_{b})\geq\frac{1}{2}\cdot\frac{2k-3}{2k-1}$.
Combined with Lemma 3.3 we have
$\epsilon^{*}(\eta_{a},\eta_{b})=\frac{1}{2}\cdot\frac{2k-3}{2k-1}$.

###### Proof.

We first show that $\nu_{u}(0,\eta_{a},\eta_{b})\leq\frac{1}{2k-1}$. For this
consider the joint distribution $\eta$ described as follows: sample
$a\sim\eta_{a}$ and let $b=(a-1)\mod k$. It is easy to see that $b$ is
distributed as $\eta_{b}$. Further, $a-b$ is $-(k-1)$ with probability
$\frac{1}{2k-1}$ and is $1$ otherwise. Thus, $q_{u}(0,\eta)=\frac{1}{2k-1}$,
and hence, $\nu_{u}(0,\eta_{a},\eta_{b})\leq\frac{1}{2k-1}$.

Next, we show that $\nu_{\ell}(1,\eta_{a},\eta_{b})\geq\frac{2k-2}{2k-1}$. For
this consider the joint distribution $\eta$ described as follows: sample
$a\sim\eta_{a}$ and if $a\neq k-1$ let $b=a$; if $a=k-1$ then let $b=0\text{
or }k-1$ with probability $0.5$. It is easy to see that $b$ is distributed as
$\eta_{b}$. Further, $a-b$ is $(k-1)$ with probability $\frac{1}{2k-1}$ and is
$0$ otherwise. Thus, $q_{\ell}(1,\eta)=\frac{2k-2}{2k-1}$, and hence,
$\nu_{\ell}(1,\eta_{a},\eta_{b})\geq\frac{2k-2}{2k-1}$.

Notice that both $\nu_{u}(r,\eta_{a},\eta_{b})$ and
$\nu_{\ell}(r,\eta_{a},\eta_{b})$ are increasing functions of $r$. Hence, for
all $r\leq 0$ we have $\nu_{u}(r,\eta_{a},\eta_{b})\leq\frac{1}{2k-1}$ and
hence for all $r\leq 0$ we have $\epsilon(r,\eta_{a},\eta_{b})\geq
0.5\cdot\frac{2k-3}{2k-1}$. Similarly, for all $r\geq 1$ we have
$\nu_{\ell}(r,\eta_{a},\eta_{b})\geq\frac{2k-2}{2k-1}$ and hence for all
$r\geq 1$ we have $\epsilon(r,\eta_{a},\eta_{b})\geq
0.5\cdot\frac{2k-3}{2k-1}$. This proves the theorem. ∎

##  4 Algorithm

In this section, we describe the Algorithm Median-Estimator that takes as
input the observed response vectors $\widehat{\mathbf{a}}$,
$\widehat{\mathbf{b}}$, and returns an estimate $\widehat{m}$ for the median
of the vector $\mathbf{a-b}$, along with an error estimate $\epsilon$; with a
guarantee that $\widehat{m}$ is within an $\epsilon$ band around the median,
i.e.,
$\widehat{m}\in\big{(}\mathbf{a-b}\big{)}_{0.5-\epsilon}^{0.5+\epsilon}.$

First, in section 4.1, we describe an algorithm Median-Estimator that uses the
algorithm Compute-Variability and returns a median estimate $\widehat{m}$ and
its $\epsilon$ band. Next, in section 4.2 we describe an algorithm Compute-
Variability that computes variability $\nu(r,\eta_{a},\eta_{b})$, given a
feasible median value $r\in\\{-(k-1),\cdots,(k-1)\\}$, and marginals
$\eta_{a}$, $\eta_{b}$.

###  4.1 Approximating the Median

In Algorithm 1, we first compute the empirical distributions
$\rho_{a},\rho_{b}$ that are close approximations to $\eta_{a},\eta_{b}$ (up
to Chernoff errors). Then, using Algorithms Compute-Variability-Lower-Quantile
and Compute-Variability-Upper-Quantile, we compute the variability for every
feasible median value in the set $\\{-(k-1),\cdots,(k-1)\\}$. The minimizer in
the set outputs the median estimate $\widehat{m}$ and the corresponding
minimum median width (see defn:2.9) $\epsilon(\widehat{m},\rho_{a},\rho_{b})$
as the approximation error.

Algorithm 1 Median-Estimator

1: Input: Treatment and control response vectors
$\widehat{\mathbf{a}},\widehat{\mathbf{b}}$ and slack parameter $\beta>0$

2: Output: Median estimate $\widehat{m}$ and the width $\epsilon$

3: Compute the empirical distributions $\rho_{a}$, $\rho_{b}$ as follows:
$\displaystyle\forall j\in\textsc{PO}_{k}\quad\rho_{a}(j)$
$\displaystyle=\frac{|\\{i|\widehat{\mathbf{a}}_{i}=j\text{ for
}i\in\\{1,2,\cdots,n\\}\\}|}{n/2}$ $\displaystyle\forall
j\in\textsc{PO}_{k}\quad\rho_{b}(j)$
$\displaystyle=\frac{|\\{i|\widehat{\mathbf{b}}_{i}=j\text{ for
}i\in\\{1,2,\cdots,n\\}\\}|}{n/2}$

4: for $r$ in $\\{-(k-1),\ldots,(k-1)\\}$ do

5: Compute $\nu_{\ell}(r,\rho_{a},\rho_{b})$ and
$\nu_{u}(r,\rho_{a},\rho_{b})$ using Algorithms 2 and 3

6: Compute
$\varepsilon(r,\rho_{a},\rho_{b}):=\max\left\\{\left(\nu_{l}(r,\rho_{a},\rho_{b})-\frac{1}{2}\right)^{+},\left(\frac{1}{2}-\nu_{u}(r,\rho_{a},\rho_{b})\right)^{+}\right\\}$

7: end for

8:
$\widehat{m}\leftarrow\operatorname*{{\text{argmin}}}\limits_{r\in\\{-(k-1),\ldots,(k-1)\\}}\varepsilon(r,\rho_{a},\rho_{b})$;
$\epsilon\leftarrow\varepsilon(\widehat{m},\rho_{a},\rho_{b})+2k\beta$.

9: return $\widehat{m},\epsilon$.

######  Theorem 4.1.

Let $\mathbf{a,b}\in\textsc{PO}_{k}^{n}$ be the treatment and control vectors.
Define the true empirical frequency vectors as
$\eta_{a}(j)=\frac{|\\{i\in[n]|{\mathbf{a}}_{i}=j\\}|}{n}$ and
$\eta_{b}(j)=\frac{|\\{i\in[n]|{\mathbf{b}}_{i}=j\\}|}{n}$. Suppose that the
_Bernoulli Design_ is used to assign the individuals to the treatment and
control groups, i.e., each individual is assigned to one of the groups
independently with probability $1/2$. Let the observed response vectors be
$\widehat{\mathbf{a}},\widehat{\mathbf{b}}$. Then, Algorithm Median-Estimator
is a $\delta$-error median estimator $\widehat{m}$ (Definition 2.16), with
width $\epsilon\leq\epsilon^{*}(\eta_{a},\eta_{b})+2k\beta$, where
$\delta=2k\cdot\exp{\left(-2\beta^{2}n\right)}$ and runs in $O(n+k^{2})$ time.

###### Proof.

Consider the observed empirical distributions $\rho_{a}$, $\rho_{b}$ (Line 3
in Alg 1) obtained from the response vectors
$\widehat{\mathbf{a}},\widehat{\mathbf{b}}$. From Chernoff and union bound, we
have the following:

$\max_{j\in\textsc{PO}_{k}}\max\\{|\rho_{a}(j)-\eta_{a}(j)|,|\rho_{b}(j)-\eta_{b}(j)|\\}\leq\beta,$

with probability at least $1-2k\cdot\exp{\left(-2\beta^{2}n\right)}$.

Therefore, we focus on the case when the distributions $\rho_{a}$ (or
$\rho_{b}$) and $\eta_{a}$ (or $\rho_{b}$) are separated by a point-wise
distance of at most $\beta$; hence the total variation (TV distance between
$\rho_{a}$ (or $\rho_{b}$) and $\eta_{a}$ (or $\rho_{b}$) is at most
$(1/2)k\beta$. Our goal is to calculate the deviation in variability (and
width) due to this distance.

Consider the comparison of $\nu_{\ell}(r,\eta_{a},\eta_{b})$ and
$\nu_{\ell}(r,\rho_{a},\rho_{b})$ for some $r$. Suppose $\eta$ is the
optimizer of the LP in Figure 1 for $\nu_{\ell}(r,\eta_{a},\eta_{b})$. Next,
we show that there is a joint distribution $\rho$ with marginals $\rho_{a}$
and $\rho_{b}$ which is within a TV distance $k\beta$ of $\eta$. Let
$\mathcal{J}_{a}$ be an optimal coupling between $\eta_{a}$ and $\rho_{a}$,
i.e., a joint distribution with marginals $\eta_{a}$ and $\rho_{a}$ such that
$\Pr\limits_{X,X^{\prime}\leftarrow\mathcal{J}_{a}}[X\neq
X^{\prime}]\leq(1/2)k\beta$. Similarly, define $\mathcal{J}_{b}$. Then, $\rho$
is described as follows: sample $(X,Y)\leftarrow\eta$ and generate
$X^{\prime}\leftarrow\mathcal{J}_{a}|X$ and
$Y^{\prime}\leftarrow\mathcal{J}_{b}|Y$ independently, and output
$(X^{\prime},Y^{\prime})$. It is easy to see that the distribution of
$(X^{\prime},Y^{\prime})$ is $\rho$ and that
$\Pr[(X,Y)\neq(X^{\prime},Y^{\prime})]\leq k\beta$. Therefore,
$|\rho-\eta|_{TV}\leq k\beta$. However, this gives that
$\Pr\limits_{(X,Y)\leftarrow\eta}[X-Y<r]$ is
$\Pr\limits_{(X,Y)\leftarrow\eta}[X-Y<r]\pm k\beta$. Together with a symmetric
argument starting with the optimizer of $\nu_{\ell}(r,\rho_{a},\rho_{b})$, we
obtain $\nu_{\ell}(r,\eta_{a},\eta_{b})$ and
$\nu_{\ell}(r,\rho_{a},\rho_{b})\pm\beta$. Hence, we also obtain
$\epsilon^{*}(\eta_{a},\eta_{b})$ is $\epsilon^{*}(\rho_{a},\rho_{b})\pm
k\beta$.

However, the above implies that the $\epsilon$ output by Algorithm 1 is at
least $\epsilon^{*}(\eta_{a},\eta_{b})$ and at most
$\epsilon^{*}(\eta_{a},\eta_{b})+2k\beta$. Also, by Proposition 2.7 let $\eta$
be the joint distribution such that $\mathbf{(a,b)}$ is typical for $\eta$.
Then, we know that for each $r\in\textsc{ITE}$,
$q_{\ell}(r,\mathbf{a-b})=q_{\ell}(r,\eta)\leq\nu_{\ell}(r,\eta_{a},\eta_{b})\leq\nu_{\ell}(r,\rho_{a},\rho_{b})+k\beta$
and
$q_{u}(r,\mathbf{a-b})=q_{u}(r,\eta)\geq\nu_{u}(r,\eta_{a},\eta_{b})\geq\nu_{u}(r,\rho_{a},\rho_{b})-k\beta$.
In particular, this is true for $r=m$, the output of Algorithm 1. Therefore,
$m\in(\mathbf{a-b})^{0.5+\epsilon}_{0.5-\epsilon}$. This completes the proof
of correctness of Algorithm 1.

For the running time analysis observe that we calculate variability values for
$2k-1$ possible median values (Line 4 in Algorithm Median-Estimator) after
computing the empirical distributions $\rho_{a},\rho_{b}$. Therefore, our
claim regarding the running time follows from Lemmas 4.2 and 4.3. ∎

###  4.2 Computing Variability

In section 3, we have shown that computing the minimum median width of
marginals (by calculating the variability for every feasible median value)
captures the inherent complexity of approximating the median. In this section,
we show that there is a natural greedy algorithm that computes the variability
exactly. To compute variability, we need to compute
$\nu_{l}(r,\eta_{a},\eta_{b})$ and $\nu_{u}(r,\eta_{a},\eta_{b})$, where $r$
is a feasible median value and $\eta_{a},\eta_{b}$ are the marginals. As
$\nu_{l}(r,\eta_{a},\eta_{b})=\max_{\eta}\sum_{x,y}\mathbf{1}\\{x-y<r\\}\
\eta_{xy}$, we need to identify an appropriate joint distribution $\eta$, that
maximizes the total probability for the lower quantile. It turns out that we
can formulate this as a linear program shown below (see Figure 1), with
marginal and non-negativity constraints.

$\displaystyle\underline{\text{Objective}}:$
$\displaystyle\quad\max\sum_{x,y}\mathbf{1}\\{x-y<r\\}\eta_{xy}$ (Marginal
Constraints)
$\displaystyle\quad\sum_{y\in\textsc{PO}_{k}}\eta_{xy}=\eta_{a}(x)\quad\forall
x\in\textsc{PO}_{k}$ (Marginal Constraints)
$\displaystyle\quad\sum_{x\in\textsc{PO}_{k}}\eta_{xy}=\eta_{b}(y)\quad\forall
y\in\textsc{PO}_{k}$ (Non-negativity) $\displaystyle\quad\eta_{xy}\geq
0\quad\forall x,y\in\textsc{PO}_{k}\times\textsc{PO}_{k}$
Figure 1: Linear Program for computing the lower quantile component of
variability Algorithm 2 Compute-Variability-Lower-Quantile

1: Input: Marginal distributions $\eta_{a},\eta_{b}$ and feasible median value
$r\in\\{-(k-1),\cdots,(k-1)\\}$.

2: Output: $\nu_{l}(r,\eta_{a},\eta_{b})$.

3: $\forall x\in\textsc{PO}_{k},y\in\textsc{PO}_{k}$, set $\eta_{xy}\leftarrow
0$

4: for $y\in\\{k-1,k-2,\cdots,0\\}$ do

5: for $x\in\\{y+r-1,y+r-2,\cdots,0\\}\cap\textsc{PO}_{k}$ do

6: Identify all the marginal constraints in LP 1 containing $\eta_{xy}$ and
then slowly increase the value of $\eta_{xy}$ until one of them becomes tight.

7: We mark the corresponding column or row with the tight constraint as
_frozen_ and do not update any values in them.

8: end for

9: end for

10: Return
$\nu_{l}(r,\eta_{a},\eta_{b})\leftarrow\sum_{x,y\in\textsc{PO}_{k}\times\textsc{PO}_{k}}\mathbf{1}\\{x-y<r\\}\eta_{xy}$.

We can view solving the linear program (in Figure 1) as essentially filling in
the entries of a 2-dimensional matrix of size
$\textsc{PO}_{k}\times\textsc{PO}_{k}$ such that the difference in indices is
strictly smaller than $r$, i.e., $x-y<r\ \forall x,y$. To do so, we employ a
greedy strategy in Algorithm Compute-Variability-Lower-Quantile (see Alg 2).
Without loss of generality, we fill the matrix, greedily from the last column
and last row and move up the column (see Lines 4-5 in Alg 2). Our greedy
strategy involves increasing the value of $\eta_{xy}$ until one of the
marginal constraints is tight and freezing the values of the row/column
corresponding to the constraint (see Lines 7-10 in Alg 2). A similar approach
results in an algorithm 3 for computing the upper quantile component of
variability. Formally:

######  Lemma 4.2.

Given a feasible value $r\in\\{-(k-1),\cdots,(k-1)\\}$ and marginals
$\eta_{a},\eta_{b}$, there is an algorithm for computing
$\nu_{l}(r,\eta_{a},\eta_{b})$ exactly, and runs in $O(k)$ time.

###### Proof.

Let $\eta^{*}$ denote the output of the greedy Algorithm and $\bar{\eta}$ be
an optimizer of Figure 1. Let us assume:

$\sum_{x,y}\mathbf{1}\\{x-y<r\\}\eta^{*}_{xy}<\sum_{x,y}\mathbf{1}\\{x-y<r\\}\bar{\eta}_{xy}$

As the greedy algorithm fills the entries of $\eta^{*}_{xy}$ by making sure
that one of the constraints is tight: we will be able to identify a pair (wlog
first pair) $x,y$ such that $\eta^{*}_{xy}>\bar{\eta}_{xy}$. Otherwise, we
would have $\bar{\eta}=\eta^{*}$. Now consider the marginal constraint

$\sum_{x\in\textsc{PO}_{k}}\eta^{*}_{xy}=\sum_{x\in\textsc{PO}_{k}}\bar{\eta}_{xy}=\eta_{b}(y).$

Therefore, we have a $u$ such that: $\eta^{*}_{uy}<\bar{\eta}_{uy}$. Similarly
the marginal constraint

$\sum_{y\in\textsc{PO}_{k}}\eta^{*}_{xy}=\sum_{y\in\textsc{PO}_{k}}\bar{\eta}_{xy}=\eta_{a}(x).$

Therefore, we have a $v$ such that: $\eta^{*}_{xv}<\bar{\eta}_{xv}$. We update
$\bar{\eta}$ by increasing the value of $\bar{\eta}_{xy}$ as follows:

  1. 1.

$\bar{\eta}^{\text{upd}}_{uy}=\max\\{0,\bar{\eta}_{uy}-(\eta^{*}_{xy}-\bar{\eta}_{xy})\\}$

  2. 2.

$\bar{\eta}^{\text{upd}}_{xv}=\max\\{0,\bar{\eta}_{xv}-(\eta^{*}_{xy}-\bar{\eta}_{xy})\\}$

  3. 3.

$\bar{\eta}^{\text{upd}}_{uv}=\bar{\eta}_{uv}+(\eta^{*}_{xy}-\bar{\eta}_{xy})$

  4. 4.

$\bar{\eta}^{\text{upd}}_{xy}=\bar{\eta}_{xy}+(\eta^{*}_{xy}-\bar{\eta}_{xy})$

Observe that such an update satisfies the marginal constraints and the
objective function remains the same for both $\bar{\eta}$ and
$\bar{\eta}^{\text{upd}}$. The updated $\bar{\eta}^{\text{upd}}$ such that
$\bar{\eta}^{\text{upd}}_{xy}=\eta^{*}_{xy}$. By repeatedly following this
updating process, we will have $\bar{\eta}^{\text{upd}}=\eta^{*}$, a
contradiction to our claim. As we eliminate a row or column by making a
marginal constraint tight in each iteration of our algorithm. The total number
of iterations and the values of $(x,y)\in\textsc{PO}_{k}\times\textsc{PO}_{k}$
pairs that we fill are at most $2k$. Hence, the running time of our algorithm
is $O(k)$. ∎

As $\nu_{u}(r,\eta_{a},\eta_{b})=\min_{\eta}\sum_{x,y}\mathbf{1}\\{x-y\leq
r\\}\ \eta_{xy}$, we need to identify an appropriate joint distribution
$\eta$, that minimizes the total probability for the upper quantile. It turns
out that we can formulate this as a linear program shown below (see Figure 2),
with marginal and non-negativity constraints. Similar to Algorithm 2, we have
Algorithm Compute-Variability-Upper-Quantile that computes the value
$\nu_{u}(r,\eta_{a},\eta_{b})$ exactly.

$\displaystyle\underline{\text{Objective}}:$
$\displaystyle\quad\min\sum_{x,y}\mathbf{1}\\{x-y\leq r\\}\eta_{xy}$
$\displaystyle\quad\sum_{y\in\textsc{PO}_{k}}\eta_{xy}=\eta_{a}(x)\quad\forall
x\in\textsc{PO}_{k}$
$\displaystyle\quad\sum_{x\in\textsc{PO}_{k}}\eta_{xy}=\eta_{b}(y)\quad\forall
y\in\textsc{PO}_{k}$ $\displaystyle\quad\eta_{xy}\geq 0\quad\forall
x,y\in\textsc{PO}_{k}\times\textsc{PO}_{k}$
Figure 2: Linear Program for computing the upper quantile component of
variability Algorithm 3 Compute-Variability-Upper-Quantile

1: Input: Marginal distributions $\eta_{a},\eta_{b}$ and feasible median value
$r\in\\{-(k-1),\cdots,(k-1)\\}$.

2: Output: $\nu_{u}(r,\eta_{a},\eta_{b})$.

3: $\forall x\in\textsc{PO}_{k},y\in\textsc{PO}_{k}$, set $\eta_{xy}\leftarrow
0$

4: for $y\in\\{k-1,k-2,\cdots,0\\}$ do

5: for $x\in\\{y+r,y+r-1,\cdots,0\\}$ do

6: Identify all the marginal constraints in LP 2 containing $\eta_{xy}$ and
then slowly increase the value of $\eta_{xy}$ until one of them becomes tight.

7: We mark the corresponding column or row with the tight constraint as
_frozen_ and do not update any values in them.

8: end for

9: end for

10: Return
$\nu_{l}(r,\eta_{a},\eta_{b})\leftarrow\sum_{x,y\in\textsc{PO}_{k}\times\textsc{PO}_{k}}\mathbf{1}\\{x-y\leq
r\\}\eta_{xy}$.

######  Lemma 4.3.

Given a feasible median value $r\in\\{-(k-1),\cdots,(k-1)\\}$ and marginals
$\eta_{a},\eta_{b}$, Algorithm 3 computes the value
$\nu_{u}(r,\eta_{a},\eta_{b})$ exactly and runs in time $O(k)$.

###### Proof.

A similar proof as that of lemma 4.2. ∎

##  5 Conclusion

In this work, we studied the task of median treatment effect estimation, to
capture the heterogeneity in the data. We argued that the task is inestimable,
provided new notions of approximations, and showed optimal algorithms with
minimum approximation error. We want to highlight that our results make no
distributional assumptions, and can be extended for any quantile $q<0.5$ (in
addition to median), by modifying the definition of variability (definition
2.8) appropriately. For future work, a potential direction is to study the
settings where the potential outcomes are continuous.

## Acknowledgments

The authors would like to thank the following individuals for many helpful
discussions: Anup Rao, Kirankumar Shiragur, Christopher Harshaw, Shiva
Kasiviswanathan, and Avi Feller. The authors also would like to thank the
organizers for the invitation to the Causality program at the Simons Institute
for the Theory of Computing, which helped facilitate early discussions of this
work.

## References

  * Addanki et al. [2022] Raghavendra Addanki, David Arbour, Tung Mai, Cameron Musco, and Anup Rao.  Sample constrained treatment effect estimation.  _Advances in Neural Information Processing Systems_ , 35:5417–5430, 2022.
  * Athey and Imbens [2016] Susan Athey and Guido Imbens.  Recursive partitioning for heterogeneous causal effects.  _Proceedings of the National Academy of Sciences_ , 113(27):7353–7360, 2016.
  * Bojinov et al. [2020] Iavor Bojinov, Guillaume Saint-Jacques, and Martin Tingley.  Avoid the pitfalls of a/b testing make sure your experiments recognize customers’ varying needs.  _Harvard Business Review_ , 98(2):48–53, 2020\.
  * Chernozhukov and Hansen [2005] Victor Chernozhukov and Christian Hansen.  An iv model of quantile treatment effects.  _Econometrica_ , 73(1):245–261, 2005.
  * Chernozhukov et al. [2018] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins.  Double/debiased machine learning for treatment and structural parameters, 2018.
  * Ding et al. [2019] Peng Ding, Avi Feller, and Luke Miratrix.  Decomposing treatment effect variation.  _Journal of the American Statistical Association_ , 114(525):304–317, 2019.
  * Harshaw et al. [2023] Christopher Harshaw, Fredrik Sävje, Daniel A Spielman, and Peng Zhang.  Balancing covariates in randomized experiments with the gram–schmidt walk design.  _Journal of the American Statistical Association_ , pages 1–27, 2023\.
  * Heckman et al. [1997] James J Heckman, Jeffrey Smith, and Nancy Clements.  Making the most out of programme evaluations and social experiments: Accounting for heterogeneity in programme impacts.  _The Review of Economic Studies_ , 64(4):487–535, 1997.
  * Howard and Ramdas [2022] Steven R Howard and Aaditya Ramdas.  Sequential estimation of quantiles with applications to a/b testing and best-arm identification.  _Bernoulli_ , 28(3):1704–1728, 2022.
  * Imai and Ratkovic [2013] Kosuke Imai and Marc Ratkovic.  Estimating treatment effect heterogeneity in randomized program evaluation.  _Annals of Applied Statistics_ , 7(1):443–470, 2013.
  * Kallus [2022] Nathan Kallus.  What’s the harm? sharp bounds on the fraction negatively affected by treatment.  _Advances in Neural Information Processing Systems_ , 35:15996–16009, 2022.
  * Kennedy [2023] Edward H Kennedy.  Towards optimal doubly robust estimation of heterogeneous causal effects.  _Electronic Journal of Statistics_ , 17(2):3008–3049, 2023.
  * Kennedy et al. [2017] Edward H Kennedy, Zongming Ma, Matthew D McHugh, and Dylan S Small.  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, 2017.
  * Leqi and Kennedy [2021] Liu Leqi and Edward H Kennedy.  Median optimal treatment regimes.  _arXiv preprint arXiv:2103.01802_ , 2021.
  * Martinez-Taboada and Kennedy [2023] Diego Martinez-Taboada and Edward H Kennedy.  Counterfactual density estimation using kernel stein discrepancies.  _arXiv preprint arXiv:2309.16129_ , 2023.
  * Pouget-Abadie et al. [2019] Jean Pouget-Abadie, Kevin Aydin, Warren Schudy, Kay Brodersen, and Vahab Mirrokni.  Variance reduction in bipartite experiments through correlation clustering.  _Advances in Neural Information Processing Systems_ , 32, 2019.