MS-0001-1922.65

Wang, Chen, and Wang

Contextual Optimization under Covariate Shift

Contextual Optimization under Covariate Shift: A Robust Approach by
Intersecting Wasserstein Balls

Tianyu Wang Department of Industrial Engineering and Operations Research,
Columbia University, New York, NY 10027<EMAIL_ADDRESS>Ningyuan
Chen Department of Management, University of Toronto Mississauga,
Rotman School of Management, University of Toronto, Canada,
<EMAIL_ADDRESS>Chun Wang School of Economics and Management,
Tsinghua University, Beijing, China<EMAIL_ADDRESS>
In contextual optimization, a decision-maker observes historical samples of
uncertain variables and associated concurrent covariates, without knowing
their joint distribution. Given an additional covariate observation, the goal
is to choose a decision that minimizes some operational costs. A prevalent
issue here is covariate shift, where the marginal distribution of the new
covariate differs from historical samples, leading to decision performance
variations with nonparametric or parametric estimators. To address this, we
propose a distributionally robust approach that uses an ambiguity set by the
intersection of two Wasserstein balls, each centered on typical nonparametric
or parametric distribution estimators. Computationally, we establish the
tractable reformulation of this distributionally robust optimization problem.
Statistically, we provide guarantees for our Wasserstein ball intersection
approach under covariate shift by analyzing the measure concentration of the
estimators. Furthermore, to reduce computational complexity, we employ a
surrogate objective that maintains similar generalization guarantees. Through
synthetic and empirical case studies on income prediction and portfolio
optimization, we demonstrate the strong empirical performance of our proposed
models.111Our previous draft (Wang et al. 2021) has a different title
“Distributionally Robust Prescriptive Analytics with Wasserstein Distance”,
which focuses on the measure concentration results of the kernel estimator
under Wasserstein distance without covariate shifts. It constitutes Section
4.1 of this paper.

distributionally robust optimization, Wasserstein distance, contextual
optimization, ambiguity intersection, covariate shift

##  1 Introduction

Most decision-making problems involve uncertainty. In an increasingly data-
rich world, the decision-maker may rely on covariates or contextual
information to predict the uncertain information first when making such a
data-driven decision. This is the central question in the so-called
_contextual (stochastic) optimization_ , which has attracted the attention of
many scholars.

##### Contextual Optimization.

Consider the following optimization problem:

$\min_{z\in\mathcal{Z}}\mathbb{E}[c(z,Y)|X=x].$ (1)

Here $z$ denotes the decision vector taken from the feasible region
$\mathcal{Z}\subseteq\mathbb{R}^{d_{z}}$. The objective function form $c(z,Y)$
is known and depends on the random vector $Y\in\mathbb{R}^{d_{y}}$, which
encodes the random environmental information and typically unobserved before
making decisions. Usually, the decision-maker has access to a covariate
$X\in\mathbb{R}^{d_{x}}$ without knowing the conditional distribution $Y|X=x$.
Besides, the decision-maker also observes $n$ samples
$\\{(x_{i},y_{i})\\}_{i=1}^{n}$ instead of the joint distribution $(X,Y)$,
presumably from the decisions made in the past. The goal of contextual
optimization is to choose a reasonable decision $z$ based on the current
covariate $X=x$ and historical data $\\{(x_{i},y_{i})\\}_{i=1}^{n}$ to
minimize the objective (1). Contextual optimization has been studied under
different settings (Ban and Rudin 2019, Bertsimas and Kallus 2020, Bertsimas
and Koduri 2022, Elmachtoub and Grigas 2022, Qi et al. 2021, Kallus and Mao
2023). Theoretical generalization bounds of specific algorithms are
established in El Balghiti et al. (2019), Hu et al. (2022). This formulation
has a wide range of applications; see Sadana et al. (2024) for a review.

Both algorithmically and theoretically, the literature has been focused on a
stationary environment for the covariate $X$. In other words, the distribution
of $X$ in the historical data $\\{(x_{i},y_{i})\\}_{i\in[n]}$ coincides with
the that of the test environment $X=x$. When this assumption fails, we refer
to the phenomenon as “covariate shift”. Needless to say, covariate shift will
bring challenge to the design of algorithms and the performance analysis. In
practice, however, covariate shift is prevalent.

##### Motivating Example.

We illustrate the idea of covariate shift through the following example.
Suppose we want to allocate a limited resource among individuals to maximize
the total social welfare, which is prevalent in contextual optimization
literature (Wilder et al. 2019). We have already done some surveys to extract
some data $\\{(x_{i},y_{i})\\}_{i=1}^{n}$ with individuals’ features $X$ and
their social outcome $Y$. However, the targeted population of individuals
often changes when we deploy our resource allocation systems in practice since
we have some unseen individuals. The survey data we have may not be
representative of all individuals (e.g. younger individuals may be more
inclined to participate in the survey), incurring a so-called sampling bias.
We may also collect data from one state and test the trained model in other
states. In both cases, the feature distribution differs between the historical
data and the new environment while the conditional distribution $Y|X$ is still
the same across the two environments. We want to understand how models would
perform under such misalignment and how to use data from historical data to
construct robust models for a different domain.

In this paper, we study a general data-driven distributionally robust
contextual optimization problem (1) and propose an algorithmic framework to
mitigate the negative impact when the new covariate $X=x$ is potentially far
from the region well covered by the historical data. More precisely, we build
a distributionally robust optimization (DRO) model that sets the ambiguity set
as the intersection of two Wasserstein balls, constructed using the parametric
and nonparametric estimators from the historical data, respectively. To see
the intuition why the intersection of Wasserstein balls can help address the
issue brought by covariate shift, note that nonparametric estimators such as
the Kernel methods usually perform well when the covariate to be predicted is
well covered by the historical data because it does not suffer from model
misspecification. On the other hand, they may falter and become unstable when
the covariate is drawn from a different distribution than the training data,
sometimes coined as “out-of-distribution” in the machine learning literature.
To the contrary, parametric estimators such as linear regression can
generalize robustly to shifted covariates but may suffer from model
misspecification. By formulating the DRO using the intersection, we can obtain
a solution that strikes a balance and enjoys the best of both worlds. We
summarize our main contribution below.

$\bullet$ Modeling framework of contextual optimization under covariate
shifts. We apply a Wasserstein-DRO framework for general contextual
optimization problems under distribution shifts. This framework is flexible to
implement with mild technical assumptions. Especially, we ensemble the power
of two learning procedures through the intersection of two Wasserstein
ambiguity sets, where both centers and their ball sizes are adapted to the new
covariate $x$. This general ambiguity set design is new to the DRO and
contextual optimization community with practical interests.

$\bullet$ Computational reformulation and interpretation. We provide a
tractable reformulation to the formulated DRO problem. The tractability of the
reformulation is a highly nontrivial extension to the classic Wasserstein DRO
formulation since it involves the intersection of Wasserstein balls, which may
be infeasible and necessitates another layer of duality. Our formulation can
be of independent interests on the computational side. However, such
reformulation leads to an optimization problem with too many decision
variables. Therefore, we also construct a surrogate objective driven by
another ambiguity to reduce such computation complexity while remaining
provable performance guarantees.

$\bullet~{}$Measure concentration and statistical guarantees. We provide
measure concentration results under the $p$-Wasserstein distance for
parametric and nonparametric estimators in the environment with covariate
shifts. In particular, we provide the finite-sample measure concentration
result under the $p$-Wasserstein distance for two important classes of
nonparametric estimators: the Nadaraya-Watson kernel and kNN estimators. To
our knowledge, this result has not been derived even without covariate shifts.
The bound further demonstrates the need to combine nonparametric and
parametric estimators in a flexible manner, as it can be better or worse than
that of the parametric estimator depending on the magnitude of covariate
shifts. Based on the measure concentration, we establish out-of-sample
performance guarantees for the solutions obtained from the DRO formulation
when the ambiguity set is constructed using the intersection.

###  1.1 Related Works

##### Optimization under Uncertainty.

In data-driven stochastic optimization, sample average approximation (SAA)
(Shapiro et al. 2014) and distributionally robust optimization (DRO) (Delage
and Ye 2010) are both popular approaches. In particular, DRO gets the
attention of many scholars in recent years due to its ability to provide
robust solutions against a class of potentially worst-case scenarios of the
unknown parameter. The class of scenarios can be considered as a set of
distributions close to the historical data under some distance metrics. This
includes the Wasserstein metric (Esfahani and Kuhn 2018, Gao and Kleywegt
2023), goodness-of-fit test (Bertsimas et al. 2018) and $\phi$-divergence
(Ben-Tal et al. 2013, Duchi et al. 2021). We base our DRO formulation under
the umbrella of Wasserstein-DRO, which has been explored extensively in recent
years throughout literature for its statistical guarantees (Blanchet et al.
2019, Gao et al. 2022, Gao and Kleywegt 2023), computational benefits
(Esfahani and Kuhn 2018, Blanchet and Murthy 2019) and real-world applications
(Shafieezadeh Abadeh et al. 2015). We refer readers to a more comprehensive
discussion in Rahimian and Mehrotra (2019).

The stationary assumption of the historical data and the test environment,
i.e., no covariate shifts, does not hold in many practical problems.
Theoretically, DRO is shown to maintain out-of-sample performance under
limited data as long as the ambiguity set is comprehensive (Lee and Raginsky
2017). However, such a guarantee hinges on a large enough ambiguity set that
takes into account potential distributional shifts, which can be overly
conservative in practice. Indeed, empirical studies have shown that a
conservative ambiguity set can lead to performance degradation (Hu et al.
2018, Liu et al. 2024) that offsets the potential benefit of protection
against covariate shifts. A formulation designed specifically for covariate
shifts in this paper offers an attractive alternative.

##### Decision-making under Distribution Shifts.

Some previous studies (Besbes et al. 2022, Bennouna and Van Parys 2022, Gupta
and Kallus 2022, Besbes et al. 2023) provide statistical guarantees for
existing models or build new models when the historical data is not sampled
i.i.d. Their setups are different from ours. Besbes et al. (2022), Gupta and
Kallus (2022), Besbes et al. (2023) investigate the rate-optimal model
performance when each sample (or some samples in one group) is generated from
a different distribution. Bennouna and Van Parys (2022) consider a DRO model
under three types of data errors including statistical estimation, finite
precision, and misspecification. In contrast, we are interested in the
distribution shifts setting (Quinonero-Candela et al. 2008), or more
specifically, the covariate shift (Shimodaira 2000) between the historical
data and the test environment. That is, the historical data is sampled i.i.d.
while the marginal distribution of covariates in the test environment shifts.
There are other types of distribution shifts studied in the contextual
optimization; for example, Johnson-Yu et al. (2023) study an end-to-end robust
method under _label shifts_ (shift in the distribution of $Y$). Because we
focus on contextual optimization, standard approaches to handle covariate
shifts such as importance weight or distribution matching (Sugiyama et al.
2007, Quinonero-Candela et al. 2008) do not apply: they require a test dataset
of a reasonable size whereas we only have one covariate sample in the test
environment to make decisions in contextual optimization.

Recent advances in DRO allow the decision-maker to incorporate richer
covariate information into the optimization problem. Different
distributionally robust frameworks have been proposed including Esteban-Pérez
and Morales (2022), Bertsimas and Van Parys (2021), Yang et al. (2022), Kannan
et al. (2023) with applications to portfolio optimization (Nguyen et al. 2021,
Pun et al. 2023) and newsvendor (Zhang et al. 2024). However, these approaches
do not assume any structure between the historical data and the test
environment. Therefore, these models only provide coverage guarantees under
worst-case distribution shifts characterized by the statistical distances,
which are sometimes too conservative in the absence of other shifts than the
covariate. Duchi et al. (2023) proposed an ambiguity set specialized for
covariate shifts but allow any arbitrary covariate shift for the test domain
within the DRO model since they do not focus on the contextual optimization
setup and not observe any new covariates.

##### Statistical Tradeoff in Estimating Distributions.

Our approach belongs to an estimate-then-optimize procedure in contextual
optimizations (Bertsimas and Kallus 2020) where we first estimate the
conditional distribution and use that estimated distribution to make
subsequent decisions. Parametric and nonparametric procedures are two families
of estimators commonly used to estimate the conditional distribution from the
historical data. Their fundamental tradeoff has been characterized in the
literature (Hastie et al. 2009, Györfi et al. 2006). Parametric estimators
enjoy a better convergence rate but require stronger assumptions on the
distribution of $(X,Y)$ to ensure consistency, as they assume a specific
parametric form for the distribution. These methods include linear models or
fixed-width neural networks. In contrast, nonparametric estimators typically
relax distribution assumptions and are shown to be minimax optimal in some
cases (Györfi et al. 2006, Wainwright 2019). These methods include kernel,
nearest neighbor, and tree-based approaches. However, nonparametric approaches
are more adversely affected by covariate shifts, especially when there are few
samples around the new covariate $X=x$, making them less favorable compared to
parametric procedure in this case. In this study, we propose to combine both
types of estimators by intersecting their respective Wasserstein balls so that
the downstream decision may enjoy the best of both worlds.

To combine the power of both parametric and nonparametric models, the so-
called _ensemble learning_ procedure can capture diverse aspects of the data
and provide more robust predictions under distributional changes (Jacobs et
al. 1991, Dietterich 2000, Lakshminarayanan et al. 2017). This refers to
bagging (Breiman 2001) (or boosting (Freund and Schapire 1997)) approaches by
averaging (or augmenting) several estimators $\hat{\mathbb{E}}[Y|x]$ to
approximate $\mathbb{E}[Y|x]$ (Hastie et al. 2009). These methods are built to
output a _single_ more robust estimate $\mathbb{E}[Y|X=x]$ or sometimes
$\mu_{Y|x}$. In the contextual optimization, the focus is on the downstream
decision $z(x)$ instead of the intermediate distribution. Taskesen et al.
(2021) develops a combination of both source and target data for transfer
learning which only requires moment construction for their Gaussian
distribution on the non-contextual side. Compared to their formulation, we
consider the intersection of two distance-based balls instead of the moment-
based balls. Except for that work, to the best of our knowledge, we are the
first to put the ensemble learning framework in the distributionally robust
decision-making setup through the intersection of two distribution estimators.

###  1.2 Organization

We introduce the problem setup and formulate our proposed ambiguity
intersection in Section 2. Under this constructed ambiguity set, we show that
the min-max optimization problem can be reformulated as convex optimization
problems to ensure computational tractability in Section 3 as long as the
original objective is convex. After that, we study the measure concentration
results of the ambiguity set, specifically focusing on both estimators with an
illustration of the magnitude of covariate shifts in Sections 4.1 and 4.2 to
discuss its statistical properties. Section 5 provides a computational
surrogate of this optimization problem with statistical performance guarantees
and compare the performance of this surrogate optimization with the original
model. Section 6 provide some discussions on such problem setup. Section 7
presents numerical results on both synthetic and real data examples of income
prediction and CVaR-based portfolio optimization to demonstrate the benefit of
our formulation. Additional discussions and experimental results, and proofs
of all theoretical results are deferred to the Appendix.

##### Technical notations

Throughout the paper, we use upper case letters to represent random quantities
and lower case letters to represent deterministic quantities or samples. We
denote $[n]\coloneqq\\{1,2,...,n\\}$, $\|\cdot\|_{p}$ as the $\ell_{p}$-norm.
We use $\delta_{x}$ to represent Dirac measure on the single point $x$ and
$\mathbb{I}_{\left\\{\cdot\right\\}}$ denote the indicator function. Denote
$\sigma_{\mathcal{Y}}(\cdot)$ as the support function of the domain
$\mathcal{Y}$, i.e.,
$\sigma_{\mathcal{Y}}(\xi)=\sup_{y\in\mathcal{Y}}\\{\xi^{T}y\\}$. Denote
$f^{*}(\cdot)$ as the convex conjugate of $f(\cdot)$ and $\|z\|_{*}$ denote
the dual norm of $z$. For a set $\mathcal{X}$, denote $\mu(\mathcal{X})$ as
the family of probability distributions on $\mathcal{X}$. We use
$\mu_{Y|x},\nu_{Y|x}$ to denote the conditional probability measure of $Y$
given $X=x$; $\hat{\mu}_{Y}$ for a generic measure of $Y$ based on the
historical samples, and specified to the kernel version $\mu_{n,x}$ in (11).
For the random sequences $\alpha_{n},\beta_{n}$, we write
$\alpha_{n}=o_{p}(\beta_{n})$ if $\alpha_{n}/\beta_{n}$ converges to zero in
probability or $\alpha_{n}=O_{p}(\beta_{n})$ if $\alpha_{n}/\beta_{n}$ is
uniformly bounded in probability. We use
$\mathcal{N}(\mu,\sigma^{2}),U(\cdot),\exp(\cdot)$ as the normal, uniform, and
exponential distribution. We use $\mathcal{P}_{p}(\mathcal{X})$ to denote the
set of all distributions with finite $p$-th moment in the domain
$\mathcal{X}$.

##  2 Setups

Recall that we are interested in solving the contextual (stochastic)
optimization problem (1) after observing $X=x$. We use two measures, $\mu$ and
$\nu$ (abbreviates for $\mu_{(X,Y)}$ and $\nu_{(X,Y)}$), to distinguish
between the distributions of $(X,Y)$ in the training (historical) data and the
test data, respectively. In particular, denote the marginal measure of $X$ and
the measure $Y$ given $X=x$ as $\mu_{X},\mu_{Y|x}$ and $\nu_{X},\nu_{Y|x}$ in
the two distributions. For any $x\in\mathcal{X}$, we do not know the true
conditional distribution, on which the random cost $c(z;Y)$ depends, but only
observe $\\{(x_{i},y_{i})\\}_{i=1}^{n}$ drawn from $\mu$ independently. The
objective is to minimize the cost $\mathbb{E}_{\nu_{Y|x}}[c(z(x);Y)]$ under
the test distribution. That is, we want to find the optimal solution to:

$J^{*}\coloneqq\inf_{z\in\mathcal{Z}}\mathbb{E}_{\nu_{Y|x}}[c(z,Y)],$ (2)

for all $x\in\mathcal{X}$.

We adopt the framework of distributionally robust optimization with
Wasserstein distance to tackle the problem. Based on the historical samples
$\left\\{(x_{i},y_{i})\right\\}_{i=1}^{n}$ and $X=x$, we obtain our decision
through the following objective:

$\hat{J}_{n}\coloneqq\inf_{z\in\mathcal{Z}}\sup_{\hat{\nu}_{Y|x}\in\mathcal{A}}\mathbb{E}_{\hat{\nu}_{Y|x}}[c(z,Y)],$
(3)

where the ambiguity set is defined as:

$\mathcal{A}=\\{\mu:\mathcal{W}_{p}(\mu,\hat{\mu}_{1})\leq\varepsilon_{1},\mathcal{W}_{p}(\mu,\hat{\mu}_{2})\leq\varepsilon_{2}\\}.$
(4)

Here $\mathcal{W}_{p}(\cdot,\cdot)$ is the $p$-Wasserstein distance:

######  Definition 2.1 (Wasserstein distance)

The type-$p$ $(1\leq p<+\infty)$ Wasserstein distance between two
distributions $\mu$ and $\nu$ supported on $\Xi$ is defined as:

$\mathcal{W}_{p}(\mu,\nu)=\inf_{\Pi\in\mathcal{M}{(\Xi\times\Xi)}}\left\\{\left(\int_{\Xi\times\Xi}\|x-y\|_{1}^{p}\Pi(dx,dy)\right)^{\frac{1}{p}}:\Pi_{x}=\mu,\Pi_{y}=\nu\right\\},$

where $\Pi_{x}$ and $\Pi_{y}$ are the marginal distribution of $\Pi$.

At this point, we keep the two reference measures $\hat{\mu}_{1}$ and
$\hat{\mu}_{2}$ generic with the following representation assumptions:
[Structure of Distributions] $\hat{\mu}_{1},\hat{\mu}_{2}$ in (3) are two
discrete distributions with their formulations being
$\hat{\mu}_{1}=\sum_{i=1}^{n}w_{1,i}\delta_{y_{1,i}}$,
$\hat{\mu}_{2}=\sum_{j=1}^{n}w_{2,j}\delta_{y_{2,j}}$ for $n$ points
$\\{y_{1,i}\\}_{i\in[n]}$ and $\\{y_{2,j}\\}_{j\in[n]}$, $w_{1,i}\geq
0,w_{2,j}\geq 0,\forall i,j\in[n]$ and
$\sum_{i\in[n]}w_{1,i}=\sum_{j\in[n]}w_{2,j}=1$. Definition 2.1 is natural in
literature since the standard reference measures constructed from data are
typically discrete (Ban and Rudin 2019, Bertsimas and Kallus 2020, Bertsimas
and Van Parys 2021). Furthermore, as we will specify in Section 4, our
proposed $\hat{\mu}_{1},\hat{\mu}_{2}$ also satisfy this assumption.

The form of the ambiguity set deviates largely from the literature, which
typically uses a Wasserstein ball centered at a single reference measure
($\\{\mu:\mathcal{W}_{p}(\mu,\hat{\mu}_{1})\leq\varepsilon\\}$) as the
ambiguity set. Instead, we use the intersection of two Wasserstein balls. The
intersection of two balls ensures that we only consider the estimated
conditional distribution $\hat{\nu}_{Y|x}$ to be close to both reference
measures. This usually leads to a more robust result when the two reference
measures, usually estimated from the training data, are largely complementary
and may outperform the other in different situations. Our method is different
from standard ensemble learning procedures in statistical learning, where the
goal of these methods is a quantity that closes to $\mathbb{E}[Y|x]$ or
$\nu_{Y|x}$ and does not consider the downstream optimization problem.

When the test distribution is identical to the training distribution, i.e.,
$\mu\equiv\nu$, the literature (Bertsimas and Van Parys 2021, Ban and Rudin
2019) has shown that using an ambiguity set centered at a nonparametric
estimator achieves good performance, especially when the actual form of
$\nu_{Y|x}$ is unknown and may be misspecified by parametric estimators.
However, in this paper, we focus on a setting where the distribution of the
covariate $X$ in the test data may shift from the training data:

######  Definition 2.2 (Covariate Shift)

We say there is a covariate shift between $\mu_{(X,Y)}$ and $\nu_{(X,Y)}$ if
$\nu_{X}\neq\mu_{X}$ but $\nu_{Y|x}=\mu_{Y|x}$ for all $x\in\mathcal{X}$.

Clearly, the covariate shift describes a shift in the marginal distribution of
the covariate $X$ but not the conditional distribution of $Y$ given $X$. As
the following example shows, in this case, nonparametric estimators may suffer
from a limited number of local samples and a large variance:

###### Example 2.3

Consider a mean prediction task where the underlying true relationship is
$Y|x\sim N(2x+0.3\sin(3\pi x/2),0.1)$. The training data
$\\{x_{i}\\}_{i\in[n]}$ consists of 100 (or 500) sample points drawn from
$U(0,1)$. The test set is from $[1,1.4]$ and out of the domain of the training
set. As illustrated in Figure 1, when the test point $x>1$, a large covariate
shift occurs since we have few training samples there and the nonparametric
regression estimator suffers high variance in that area.

(a) $n=100$

(b) $n=500$

Figure 1: Larger variance of the nonparametric estimator (Nadaraya-Watson
regression) versus the parametric estimator (linear regression)

Summarizing the insights from the previous literature and the example, we can
see that using a single parametric or nonparametric estimator as the reference
measure in the ambiguity set may lead to poor performance, due to model
misspecification or covariate shift.

When we consider the ambiguity set (4), we are essentially constructing a
“joint estimator” for the conditional distribution. This approach ensures that
the performance guarantee of the decision is upper bounded by the radius of
the better Wasserstein ball while being less affected by the inferior
Wasserstein ball, which typically has a larger radius. Meanwhile, intersecting
with the inferior ball can eliminate some peculiar distributions present in
the better ball while keeping the true distribution inside, thus producing a
more reasonable ambiguity set, as demonstrated in Section 6.1. Therefore this
intersection can achieve the “best of the two worlds”, with and without
covariate shifts. Furthermore, due to the computational complexity this “joint
estimator” suffers from, we suggest a computationally efficient surrogate,
referred to as the “interpolate estimator”, which offers comparable
performance guarantees. We summarize the roadmap of our subsequent analysis in
Figure 2. We denote the optimization model under the “joint estimator” as
2W-DRO and “interpolate estimator” as IW-DRO respectively.

Figure 2: Roadmap of Theoretical Guarantees

##  3 2W-DRO: Feasibility and Tractable Reformulation

In this section, we focus on the computational part, i.e., problem
reformulation in (3). In practice, the intersection of the two Wasserstein
balls may be empty if $\varepsilon_{1}$ and $\varepsilon_{2}$ are not chosen
wisely. We provide the following sufficient and necessary conditions with a
computationally tractable certificate to indicate when $\mathcal{A}$ is empty
or not.

######  Proposition 3.1 (Feasibility of 2W-DRO)

The ambiguity set
$\mathcal{A}=\\{\mu:\mathcal{W}_{p}(\mu,\hat{\mu}_{1})\leq\varepsilon_{1},\mathcal{W}_{p}(\mu,\hat{\mu}_{2})\leq\varepsilon_{2}\\}$
is not empty if and only if
$\mathcal{W}_{p}(\hat{\mu}_{1},\hat{\mu}_{2})\leq\varepsilon_{1}+\varepsilon_{2}$.

This result follows the same argument as Proposition 5.3 in Taskesen et al.
(2021) by noticing that $(\mathcal{Y},\mathcal{W}_{p})$ is a geodesic metric
space (Villani 2009). Since $\hat{\mu}_{1},\hat{\mu}_{2}$ from Definition 4.1
both have finite support, we can compute
$\mathcal{W}_{p}(\hat{\mu}_{1},\hat{\mu}_{2})$ through a linear optimization
problem by the optimal transport. Once we have checked that the ambiguity set
is not empty, we next provide a tractable reformulation of (3). We first point
out that the strong duality still holds, similar to the case where the
ambiguity set is a single Wasserstein ball.

[Regularity of Cost Function] The expected cost is bounded such that
$|\mathbb{E}_{Y\sim\hat{\mu}_{1}}[c(z,Y)]|<\infty,|\mathbb{E}_{Y\sim\hat{\mu}_{2}}[c(z,Y)]|<\infty,\forall
z\in\mathcal{Z}$.

######  Theorem 3.2 (Strong Duality of 2W-DRO)

Suppose Assumptions 2.1 and 3.1 holds. Then for
$\mathcal{A}=\\{\mu:\mathcal{W}_{p}(\mu,\hat{\mu}_{1})\leq\varepsilon_{1},\mathcal{W}_{p}(\mu,\hat{\mu}_{2})\leq\varepsilon_{2}\\}$,
we have $c_{P}(z)=c_{D}(z)$ with:

$\displaystyle c_{P}(z)$
$\displaystyle\coloneqq\sup_{\mu\in\mathcal{A}}\mathbb{E}_{Y\sim\mu}[c(z,Y)],$
(5) $\displaystyle c_{D}(z)$
$\displaystyle\coloneqq\inf_{\lambda_{1},\lambda_{2}\geq
0}\left\\{\lambda_{1}\varepsilon_{1}^{p}+\lambda_{2}\varepsilon_{2}^{p}+\sup_{\mathbb{P}\in\mathcal{P}(\hat{\mu}_{1},\hat{\mu}_{2})}\mathbb{E}_{\mathbb{P}}\left[\sup_{y\in\mathcal{Y}}\left\\{c(z,y)-\lambda_{1}\|\tilde{Y}_{1}-y\|_{1}^{p}-\lambda_{2}\|\tilde{Y}_{2}-y\|_{1}^{p}\right\\}\right]\right\\},$

where in $c_{D}(z)$, the first suprema set
$\mathcal{P}(\mu_{1},\mu_{2})=\\{\mathbb{P}_{(\tilde{Y}_{1},\tilde{Y}_{2})}\in\mathcal{P}_{1}(\mathbb{R}^{2d_{y}})\mid\mathbb{P}_{\tilde{Y}_{1}}=\hat{\mu}_{1},\mathbb{P}_{\tilde{Y}_{2}}=\hat{\mu}_{2}\\}$,
i.e., the set of all potential joint distribution supported on
$\mathbb{R}^{2d_{y}}$ where the marginal distribution of the first $d_{y}$
components is $\hat{\mu}_{1}$, and the marginal distribution of the last
$d_{y}$ components is $\hat{\mu}_{2}$.

When the ambiguity set is a single Wasserstein ball, we recover the standard
Wasserstein duality result. For this generalized duality result in (5), we
first verify that $\tilde{c}_{z}((\tilde{Y}_{1},\tilde{Y}_{2}),y)\coloneqq
c(z,Y)-\lambda_{1}\|Y_{1}-y\|_{1}^{p}-\lambda_{2}\|\tilde{Y}_{2}-y\|_{1}^{p}$
satisfies an augmented version of interexchangebility principle in Zhang et
al. (2022), which gives us the suprema of probability distributions over the
two marginal probability couplings. Then we demonstrate that this maximization
over the two marginal probability couplings is equivalent to the maximization
over the joint distribution (i.e., $\mathcal{P}(\hat{\mu}_{1},\hat{\mu}_{2})$)
in (5).

######  Remark 3.3 (Generalization of Strong Duality)

We can generalize the duality result in (5) with the intersection of more than
two Wasserstein balls under a generalized definition of Wasserstein distance
compared with Definition 2.1. We leave the discussion in Appendix 9.2.

Note that the current formulation (5) is still not tractable due to supremas
of both $\mathbb{P}\in\mathcal{P}(\hat{\mu}_{1},\hat{\mu}_{2})$ and
$y\in\mathcal{Y}$ inside the minimization problem. Therefore, we need to
remove this suprema to transform it into a tractable optimization problem. We
assume: [Structure of Cost Function] $c(z,y)=\max_{k\in[K]}c_{k}(z,y)$, where
each component $c_{k}(z,y)$ (abbreviated as $c_{z,k}(y)$) is proper, convex
and upper semi-continuous in $y$. Formally, when we use 1-Wasserstein distance
in $\mathcal{A}$, our general theorem of a tractable reformulation is stated
as follows:

######  Theorem 3.4 (Convex Reformulation of 2W-DRO)

Suppose Assumptions 2.1, 3.1 and 3.3 hold and $p=1$ in (3), we have:

$\displaystyle c_{P}(z)=\inf_{\lambda_{1},\lambda_{2}\geq
0,\bm{u},\bm{v},\bm{s},\bm{\alpha},\bm{\beta}}{}$
$\displaystyle\lambda_{1}\varepsilon_{1}+\lambda_{2}\varepsilon_{2}-\sum_{i=1}^{n}w_{1,i}u_{i}-\sum_{j=1}^{n}w_{2,j}v_{j}$
(6) s.t.
$\displaystyle-(u_{i}+v_{j})\geq[-c_{z,k}+\mathbf{1}_{\mathcal{Y}}]^{*}(\alpha_{i,k}+\beta_{j,k})-\alpha_{i,k}^{\top}y_{1,i}-\beta_{j,k}^{\top}y_{2,j},\forall
i,j\in[n],k\in[K];$
$\displaystyle\|\alpha_{i,k}\|_{\infty}\leq\lambda_{1},\|\beta_{j,k}\|_{\infty}\leq\lambda_{2},\forall
i,j\in[n],k\in[K].$

Compared with the standard Wasserstein DRO duality, our proof requires another
layer of duality to resolve the second suprema, i.e.
$\mathcal{P}(\hat{\mu}_{1},\hat{\mu}_{2})$ in (5), which is new to literature.
This duality is difficult for general $\hat{\mu}_{1},\hat{\mu}_{2}$ while the
transformations are possible in our case, when the two distributions
$\hat{\mu}_{1},\hat{\mu}_{2}$ have finite discrete support. We observe that
the joint distribution in $\mathcal{P}(\hat{\mu}_{1},\hat{\mu}_{2})$ can be
encoded as a transportation plan (joint probability for each pair)
$\\{\gamma_{i,j}\\}_{i,j\in[n]}$ where each $\gamma_{i,j}$ denotes the
transportation amount (size of probability)
$\Pr_{\mathbb{P}}[(\tilde{Y}_{1},\tilde{Y}_{2})=(y_{1,i},y_{2,j})]=\gamma_{i,j}$
and the cost of each pair denotes the perturbed loss
$s_{i,j}\coloneqq\sup_{y\in\mathcal{Y}}\left[c(z,y)-\lambda_{1}\|y_{1,i}-y\|_{1}-\lambda_{2}\|y_{2,j}-y\|_{1}\right]$.
That is:

$\sup_{\mathbb{P}\in\mathcal{P}(\hat{\mu}_{1},\hat{\mu}_{2})}\mathbb{E}_{\mathbb{P}}\left[\sup_{y\in\mathcal{Y}}\left\\{c(z,y)-\lambda_{1}\|\tilde{Y}_{1}-y\|_{1}^{p}-\lambda_{2}\|\tilde{Y}_{2}-y\|_{1}^{p}\right\\}\right]=-\inf_{\bm{\gamma}}\sum_{i,j\in[n]}\gamma_{i,j}(-s_{i,j}).$

Then we can use the duality of the classical transportation problem with dual
variables $u_{i}$ and $v_{j}$ and obtain (6).

###### Remark 3.5

In terms of the additional theoretical analysis under multiple constraints,
our analysis is in sharp contrast to the naive case of $\mu$ being a _finite-
dimensional variable_ , i.e., $\\{\mu\in\mathbb{R}^{d_{x}}:A\mu\leq b,C\mu\leq
d\\}\Leftrightarrow\\{\mu\in\mathbb{R}^{d_{x}}:[A;C]^{\top}\mu\leq[b;d]^{\top}\\}$)
since the decision variable $\mu$ is a _distribution of infinite-dimensional_.
Then as we mentioned above, obtaining the tractable reformulation of the
intersection of two distribution constraints is more difficult than the one
Wasserstein case, which involves another layer of suprema.

More specifically, when the cost function is piecewise linear with uncertainty
region $\mathcal{Y}$ being a polyhedron, we can obtain the following linear
optimization reformulation:

######  Corollary 3.6 (Linear Reformulation of 2W-DRO)

Suppose Assumptions 2.1, 3.1 and 3.3 hold, where each objective component
$c_{z,k}=a_{z,k}^{\top}y+b_{z,k}$ with $a_{z,k}$ and $b_{z,k}$ being linear in
$z$, and the uncertainty region is
$\mathcal{Y}=\\{y\in\mathbb{R}^{d_{y}}:Ey\leq f\\}$ for some $E,f$. When $p=1$
in (3), the optimization problem (6) can be further simplified to:

$\displaystyle c_{P}(z)=\inf_{\lambda_{1},\lambda_{2}\geq
0,\bm{u},\bm{v},\bm{\gamma},\bm{\alpha},\bm{\beta}}{}$
$\displaystyle\lambda_{1}\varepsilon_{1}+\lambda_{2}\varepsilon_{2}-\sum_{i=1}^{n}w_{1,i}u_{i}-\sum_{j=1}^{n}w_{2,j}v_{j}$
(7) s.t. $\displaystyle-(u_{i}+v_{j})\geq
b_{z,k}-(\alpha_{i,k}^{\top}y_{1,i}+\beta_{j,k}^{\top}y_{2,j})+\gamma_{i,j,k}^{\top}f,\forall
i,j\in[n],k\in[K];$
$\displaystyle\alpha_{i,k}+\beta_{j,k}=E^{\top}\gamma_{i,j,k}-a_{z,k},\forall
i,j\in[n],k\in[K];$
$\displaystyle\|\alpha_{i,k}\|_{\infty}\leq\lambda_{1},\|\beta_{j,k}\|_{\infty}\leq\lambda_{2},\gamma_{i,j,k}\geq
0,\forall i,j\in[n],k\in[K].$

As a consistency check, we see that when one ambiguity size
$\varepsilon_{2}\to\infty$, the intersection of the two Wasserstein balls
reduces to the standard result of one Wasserstein ball in Esfahani and Kuhn
(2018). Then we present problem reformulations of piecewise linear objectives
that will be used in our numerical studies as special cases in (7).

######  Example 3.7 (Mean-CVaR Portfolio Selection)

We consider the mean-CVaR portfolio optimization, where the measure CVaR
captures the trade-off between risk and the expected utility. Here
$c(z,y)=\text{CVaR}_{\rho}(y^{\top}z)-y^{\top}z$, where
$Y\in\mathbb{R}^{d_{y}}$ represents the random asset returns and
$z\in\mathcal{Z}=\\{z\in\mathbb{R}^{d_{z}}|\sum_{i\in[d_{z}]}z_{i}=1,z_{i}\geq
0\\}$ represents the allocation vector with $d_{y}=d_{z}$. Using the linear
reformulation of CVaR from Rockafellar and Uryasev (2002), the portfolio
framework can be reformulated with the augmented objective through an
augmented $(z,r)$ with $r$ being the $\eta$-quantile VaR that fits into the
framework (7):

$\displaystyle\inf_{z}\sup_{\mu_{y}\in\mathcal{A}}\mathbb{E}_{y}[c(z,y)]$
$\displaystyle=\min_{z}\max_{\mu_{y}\in\mathcal{A}}\mathbb{E}_{y}[\text{CVaR}_{\rho}(y^{\top}z)-y^{\top}z]$
(8)
$\displaystyle=\inf_{z,r}\left\\{r-y^{\top}z+\frac{1}{\rho}\max_{\mu_{y}\in\mathcal{A}}\mathbb{E}_{y}[\max\\{-y^{\top}z-r,0\\}]\right\\}$
$\displaystyle=\min_{z,r}\max_{\mu_{y}\in\mathcal{A}}\mathbb{E}_{y}\left[\max\left\\{-\left(1+\frac{1}{\rho}\right)y^{\top}z+\left(1-\frac{1}{\rho}\right)r,-y^{\top}z+r\right\\}\right].$

Here we do not impose linear constraints for the distribution of assets’
returns. That is $E=0\in\mathbb{R}^{d_{1}\times d_{y}}$, $f=0$,
$\mathcal{Y}=\mathbb{R}^{d_{y}}$. Then the optimization problem is reduced to:

$\displaystyle c_{P}(z)=\inf_{\lambda_{1},\lambda_{2}\geq
0,r,\bm{u},\bm{v},\bm{\alpha},\bm{\beta}}{}$
$\displaystyle\lambda_{1}\varepsilon_{1}+\lambda_{2}\varepsilon_{2}-\sum_{i=1}^{n}w_{1,i}u_{i}-\sum_{j=1}^{n}w_{2,j}v_{j}$
(9) s.t.
$\displaystyle-(u_{i}+v_{j})\geq(1-\frac{1}{\rho})r-(\alpha_{i,1}^{\top}y_{1,i}+\beta_{j,1}^{\top}y_{2,j}),\forall
i,j\in[n];$ $\displaystyle-(u_{i}+v_{j})\geq
r-(\alpha_{i,2}^{\top}y_{1,i}+\beta_{j,2}^{\top}y_{2,j}),\forall i,j\in[n];$
$\displaystyle\alpha_{i,1}+\beta_{j,1}=(1+\frac{1}{\rho})z,\forall i,j\in[n];$
$\displaystyle\alpha_{i,2}+\beta_{j,2}=z,\forall i,j\in[n];$
$\displaystyle\|\alpha_{i,1}\|_{\infty},\|\alpha_{i,2}\|_{\infty}\leq\lambda_{1},\|\beta_{j,1}\|_{\infty},\|\beta_{j,2}\|_{\infty}\leq\lambda_{2},\forall
i,j\in[n].$

######  Example 3.8 (Newsvendor Problem)

Consider the newsvendor loss $c(z,y)=\max\\{b\cdot y-bz,-h\cdot y+hz\\}$, and
$E=[-1;1]^{\top},f=(-L,U)^{\top}$ in (7), where the feasible region of demand
$Y$ is $\mathcal{Y}=\\{Y\in\mathbb{R}:L\leq Y\leq U\\}.$ Therefore, the primal
problem can be reformulated as:

$\displaystyle c_{P}(z)=\min_{\lambda,s_{i},\bm{\gamma},\bm{u},\bm{v}}$
$\displaystyle~{}\lambda_{1}\varepsilon_{1}+\lambda_{2}\varepsilon_{2}-\sum_{i=1}^{n}w_{1,i}u_{i}-\sum_{j=1}^{n}w_{2,j}v_{j}$
(10) s.t. $\displaystyle~{}-(u_{i}+v_{j})\geq-
bz-(\alpha_{i,1}^{\top}y_{1,i}+\beta_{j,1}^{\top}y_{2,j})-\gamma_{i,j,1}L+\gamma_{i,j,2}U,\forall
i,j\in[n];$ $\displaystyle~{}-(u_{i}+v_{j})\geq
hz-(\alpha_{i,2}^{\top}y_{1,i}+\beta_{j,2}^{\top}y_{2,j})-\gamma_{i,j,3}L+\gamma_{i,j,4}U,\forall
i,j\in[n];$
$\displaystyle~{}\alpha_{i,1}+\beta_{j,1}=-b-\gamma_{i,j,1}+\gamma_{i,j,2},\forall
i,j\in[n];$
$\displaystyle~{}\alpha_{i,2}+\beta_{j,2}=h-\gamma_{i,j,3}+\gamma_{i,j,4},\forall
i,j\in[n];$
$\displaystyle~{}\|\alpha_{i,k}\|_{\infty}\leq\lambda_{1},\|\beta_{j,k}\|\leq\lambda_{2},\gamma_{i,j,k}\geq
0,\forall i,j\in[n],k\in[4].$

The intersection of two Wasserstein balls $\mathcal{A}$ constructed here is
new in the stochastic optimization literature. Although some previous works
consider the intersection of distribution constraints in the ambiguity set,
instead of ensembling the power of two reference measures in downstream
decision-making in our setup, they consider some different formulations. For
example, moment ambiguity (Chen et al. 2019), Wasserstein and Goodness-of-
Fitness test ambiguity (Tanoumand et al. 2023).

##  4 Statistical Analysis of Two Reference Measures

In this section, we discuss the statistical properties of 2W-DRO for (4).
First, we focus on the two specific $\hat{\mu}_{1},\hat{\mu}_{2}$ that are
used in the without and with covariate shifts and the reasons for using them
as the intersection of two Wasserstein balls. Before proceeding to the
statistical analysis, we introduce a number of technical conditions and
assumptions.

######  Definition 4.1 (Regularity of the Conditional Distribution)

The conditional distribution $Y|X$ satisfies the regularity condition if
$Y|X=x\sim f(x)+\eta$ for all $x$, where

  * •

there exist $L>0,\beta\in(0,1]$ so that $f$ satisfies the $(L,\beta)$-Holder
condition under the Euclidean norm, i.e. $|f(x_{1})-f(x_{2})|\leq
L\|x_{1}-x_{2}\|_{2}^{\beta}$ for all $x_{1},x_{2}\in\mathcal{X}$;

  * •

$\mathbb{E}[\eta]=0,\mathbb{E}[\eta\eta^{\top}]\leq\sigma^{2}\mathbf{I},\exp(\zeta\|\eta\|_{2}^{\alpha})\in\mathcal{P}_{1}(\mathbb{R})$
for some constants $\zeta>0$ and $\alpha>p$ independent of $x$, where $p$ is
the order of Wasserstein distance in Definition 2.1.

We assume $\mu_{(X,Y)}$ satisfies Definition 4.1 and the marginal distribution
$\mu_{X}(\cdot)\geq c>0$ for some $c>0$ and all
$x\in[0,1]^{d_{x}}\subseteq\mathcal{X}$. Assumption 4 specifies the data
generating process of the historical training data, which is standard in the
literature. The lower bound of the measure $\mu_{x}(\cdot)$ guarantees that
for any covariate $x$, there are enough samples around $x$ when the sample
size tends to infinity under $\mu_{(X,Y)}$. The light-tail assumption
$\exp(\zeta\|\eta\|^{\alpha})\in\mathcal{P}_{1}(\mathbb{R})$ may appear
technical, but it is common in the literature to serve as a necessary
condition to bound the Wasserstein distance (Bolley et al. 2007, Fournier and
Guillin 2015).

###  4.1 Analysis without Covariate Shift

In this section, we provide the theoretical analysis of the model formulation
in (4), without covariate shifts, i.e. $\nu_{(X,Y)}=\mu_{(X,Y)}$. Since we do
not assume the form of the parametric function class $f\in\mathcal{F}$ in
Definition 4.1, in order for consistency for $\nu_{Y|x}$, we focus on the
statistical property of $\hat{\mu}_{1}$ and construct the conditional
estimator using nonparametric estimators (Györfi et al. 2006).

Among the nonparametric estimators, the Nadaraya-Watson kernel estimator
(Nadaraya 1964, Watson 1964) has been used widely in practice including
stochastic optimization (Hannah et al. 2010), dynamic programming (Hanasusanto
and Kuhn 2013) and newsvendor problem (Ban and Rudin 2019). We consider the
following $\hat{\mu}_{1}$ as the reference measure:

$\hat{\mu}_{1}=\mu_{n,x}\coloneqq\sum_{i=1}^{n}\frac{K(\frac{x-x_{i}}{h_{n}})}{\sum_{j=1}^{n}K(\frac{x-x_{j}}{h_{n}})}\delta_{y_{i}},$
(11)

where $K(\cdot):\mathbb{R}^{d_{x}}\mapsto[0,+\infty)$ is the kernel function
and $h_{n}\in\mathbb{R}^{+}$ is the bandwidth that depends on the sample size
$n$. Typically $K(x)$ decreases in $\|x\|_{2}$, allocating more weight to
local samples. For example, common choices of $K(x)$ include
$\mathbb{I}_{\left\\{\|x\|_{2}\leq 1\right\\}}$, $(1-\|x\|_{2}^{2})_{+}$ and
$\exp(-\|x\|_{2}^{2})$. By using (11) as the center of the Wasserstein ball,
we capture the following pattern of the historical data: when a sample
$(x_{i},y_{i})$ has a covariate $x_{i}$ close to $x$, the kernel function
$K((x-x_{i})/h_{n})$ tends to be large and the distribution $\hat{\mu}_{1}$
assigns more weight to $y_{i}$. And we weigh the historical samples more if
the covariate is closer to the current covariate $x$.

The kernel estimator is a nonparametric way to capture the joint distribution
of $(X,Y)$ in the historical data and leverage the information $X=x$. The
statistical properties of the kernel estimates in estimating
$\mathbb{E}[Y|X=x]$ have been documented in textbooks such as Györfi et al.
(2006), Tsybakov (2008). In contrast to existing theoretical analysis
primarily analyzing the properties of estimating $\mathbb{E}[Y|x]$ using
$\hat{\mathbb{E}}[Y|X=x]\coloneqq\sum_{i=1}^{n}\frac{K((x-x_{i})/h_{n})}{\sum_{j=1}^{n}K((x-x_{j})/h_{n})}y_{i}$,
our new technical contribution here is to provide general finite-sample
measure concentration results of estimating the entire distribution
$\nu_{Y|x}$ using (11).

The next assumption is concerned with the kernel function for the Nadaraya-
Watson estimator. For the kernel function $K(\cdot)$ in (11), there exist
$R\geq r>0$ and $b_{R}\geq b_{r}>0$ with some $C,\gamma>0$ such that
$K(\cdot)$ satisfies $b_{r}\mathbb{I}_{\left\\{\|x\|_{2}\leq r\right\\}}\leq
K(x)\leq b_{R}$, $\mathbb{E}_{z\sim\mu_{x}}[K(n^{\alpha}z)]\leq
Cn^{-\alpha\gamma}$ for all $\alpha>0$, and $\int_{0}^{\infty}t^{\beta
p}\tilde{K}(t)dt\leq C_{\beta,p}<\infty$, where
$\tilde{K}(\cdot):[0,\infty]\mapsto\mathbb{R}^{+}$ is an induced kernel
function defined as $\tilde{K}(x)=K(\|x\|_{2})$. In other words, the kernel
function is lower bounded by a naïve kernel. Besides, the decreasing rate
cannot be too slow. We first establish the measure concentration for
$\mu_{Y|x}$ under the $p$-Wasserstein distance.

######  Theorem 4.2 (Measure Concentration of the Kernel Estimator)

Suppose $\mu_{(X,Y)}$ satisfies Assumption 4 and the kernel function in (11)
satisfies Assumption 4.1. If $d_{y}=1$, for $x\in[0,1]^{d_{x}}$, $t\in[0,1]$
and a sequence of $\\{h_{n}\\}_{n=1}^{+\infty}$ such that $h_{n}^{\beta}\to 0$
and $nh_{n}^{d_{x}}\to\infty$ as $n\to\infty$, we have:

$\mathbb{P}(\mathcal{W}_{p}(\mu_{Y|x},\mu_{n,x})>t+c_{0}Lh_{n}^{\beta})\leq
c_{1}\exp(-c_{2}nh_{n}^{d_{x}}t^{2}),$ (12)

where $c_{0},c_{1},c_{2}$ are independent of $n$, $t$ and $x$.

As far as we know, this is the first result that establishes the nonasymptotic
rate of the Nadaraya-Watson kernel estimator using the Wasserstein distance.
We also generalize our results to other nonparametric estimators including the
k-Nearest Neighbor (kNN) as an alternative nonparametric formulation in
Appendix LABEL:app:kNN.

To prove the theorem, we first investigate the case when $K(x)$ represents the
naïve kernel with $K(x)=\mathbb{I}_{\left\\{\|x\|_{2}\leq r\right\\}}$. That
is,

$\mu^{r}_{n,x}\coloneqq\sum_{i=1}^{n}\frac{\mathbb{I}_{\left\\{\|x-x_{i}\|_{2}\leq
rh_{n}\right\\}}}{\sum_{j=1}^{n}\mathbb{I}_{\left\\{\|x-x_{j}\|_{2}\leq
rh_{n}\right\\}}}\delta_{y_{i}}.$ (13)

Because of Assumption 4, we can consider the shifted measure

$\bar{\mu}^{r}_{n,x}:=\sum_{i=1}^{n}\frac{\mathbb{I}_{\left\\{\|x-x_{i}\|_{2}\leq
rh_{n}\right\\}}}{\sum_{j=1}^{n}\mathbb{I}_{\left\\{\|x-x_{j}\|_{2}\leq
rh_{n}\right\\}}}\delta_{\\{y_{i}-f(x_{i})+f(x)\\}}.$ (14)

We then control the distance $\mathcal{W}_{p}(\mu_{Y|x},\bar{\mu}^{r}_{n,x})$
and $\mathcal{W}_{p}(\bar{\mu}^{r}_{n,x},\mu_{n,x}^{r})$. For the first term
$\mathcal{W}_{p}(\mu_{Y|x},\bar{\mu}^{r}_{n,x})$, the benefit of analyzing
$\bar{\mu}_{n,x}^{r}$ is that after the shift in $y_{i}$,
$\\{y_{i}-f(x_{i})+f(x)\\}_{i=1}^{n}$ can be regarded as i.i.d. samples. This
allows us to leverage the results in Fournier and Guillin (2015) to control
that first term by focusing on the samples in the ball of radius $rh_{n}$.
Bounding the second term $\mathcal{W}_{p}(\bar{\mu}^{r}_{n,x},\mu_{n,x}^{r})$
is straightforward, only requiring to use a probability coupling
$\hat{\Pi}=\sum_{i=1}^{n}\frac{\mathbb{I}_{\left\\{\|x-x_{i}\|_{2}\leq
rh_{n}\right\\}}}{\sum_{j=1}^{n}\mathbb{I}_{\left\\{\|x-x_{j}\|_{2}\leq
rh_{n}\right\\}}}\delta_{(y_{i}-f(x_{i})+f(x),y_{i})}$ and apply the Hölder
continuity of $f(x)$ from Definition 4.1.

Then we extend the result from the naïve kernel to other kernels like
$K(\cdot)$ being the boxed kernel such that
$b_{r}\mathbb{I}_{\left\\{\|x\|_{2}\leq R\right\\}}\leq K(x)\leq
b_{R}\mathbb{I}_{\left\\{\|x\|_{2}\leq r\right\\}}$ for some $R>r>0$ and more
general kernels satisfying Section 4.1. Correspondingly,
$\bar{\mu}_{n,x}=\sum_{i=1}^{n}\frac{K((x-x_{i})/h_{n})}{\sum_{j=1}^{n}K((x-x_{i})/h_{n})}\delta_{\\{y_{i}-f(x_{i})+f(x)\\}}$
and the problem still reduces to controlling
$\mathcal{W}_{p}(\bar{\mu}_{n,x},\mu_{n,x})$ and
$\mathcal{W}_{p}(\mu_{Y|x},\bar{\mu}^{r}_{n,x})$. We bound the term
$\mathcal{W}_{p}(\mu_{Y|x},\bar{\mu}_{n,x}^{r})$ by the same probability
coupling as before and the properties from Section 4.1. We then obtain the
term $h_{n}^{\beta}$ related to the bias term in the standard kernel
estimator. We bound the second term
$\mathcal{W}_{p}(\bar{\mu}^{r}_{n,x},\mu_{Y|x})$ through a bridge
$\bar{\mu}_{n,x}^{(1)}:=\sum_{i\in I}\frac{K((x-x_{i})/h_{n})}{\sum_{j\in
I}K((x-x_{j})/h_{n})}\delta_{f(x)+y_{i}-f(x_{i})}$ where
$I=\\{i\in[n]:\|x-x_{i}\|_{2}\leq rh_{n}\\}$. The complete proof can be found
in Appendix LABEL:app:proof-ws-concentrate.

######  Remark 4.3 (Difference from Kannan et al. (2023))

The formulation of our $\mu_{n,x}$ in (11) looks very similar to that in
Kannan et al. (2023) if they set
$\hat{\mu}_{n,x}=\frac{1}{n}\sum_{i=1}^{n}\delta_{\\{y_{i}-\hat{f}(x_{i})+\hat{f}(x)\\}}$
with
$\hat{f}(x)=\sum_{i\in[n]}\frac{K((x-x_{i})/h_{n})}{\sum_{j\in[n]}K((x-x_{j})/h_{n})}y_{i}$.
However, the theoretical guarantee of our proposed $\mu_{n,x}$ can be derived
under milder assumptions than (4.1). For example, if we use the naïve kernel
$\mu_{n,x}$ in (13), the convergence of our proposed $\mu_{n,x}$ still holds
even when the noise $\eta_{x}$ depend on $x$ as long as
$\mathcal{W}(\eta_{x_{1}},\eta_{x_{2}})\to 0$ when $\|x_{1}-x_{2}\|$ is small.
This is because we do not take account of historical samples with their
$x_{i}$ that are not near the new covariate $x$. In contrast, such
heterogeneity has more influence on the $\hat{\mu}_{n,x}$ of Kannan et al.
(2023) since it estimates the noise for all samples by
$\\{y_{i}-\hat{f}(x_{i})\\}_{i=1}^{n}$ instead of samples near the new
covariate $x$, which requires that the distribution of $\eta$ is the same
across all samples.

Theorem 4.2 shows the result for $d_{y}=1$. For $d_{y}>1$, through a direct
application of the Wasserstein distance measure concentration result (i.e.,
Theorem 2 in Fournier and Guillin (2015)), we have:

###### Corollary 4.4

When $d_{y}\geq 2$ and other conditions are the same as in Theorem 4.2, we
have:

$\mathbb{P}(\mathcal{W}_{p}(\mu_{Y|x},\mu_{n,x})>t+c_{0}Lh_{n}^{\beta})\leq\begin{cases}c_{1}\exp(-c_{2}nh_{n}^{d_{x}}t^{2}),&\text{
if }p>d_{y}/2\\\ c_{1}\exp(-c_{2}nh_{n}^{d_{x}}(t/\log(2+1/t))^{2}),&\text{ if
}p=d_{y}/2\\\ c_{1}\exp(-c_{2}nh_{n}^{d_{x}}t^{d_{y}/p}),&\text{ if
}p\in[1,d_{y}/2)\\\ \end{cases}.$ (15)

The numerical results of $d_{y}>1$ can be found in Section 7.3.2.

For simplicity, between the following analysis and the end of Section 3, we
assume $p=1$ to facilitate our analysis. From Corollary 4.4, the concentration
result applies to $p\geq 2$ through a simple application of Theorem 2 in
Fournier and Guillin (2015) and Corollary 4.4.

We are interested in the measure concentration result because it informs us of
the size of the preferred ambiguity set directly and connects to the finite-
sample performance guarantee of the objective. One proper choice of
$\varepsilon$ in the ambiguity set $\mathcal{A}$ is to cover the true
conditional distribution $\nu_{Y|x}$ with a high probability (Esfahani and
Kuhn 2018). Setting $\varepsilon_{1}\geq
c_{0}Lh_{n}^{\beta}+c_{1}\left(\frac{\log(1/\alpha)}{nh_{n}^{d_{x}}}\right)^{\frac{1}{\max\\{d_{y},2\\}}}$
guarantees that $\mathcal{W}_{p}(\mu,\mu_{n,x})\leq\varepsilon_{1}$ with
probability at least $1-\alpha$ from Theorem 4.2. If we set
$\varepsilon_{2}\to\infty$ such that the ambiguity set $\mathcal{A}$ reduces
to $\\{\mu:\mathcal{W}_{p}(\mu,\mu_{n,x})\leq\varepsilon_{1}\\}$, then
$\mathcal{A}$ contains the true measure $\nu_{Y|x}$ with probability at least
$1-\alpha$. Let $\hat{z}_{n}$ be the optimal solution to (3), then this
measure concentration result implies that the objective $\hat{J}_{n}$ in (3)
provides an upper bound for the true performance
$\mathbb{E}_{\nu_{Y|x}}[c(\hat{z}_{n};Y)]$ as follows:

######  Proposition 4.5 (Performance Guarantee)

Suppose Assumptions 4 and 4.1 hold. Consider the ambiguity set $\mathcal{A}$
where $\hat{\mu}_{1}=\mu_{n,x}$ with
$h_{n}=\Theta(n^{-\frac{1}{\beta\max\\{d_{y},2\\}+d_{x}}})$,
$\varepsilon_{1}(\alpha)=\Omega\left(L\log(1/\alpha)n^{-\frac{\beta}{\beta\max\\{d_{y},2\\}+d_{x}}}\right)$,
and $\hat{\mu}_{2}$ is arbitrary with $\varepsilon_{2}=\infty$ in (3), then we
have:

$\mathbb{P}\left(\mathbb{E}_{\nu_{Y|x}}[c(\hat{z}_{n},Y)]\leq\hat{J}_{n}\right)\geq
1-\alpha.$ (16)

We set the bandwidth $h_{n}$ above to minimize the size of
$\varepsilon_{1}(\alpha)$ with respect to $n$. For such particular choice of
$h_{n}$, the corresponding $\varepsilon_{1}(\alpha)\to 0$ as $n$ goes to
infinity. This means that without covariate shifts, building ambiguity sets
centered at nonparametric estimators allows the solution to converge to the
optimal ones as $\varepsilon_{1}\to 0$. Given the same probability $1-\alpha$,
if the true distribution $\nu_{Y|x}$ is contained in the Wasserstein ball
$\\{\mu:\mathcal{W}_{p}(\mu,\mu_{n,x})\leq\varepsilon_{1}\\}$ with a smaller
$\varepsilon_{1}$, then we have a tighter performance upper bound from
Proposition 4.5 since the value of the upper bound $\hat{J}_{n}$ increases
with $\varepsilon_{1}$. This chain of argument provides guarantees why we are
investigating a Wasserstein ball with smaller $\varepsilon$. And we would
defer more discussion on such theoretical connections later in Section 5.2.

In summary, when there is no distribution shift, nonparametric estimators such
as the Nadaraya-Watson kernel estimation have been shown to achieve good
performance asymptotically. Next, we focus on the case with covariate shift,
in which case parametric models may provide additional benefits.

###  4.2 Analysis with Covariate Shift

When there is a covariate shift as in Definition 2.2, we first show that
nonparametric approaches such as the Nadaraya-Watson estimator can lead to
erroneous outcomes. The main intuition is that the shift may lead to a lack of
close “neighbors” belonging to the training set on which the kernel weights of
the test points are derived. For the exposition, we set $p=\beta=1$ throughout
this subsection but similar results can be derived for other $p,\beta$.

######  Proposition 4.6 (Fragility of Nonparametric Approaches under Covariate
Shift)

Suppose Assumptions 4 and 4.1 hold. Consider $\nu_{(X,Y)}$ in Definition 2.2
and the new covariate $x\in\mathcal{X}$ with $\nu_{X}(x)\in(\nu_{L},\nu_{U})$
for some $\nu_{L},\nu_{U}>0$. For the Nadaraya-Watson kernel estimator
$\mu_{n,x}$, with $h_{n}=\Theta(n^{-\frac{1}{d_{x}+\max\\{d_{y},2\\}}})$, we
have:

$\mathcal{W}_{1}(\nu_{Y|x},\mu_{n,x})=O_{p}\left(\left(\frac{c_{x}}{n}\right)^{\frac{1}{d_{x}+\max\\{d_{y},2\\}}}\right),$
(17)

where $c_{x}:=\frac{\nu_{X}(x)}{\mu_{X}(x)}\in[0,\infty]$.

The key quantity in the result is $c_{x}$, the Radon-Nykodym derivative of the
test measure of the covariate with respect to the measure in the training
data. This term $c_{x}$ emerges from the original right-hand side of (17),
which is $O_{p}((n\mu_{X}(x))^{-1/(d_{x}+\max\\{d_{y},2\\}})$, while
$\nu_{X}(x)$ remains constant. When the two covariate distributions diverge
significantly (i.e., a substantial covariate shift), especially when
$\nu_{X}(x)\gg\mu_{x}(x)$, the concentration of the Wasserstein distance is a
lot weaker because of the large $c_{x}$. This is because of the “local
smoothing” feature of many nonparametric estimators. If there are no neighbors
to learn from, such estimators may perform poorly. In addition, we can show
that $\left(\frac{c_{x}}{n}\right)^{\frac{1}{d_{x}+\max\\{d_{y},2\\}}}$ is the
lower bound of the convergence rate of the Wasserstein distance, which makes
the dependence on $c_{x}$ inevitable.

######  Remark 4.7 (Tightness of the Dependence on $c_{x}$)

We consider a setup with $d_{x}=1,d_{y}=2$ and use the naïve kernel
$K(x)=\mathbb{I}_{\left\\{\|x\|_{2}\leq r\right\\}}$ in the estimation of
$\nu_{n,x}$. Then for any covariate $x\sim\nu_{X}$, denote
$\tilde{\nu}_{Y|x}^{n}$ as the local weighted average distribution around
$\nu_{Y|x}$, i.e.
$\tilde{\nu}_{Y|x}^{n}:=\frac{\mathbb{E}_{\tilde{x}\sim\nu_{X}}[\nu_{Y|\tilde{x}}\mathbb{I}_{\left\\{\|\tilde{x}-x\|_{2}\leq
rh_{n}\right\\}}]}{\mathbb{P}(\tilde{x}:\|\tilde{x}-x\|_{2}\leq rh_{n})}$.
Conditioned on the event
$N_{n}(:=\sum_{i=1}^{n}\mathbb{I}_{\left\\{\|x-x_{i}\|\leq{rh_{n}}\right\\}})\geq\frac{rnh_{n}}{2c_{x}}$,
which occurs with a high probability, then by the triangle inequality, we
have:

$\mathbb{E}[\mathcal{W}_{1}(\nu_{Y|x},\mu_{n,x})]\geq\max\\{\mathbb{E}[\mathcal{W}_{1}(\tilde{\nu}_{Y|x}^{n},\mu_{n,x})]-\mathcal{W}_{1}(\nu_{Y|x},\tilde{\nu}_{Y|x}^{n}),\mathcal{W}_{1}(\nu_{Y|x},\tilde{\nu}_{Y|x}^{n})-\mathbb{E}[\mathcal{W}_{1}(\tilde{\nu}_{Y|x}^{n},\mu_{n,x})]\\},$

where we take the expectation over the randomness of $Y$ above. Note that from
Dudley (1969), Weed and Bach (2019),
$\mathbb{E}[\mathcal{W}_{1}(\tilde{\nu}_{Y|x}^{n},\nu_{n,x})]=\Omega\left(\sqrt{\frac{c_{x}}{nh_{n}}}\right)$
since $\nu_{n,x}$ can be regarded as the empirical distribution with $N_{n}$
points corresponding to $\tilde{\nu}_{Y|x}^{n}$. And
$\mathcal{W}_{1}(\nu_{Y|x},\tilde{\nu}_{Y|x}^{n})\asymp rh_{n}$ when $f(x)$ is
locally linear in terms of $x$ for $\|\tilde{x}-x\|_{2}\approx rh_{n}$.
Therefore, we have:
$\mathbb{E}[\mathcal{W}_{1}(\nu_{Y|x},\mu_{n,x})]=\Omega\left(\sqrt{\frac{c_{x}}{nh_{n}}}\right)=\Omega((\frac{c_{x}}{n})^{\frac{1}{3}})$
by choosing $h_{n}=\Theta\left((\frac{c_{x}}{n})^{\frac{1}{3}}\right)$, which
matches the upper bound result in Proposition 4.6.

Such a dependence makes the nonparametric estimator and the contextual
stochastic optimization framework unreliable when a new covariate $x$ happens
to be drawn from less covered regions by the training set, particularly in the
case when $c_{x}$ is significantly larger than $nh_{n}^{d_{x}}$. For example,
it can happen when the training data is collected based on customers of
certain demographics and the decision needs to be applied on a different
group. In this case, parametric estimators may offer benefits. In particular,
although parametric estimators (such as linear models) can suffer from model
misspecification, they provide global estimation that doesn’t rely on the
presence of “neighbors”. Using parametric estimators as the base for
contextual optimization has been proposed in Kannan et al. (2023) with the
following assumption: [Parametric Regression Estimator] There exists a
parametric regression that outputs a conditional distribution
$\hat{\mu}_{n,x}$ such that, $\forall t\in[0,1]$, we have:

$\mathbb{P}(\mathcal{W}_{p}(\mu_{Y|x},\hat{\mu}_{n,x})\geq\mathcal{E}_{apx}+t)\leq
c_{1}\exp(-c_{2}nt^{\max\\{d_{y},2\\}}),\forall x\in\mathcal{X}.$ (18)

Examples of such conditional distributions and the theoretical guarantee are
presented in Kannan et al. (2023), Iyengar et al. (2023), where
$\mathcal{E}_{apx}>0$ is a uniform constant over the covariate domain
$\mathcal{X}$ incorporating the model misspecification. We simply impose it as
an assumption and will show one example of parametric estimators that satisfy
the condition. This implies that parametric regression estimators may be less
affected by covariate shift compared with nonparametric estimators since
$\mathcal{E}_{apx}$ is already a uniform term over $\mathcal{X}$.

###### Example 4.8

Consider $\mu_{(X,Y)}$ in Section 4, where
$f(x)=B^{*}x+g(x)\mathbf{1}_{d_{y}}$ in Definition 4.1 with
$B^{*}\in\argmin_{\beta}\mathbb{E}_{\mu_{(X,Y)}}[\|y-Bx\|^{2}]$. Denote
$\hat{B}$ as the OLS estimator
$\hat{B}\in\argmin_{B}\frac{1}{n}\sum_{i=1}^{n}\|y_{i}-Bx_{i}\|^{2}$, then the
distribution estimator
$\hat{\mu}_{n,x}=\frac{1}{n}\sum_{i=1}^{n}\delta_{\\{\hat{B}x+y_{i}-\hat{B}x_{i}\\}}$
satisfies Assumption 4.7 under some additional technical conditions listed in
Appendix 10.

In this case, $\mathcal{E}_{apx}$ exists due to the ignorance of the
misspecification term $g(x)$ here. In reality, for the stochastic error term
in Proposition 4.8, we are estimating the “noise” $g(x)+\eta$ instead of
$\eta$. Note that the above result (e.g. Assumption 4.7) holds for any $x$,
and the error rate does not depend on the nearby samples around the covariate
$x\in\mathcal{X}$. This makes the bound independent of the effects from any
covariate shift. Due to the uniform bound in $\mathcal{E}_{apx}$ for any $x$,
we obtain the following results when we have covariate shifts, meaning that
the upper bound of the estimation error does not get larger compared with the
problem without covariate shifts.

######  Corollary 4.9 (Robustness of Parametric Approaches under Covariate
Shift)

For the conditional distribution estimator satisfying Assumption 4.7, we have:

$\mathbb{P}(\mathcal{W}_{p}(\nu_{Y|x},\hat{\mu}_{n,x})\geq\mathcal{E}_{apx}+t)\leq
c_{1}\exp(-c_{2}nt^{\max\\{d_{y},2\\}}),\forall x\in\mathcal{X}.$ (19)

This result directly follows from the fact that we take a uniform argument of
$\mathcal{E}_{apx}$ over the same covariate domain and $\nu_{Y|x}=\mu_{Y|x}$.

Finally, based on these theoretical measure concentration results for the two
nonparametric and parametric distribution estimators, we discuss the
statistical measure concentration result of 2W-DRO. That is, in (4), we set
the expressions of these two distribution estimators as:

$\displaystyle\hat{\mu}_{1}$
$\displaystyle=\mu_{n,x}^{NP}\coloneqq\sum_{i=1}^{n}\frac{K((x-x_{i})/h_{n})}{\sum_{j=1}^{n}K((x-x_{j})/h_{n})}\delta_{y_{i}},$
(20) $\displaystyle\hat{\mu}_{2}$
$\displaystyle=\mu_{n,x}^{P}\coloneqq\frac{1}{n}\sum_{i=1}^{n}\delta_{\\{\hat{f}(x)+y_{i}-\hat{f}(x_{i})\\}}.$

Recall our main model in Section 2, we detail the ambiguity formulation of
2W-DRO as follows: i.e., where the size of $\varepsilon_{1},\varepsilon_{2}$
is set to contain the true distribution $\nu_{Y|x}$ with high probability from
the discussion in Theorem 4.10. Intuitively, $\mathcal{A}$ contains the true
distribution $\nu_{Y|x}$ with high probability. We measure the statistical
performance of the two estimators through the size (or the convergence rate)
of the Wasserstein ball to cover $\nu_{Y|x}$:

######  Theorem 4.10 (Choice of Ambiguity Size in 2W-DRO)

When we replace $\hat{\mu}_{1},\hat{\mu}_{2}$ in (4) with $\mu_{n,x}^{NP}$ and
$\mu_{n,x}^{P}$ in (20) respectively, if we set the corresponding
$\varepsilon_{1},\varepsilon_{2}$ such that:

$\varepsilon_{1}\geq\varepsilon_{1,lw}\left(:=\Theta\left(\left(\frac{c_{x}}{n}\right)^{\frac{1}{d_{x}+\max\\{d_{y},2\\}}}\right)\right),\quad\varepsilon_{2}\geq\varepsilon_{2,lw}\left(:=\mathcal{E}_{apx}+\Theta\left(\left(\frac{1}{n}\right)^{\frac{1}{\max\\{d_{y},2\\}}}\right)\right).$
(21)

Then with a high probability (e.g. $1-\alpha$ for some constant $\alpha$), the
true conditional distribution $\nu_{Y|x}\in\mathcal{A}$ for $\mathcal{A}$ in
(3).

Note that in (21), $\varepsilon_{1,lw}$ converges to 0 only if $c_{x}=o(n)$.
When $c_{x}$ is of the constant level, we have the nonparametric convergence
rate $O(n^{-1/(d_{x}+\max\\{d_{y},2\\})})$; And the lower bound of
$\varepsilon_{2,lw}$ is $\mathcal{E}_{apx}>0$ and but the convergence with
respect to $n$ not depend on the covariate shift size $c_{x}$. From the
discussion above, under covariate shifts with large $c_{x}$, the nonparametric
estimator $\hat{\mu}_{1}$ requires more samples to converge to ensure a
limited estimation error. For the parametric regression-based estimator
$\hat{\mu}_{2}$ mentioned above, we may suffer from some bounded approximation
error but the error convergence rate can still be good under finite samples in
terms of $n$. In practice, we do not observe the true scale of
$\mathcal{E}_{apx}$ and constants but are partially aware of the magnitude of
the covariate shift. The intersection of two ambiguity sets $\mathcal{A}$ can
give a tighter approximation for the ambiguity set $\mathcal{A}$ compared with
the ambiguity set being only one Wasserstein ball. Therefore, if either
$\varepsilon_{1,lw}$ or $\varepsilon_{2,lw}$ is large, the size of
$\mathcal{A}$ still depends on the smaller one, thus being not so
conservative. We defer the formal generalization statistical properties of
2W-DRO to Section 5.2.

In the computational perspective, our constructed $\mu_{n,x}^{NP}$ and
$\mu_{n,x}^{P}$ here satisfy Assumptions 2.1 and 3.1 in Section 3. Therefore,
we can borrow the result in Theorem 3.4 to reformulate the optimization
problem of (3) to a convex optimization problem when $p=1$.

We conclude this part with a high-level discussion of other useful scenarios
in the ambiguity intersection. Note that the contextual optimization problem
under covariate shift is a natural setup that intersects different estimators
in the DRO problem. However, our computational formulation derived here should
provide insights for other scenarios beyond contextual optimization. This
encompasses scenarios when we make decisions under the transfer learning setup
(Zhuang et al. 2020), i.e. we want to learn some parameters in the target
distribution while we have only a few samples from the target distribution but
lots more samples from the source distribution. Therefore, the parameter
estimation is of high variance if we only use the target data while the
parameter estimation is biased when learning solely from the source data. Then
the approach of intersecting two Wasserstein balls while setting each
reference measure as the empirical distribution of the source or target
distribution respectively would balance the risk well.

##  5 Computational Surrogate and Performance Comparison

Our computational framework derived in Theorem 3.2 is general for encompassing
scenarios of Wasserstein ambiguity intersection. However, even if the
corresponding optimization problem under $\mathcal{A}$ is tractable from
Theorem 3.4, this reformulated optimization problem is still inefficient since
the number of decision variables in (6) is $\Theta(n^{2})$. In contrast, the
number of decision variables in the standard Wasserstein-DRO problem is at
most $\Theta(n)$ (Esfahani and Kuhn 2018). This motivates us to consider
approaches to reduce computational costs. In this section, we propose one
special computational reduction approach by using a surrogate objective, i.e.,
IW-DRO and compare its performance guarantees with 2W-DRO.

###  5.1 IW-DRO: Surrogate Objective with an Interpolate Estimator

We detail the construction of the surrogate objective as follows.
Specifically, we set $\mathcal{A}$ as one Wasserstein ball with the reference
measure as the interpolation of these two distribution estimators
$\mu_{n,x}^{NP},\mu_{n,x}^{P}$. That is, we set the ambiguity set in (3) as
$\mathcal{A}^{IW}\coloneqq\\{\mu:\mathcal{W}_{p}(\mu,\hat{\mu}_{n,x})\leq\varepsilon\\}$
in (3) with:

$\hat{\mu}_{n,x}\coloneqq\kappa_{x}\mu_{n,x}^{NP}+(1-\kappa_{x})\mu_{n,x}^{P},$
(22)

where

$\kappa_{x}\coloneqq\max\left\\{1-\tau\left(\sum_{i=1}^{n}K((x-x_{i})/h_{n})\right)^{-\frac{1}{d_{y}}},0\right\\},$
(23)

with some threshold $\tau\in(0,\infty)$. Note that $\kappa_{x}$ depends on the
magnitude of covariate shifts. We first demonstrate the measure concentration
result of $\hat{\mu}_{n,x}$ as follows:

######  Theorem 5.1 (Measure Concentration of IW-DRO)

Suppose Assumptions 4, 4.1 and 4.7 hold. Consider $\nu_{(X,Y)}$ in Definition
2.2 and the new covariate $x\in\mathcal{X}$ with
$\nu_{X}(x)\in(\nu_{L},\nu_{U})$ for some $\nu_{L},\nu_{U}>0$. For the
interpolate estimator (22), following the same choice of $h_{n}$ as in
Proposition 4.6, we have:

$\mathcal{W}_{1}(\nu_{Y|x},\hat{\mu}_{n,x})\leq\min\left\\{O_{p}\left(\left(\frac{c_{x}}{n}\right)^{\frac{1}{d_{x}+\max\\{d_{y},2\\}}}\right),O_{p}\left(\left(\frac{1}{n}\right)^{\frac{1}{\max\\{d_{y},2\\}}}\right)+\mathcal{E}_{apx}\right\\}.$
(24)

This theorem shows that the interpolate estimator $\hat{\mu}_{n,x}$ nearly
serves as a smaller upper bound among
$\mathcal{W}_{1}(\mu_{Y|x},\mu_{n,x}^{NP})$ and
$\mathcal{W}_{1}(\mu_{Y|x},\mu_{n,x}^{P})$ up to a constant compared with
Theorem 4.10. And if we choose

$\hat{\varepsilon}\coloneqq\kappa_{x}\varepsilon_{1}+(1-\kappa_{x})\varepsilon_{2},$
(25)

where $\varepsilon_{1},\varepsilon_{2}$ are chosen in Theorem 4.10, the
corresponding objective which maximizes all the possible distributions over
$\mathcal{A}^{IW}(\supseteq\mathcal{A})$ is an upper bound for the true risk
following the same arguments as Proposition 4.5.

Computationally speaking, IW-DRO is easier than 2W-DRO. Given
$\hat{\mu}_{n,x}$ in (22), we set:

$w_{i}=\begin{cases}\frac{\kappa_{x}K((x-x_{i})/h_{n})}{\sum_{j=1}^{n}K((x-x_{j})/h_{n})},&i\in[n]\\\
\frac{1-\kappa_{x}}{n},&i\in\\{n+1,\ldots,2n\\}\end{cases}.$ (26)

Then we can represent
$\hat{\mu}_{n,x}=\sum_{i=1}^{n}w_{i}\delta_{y_{i}}+\sum_{i=n+1}^{2n}w_{i}\delta_{\\{\hat{f}(x)+y_{i}-\hat{f}(x_{i})\\}}$
with weights defined above in (26). Therefore, IW-DRO becomes the standard
Wasserstein DRO problem where the ambiguity set there is one single
Wasserstein ball with a reference measure being a discrete distribution with
$2n$ support points. Especially, under a large covariate shift where
$\kappa_{x}=0$, we have $w_{i}=0,\forall i\in[n]$; When there are no covariate
shifts such that enough samples are near the current covariate $x$ with
$\kappa_{x}=1$, the problem reduces to the classical kernel-DRO problem with
$w_{i}=\frac{K((x-x_{i})/h_{n})}{\sum_{j=1}^{n}K((x-x_{j})/h_{n})},i\in[n]$.
Then if we solve this computational surrogate
$\hat{J}_{n}=\inf_{z\in\mathcal{Z}}\sup_{\mu:\mathcal{W}_{1}(\mu,\hat{\mu}_{n,x})\leq\hat{\varepsilon}}\mathbb{E}_{\mu}[c(z,Y)]$,
we would obtain an optimization problem with the number of decision variables
being $O(n)$ following the same step as in Esfahani and Kuhn (2018), which
reduces the computational complexity compared with 2W-DRO.

###  5.2 Generalization Guarantees between IW-DRO and 2W-DRO

We denote $\hat{z}_{1},\hat{z}_{2}$ as the solution of the 2W-DRO and IW-DRO
model with $\mathcal{A},\mathcal{A}_{IW}$ representing the ambiguity set in
the 2W-DRO and IW-DRO models respectively. We compare the solution through the
lens of generalization performance as follows:

######  Corollary 5.2 (Objective Approximation)

Suppose Assumptions 4, 4.1 and 4.7 hold. When
$\varepsilon_{1},\varepsilon_{2},\kappa_{x},\hat{\varepsilon}$ are chosen
according to (21), (23) and (25) for $\mathcal{A}$ and $\mathcal{A}_{IW}$, the
following inequality holds with probability at least $1-\alpha$:

$\mathbb{E}_{\nu_{Y|x}}[c(z,Y)]\leq\sup_{\nu\in\mathcal{A}}\mathbb{E}_{\nu}[c(z,Y)]\leq\sup_{\nu\in\mathcal{A}_{IW}}\mathbb{E}_{\nu}[c(z,Y)],\forall
z\in\mathcal{Z}.$

As known in literature (Esfahani and Kuhn 2018, Van Parys et al. 2021), one
advantage of DRO models is that the objective can provide an upper bound of
true objective performance $\mathbb{E}[c(z,Y)]$ through the construction of
the ambiguity set. However, we may want a tighter upper bound through
different ambiguity sets under the same probability level. From Corollary 5.2,
compared with IW-DRO, 2W-DRO estimator is more accurate in the sense of
providing a tighter upper bound for $\mathbb{E}_{\nu_{Y|x}}[c(z,Y)]$ for all
$z\in\mathcal{Z}$ in the objective due to a smaller ambiguity set in 2W-DRO
while keeping the true distribution inside.

Furthermore, we can build finite-sample guarantees for the excess risk of the
proposed 2W-DRO and IW-DRO. If we choose $\varepsilon_{1},\varepsilon_{2}$
properly, $\forall x\in\mathcal{X}$, we obtain the following excess risk
bounding the effects of solution difference to the true cost function under
the true measure $\nu_{Y|x}$:

######  Theorem 5.3 (Excess Risk)

Suppose Assumptions 3.3, 4, 4.1 and 4.7 hold. When
$\varepsilon_{1},\varepsilon_{2},\kappa_{x},\hat{\varepsilon}$ are chosen
according to (21), (23) and (25) for $\mathcal{A}$ and $\mathcal{A}_{IW}$,
then $\forall x\in\mathcal{X}$, with probability at least $1-\alpha$,

$\displaystyle\mathbb{E}_{\nu_{Y|x}}[c(\hat{z}_{1},Y)]-\mathbb{E}_{\nu_{Y|x}}[c(z^{*},Y)]$
$\displaystyle\leq 2L_{z^{*}}\min\\{\varepsilon_{1},\varepsilon_{2}\\},$
$\displaystyle\mathbb{E}_{\nu_{Y|x}}[c(\hat{z}_{2},Y)]-\mathbb{E}_{\nu_{Y|x}}[c(z^{*},Y)]$
$\displaystyle\leq
2L_{z^{*}}(\kappa_{x}\varepsilon_{1}+(1-\kappa_{x})\varepsilon_{2})\leq
CL_{z}^{*}\min\\{\varepsilon_{1},\varepsilon_{2}\\},$

where $z^{*}\in\argmin_{z\in\mathcal{Z}}\mathbb{E}_{\nu_{Y|x}}[c(z,Y)]$ and
$L_{z^{*}}$ is the local Lipschitz constant of $c(z,Y)$ under $z^{*}$. And $C$
is some constant with $C>2$.

The ambiguity size proposed in the measure concentration part is the key to
the upper bound of the excess risk even with finite samples, leading to
reasons why we minimize the lower bound of the ambiguity size with and without
covariate shift. This proof is based on the arguments by combining the
coverage guarantee and regularization effects of Wasserstein DRO models and
realizing that we attain a double distributionally robust effect from that. As
a special case, we have the generalization error with the ambiguity set
$\mathcal{A}_{i}=\\{\mu:\mathcal{W}(\mu,\hat{\mu}_{i})\leq\varepsilon_{i}\\}$
since it directly follows from controlling the size of
$\varepsilon_{1},\varepsilon_{2}$ in $\hat{z}_{1}$. We see that the (upper
bound) of the excess risk of the 2W-DRO solution is smaller than the one using
only one Wasserstein ball $\mathcal{A}_{i}$, indicating that the 2W-DRO
solution is robust and achieves “best of both worlds” regardless of covariate
shifts.

We also provide guarantees for the IW-DRO estimator $\hat{z}_{2}$. Compared
with the IW-DRO solution $\hat{z}_{2}$, the 2W-DRO solution $\hat{z}_{1}$
offers better excess risk guarantees up to some constants. Theoretically, IW-
DRO only worsens the 2W-DRO model performance by a constant level from the
excess risk upper bound, which can be ignored if we are more interested in the
rate concerning the sample size $n$. Therefore, IW-DRO may be preferred under
large-scale computational problems.

##  6 Discussions

###  6.1 Calibration through Constraining the Ambiguity Set: An Illustrative
Example

We highlight the utility of 2W-DRO through the following illustrative example
to show how $\mathcal{A}$ reduces conservativeness by making some unrealistic
distributions contained in one Wasserstein ball infeasible. Recall the setup
in Section 2, under a new covariate $X=x$ in the ambiguity set
$\mathcal{A}=\\{\mathcal{W}_{p}(\mu,\hat{\mu}_{1})\leq\varepsilon_{1},\mathcal{W}_{p}(\mu,\hat{\mu}_{2})\leq\varepsilon_{2}\\}$,
we set $\hat{\mu}_{2}$ as the empirical approximation of
$N(\hat{\theta}^{\top}x,\hat{\sigma}^{2})$, where $\hat{\theta},\hat{\sigma}$
are obtained through parametric estimation from $\mathcal{D}_{n}$ and
$\hat{\mu}_{2}$ satisfies Remark 4.7.

(a) Density Plots of Distributions without Covariate Shift

(b) Density Plots of Distributions with Covariate Shift

Figure 3: Illustration with Covariate Shift Problems. D2-P1 is the
distribution $P_{1}:N(\hat{\theta}^{\top}x-\varepsilon_{2},\hat{\sigma}^{2})$
and D2-P2 is the distribution
$P_{2}:N(\hat{\theta}^{\top}x+\varepsilon_{2},\hat{\sigma}^{2})$. D1-center
and D2-center denote $\hat{\mu}_{1},\hat{\mu}_{2}$, where $\hat{\mu}_{1}$ is a
weighted empirical distribution. D1-center is the kernel-weighted empirical
distribution $\hat{\mu}_{1}$ with 1000 samples belonging to
$\\{x:\|x\|_{2}\leq r\\}$ without covariate shifts while only with 4 samples
belonging to $\\{x:\|x\|_{2}\leq r\\}$ with covariate shifts.

We analyze the effects of imposing the other distribution constraint
$\\{\mu:\mathcal{W}_{p}(\mu,\hat{\mu}_{1})\leq\varepsilon_{1}\\}$ under
covariate shift. Here, we set $\hat{\mu}_{1}$ as the kernel estimator in (11)
with $K(x)=\mathbb{I}_{\left\\{\|x\|_{2}\leq r\right\\}}$ and the bandwidth
$h_{n}=O(n^{-1/(d_{x}+2)})$. We set
$\varepsilon_{1}=\Theta\left((\frac{c_{x}}{n})^{\frac{1}{d_{x}+2}}\right)$,
$\varepsilon_{2}=\Theta(n^{-1/2})+\mathcal{E}_{apx}$ such that the true
distribution $\nu_{Y|x}$ is contained in $\mathcal{A}$ with a high probability
from Theorem 4.10. We illustrate this through Figure 3, where D1-center and
D2-center distribution correspond to $\hat{\mu}_{1},\hat{\mu}_{2}$
respectively. And $D_{2}-P_{1},D_{2}-P_{2}$ represent two distributions
$P_{1},P_{2}$ belonging to
$\\{\mu:\mathcal{W}_{p}(\mu,\hat{\mu}_{2})\leq\varepsilon_{2}\\}$. More
specifically, we illustrate how distribution constraint
$\mathcal{W}_{p}(\mu,\hat{\mu}_{1})\leq\varepsilon_{1}$ makes $P_{1}$ that is
contained in $\\{\mu:\mathcal{W}_{p}(\mu,\hat{\mu}_{2})\leq\varepsilon_{2}\\}$
infeasible in $\mathcal{A}$.

If we do not have much covariate shift (Figure 3$(a)$), this additional
constraint $\mathcal{W}_{p}(\mu,\hat{\mu}_{1})\leq\varepsilon_{1}$ does not
make $P_{1}$ or $P_{2}$ infeasible in $\mathcal{A}$ since both
$\hat{\mu}_{1},\hat{\mu}_{2}$ are close to the true distribution and the
intersection of the two distribution balls then covers both $P_{1},P_{2}$.
However, since $\varepsilon_{1}$ is not large without large covariate shifts,
the downstream decision is not too pessimistic.

If we have much covariate shift (Figure 3$(b)$), this additional constraint
$\mathcal{W}_{p}(\mu,\hat{\mu}_{1})\leq\varepsilon_{1}$ makes $P_{1}$ feasible
in $\mathcal{A}$ since $P_{1}$ is far away from $\hat{\mu}_{1}$ (i.e.
$\mathcal{W}_{p}(P_{1},\hat{\mu}_{1})>\varepsilon_{1}$). In this case, the
first distribution constraint helps calibrate the ambiguity set from the
second constraint by shrinking the ambiguity set and make some unrealistic
distributions (i.e., $P_{1}$ here) infeasible in $\mathcal{A}$. And this makes
the downstream decision not affected by such $P_{1}$ and thus not being too
conservative.

Similarly, we can see incorporating the second Wasserstein ball with a
reference measure being $\hat{\mu}_{2}$ can make some unrealistic
distributions that satisfy
$\\{\mu:\mathcal{W}_{p}(\mu,\hat{\mu}_{1})\leq\varepsilon_{1}\\}$ infeasible
in $\mathcal{A}$. In all, the nature of intersecting Wasserstein balls in
2W-DRO reduces over-conservativeness of downstream decisions by removing
peculiar distributions through the ensemble learning of the nonparametric and
parametric estimators and ensures robust but not conservative decisions
regardless of covariate shifts.

###  6.2 Other Modeling Choice of Ambiguity under Covariate Shifts.

We optimize the downstream decision instead of just estimating
$\mathbb{E}[Y|x]$ as the statistical learning. If we directly interpolate
between the two DRO solutions obtained from a single Wasserstein ball, that
solution may be infeasible and less transparent. Some other alternatives to
reduce computational challenges can be developed in parallel. For example, we
can consider the following _switching ambiguity set_
$\mathcal{A}_{S}=\\{\mu:\mathcal{W}(\mu,\hat{\mu}_{i})\leq\varepsilon_{i}\\}$,
where $i=\begin{cases}1,~{}\text{if}~{}{\kappa_{x}\geq 1/2}\\\
2,~{}\text{otherwise}\end{cases}$. $\mathcal{A}_{S}$ switches the Wasserstein
ball depending on the magnitude of covariate shifts and can be seen as a
computational surrogate like IW-DRO. Compared to that, our 2W-DRO model is
still more general and flexible. Both theoretically and empirically, setting
the ambiguity set as one Wasserstein ball may be sometimes conservative in
practice especially when the two Wasserstein balls are both large under
covariate shifts. And our proposed intersection method can make some
“peculiar” distributions infeasible in $\mathcal{A}$ by combining the power of
both estimators and make our ambiguity set less conservative from Section 6.1.

###  6.3 Why we consider Wasserstein Distance in Robust Formulations.

In this paper, we propose a DRO model where the ambiguity set is the
intersection of two distribution balls with the probability “distance” being
the Wasserstein distance. Changing Wasserstein Distance to some standard
$\phi$-divergences (Ben-Tal et al. 2013) may be impossible since the
intersection of two $\phi$-divergence balls is always empty unless the two
reference measures have the same support, which is not the case in our
proposed $\mu_{n,x}^{NP}$ and $\mu_{n,x}^{P}$. However, some recently proposed
ambiguity sets characterizing the probability distance $D$ by some smoothed
$\phi$-divergence or sinkhorn distance may be useful to ensure nonemptyness.

##  7 Numerical Studies

We demonstrate the performance of our proposed models, i.e., 2W-DRO in Section
3 and the surrogate model IW-DRO in Section 5, in the two tasks of income
predictions and portfolio optimization, and on both synthetic and empirical
experiments. Since these two are both piecewise linear optimization problems,
we set $p=1$ in (4) and the problem reduces to linear optimization problems.
All optimization problems are implemented in Python via Gurobi on a PC.

###  7.1 Benchmarks

In this part, we only compare them with DRO models where $\mathcal{A}$ in (3)
is a single Wasserstein ball (i.e.,
$\mathcal{A}=\\{\mu:\mathcal{W}_{1}(\mu,\hat{\mu})\leq\varepsilon\\})$, which
helps assess the utility of in our proposed model. Specifically, NP-DRO sets
$\hat{\mu}$ as the Nadaraya-Watson kernel nonparametric estimator in (11) and
P-DRO sets $\hat{\mu}$ as a parametric estimator, i.e. the linear residual
model used in Example 4.8. We fix the ambiguity size $\varepsilon$ in NP-DRO,
P-DRO, and IW-DRO models to focus on the effects of the reference measures.
And for 2W-DRO, give these fixed $\varepsilon$, recall the definition of
$\kappa_{x}$ in (23), we set the two ambiguity sizes there
$(\varepsilon_{1},\varepsilon_{2})=(\mathcal{W}_{1}(\hat{\mu}_{1},\hat{\mu}_{2})-\varepsilon,\varepsilon)$
if $\kappa_{x}>0.5$ and
$(\varepsilon,\mathcal{W}_{1}(\hat{\mu}_{1},\hat{\mu}_{2})-\varepsilon)$
otherwise. Here we design
$\varepsilon_{1}+\varepsilon_{2}=\mathcal{W}_{1}(\hat{\mu}_{1},\hat{\mu}_{2})$
to ensure the feasibility of $\mathcal{A}$ from Proposition 3.1. We introduce
other benchmarks used in contextual optimization in Section 7.3.2.

###  7.2 Income Prediction

We first apply our framework to an income prediction task with the cost
$c(z;Y)=|z-Y|$, which is a special case of (10) with $b=h=1$, where $z$ and
$y$ denote the predicted and true realized income given the individual’s
feature $x$. We consider a real-world ACS Income dataset for the state of
California state (Ding et al. 2021, Liu et al. 2024), which includes 195,665
samples and categorize individuals in this dataset as either young or old
based on a threshold of 25 years of age. In the historical data $\mu_{(X,Y)}$,
we include $\alpha\times 100\%$ old people and $(1-\alpha)\times 100\%$ young
people. However in the test environment, we reverse the combination to
$\alpha\times 100\%$ old people and $(1-\alpha)\times 100\%$ young people. We
adjust the value of $\alpha$ to indicate the magnitude of covariate shifts.

##### Results.

We randomly generate 100 problem instances, each sampling a different set of
historical data. We assess the performance quality of each solution
$\hat{z}(X)$ through the metric _coefficient of prescriptiveness_ from
(Bertsimas and Kallus 2020), which is defined by:

$\text{coefficient of prescriptiveness}\coloneqq
1-\frac{\mathbb{E}_{\nu_{(X,Y)}}[c(\hat{z}(X);Y)]-\mathbb{E}_{\nu_{(X,Y)}}[c(z^{*}(X);Y)]}{\mathbb{E}_{\nu_{(X,Y)}}[c(\hat{z}^{SAA};Y)]-\mathbb{E}_{\nu_{(X,Y)}}[c(z^{*}(X);Y)]},$
(27)

where $z^{*}(X)$ is the optimal solution to (2) for each $x$ and
$\hat{z}^{SAA}=\argmin_{z\in\mathcal{Z}}\sum_{i=1}^{n}c(z;Y_{i})$. For each
term $\mathbb{E}_{\nu_{(X,Y)}}[c(z(X);Y)]$ in (27), we approximate their value
by $\frac{1}{50}\sum_{j\in[50]}c(z(\tilde{x}_{j});\tilde{y}_{j})$ with
$\\{(\tilde{x}_{j},\tilde{y}_{j})\\}_{j\in[50]}$ from the test environment.
This metric _coefficient of prescriptiveness_ evaluates the relative
performance of the model compared with the SAA solution to show the value of
covariates used in this model, where a high value implies a relatively good
model performance. We report our results in Figure 4. Across all different
scenarios, 2W-DRO has a better performance on average while the performance of
IW-DRO is slightly worse but comparable. Besides, the relative performance
gain of 2W-DRO is more stable due to a smaller standard deviation there. This
is due to constraining the distribution space from the intersection of two
ambiguity sets. Furthermore, if we only compare NP-DRO and IW-DRO in Figure 4
$(a)(b)$, when $\alpha$ is relative large, NP-DRO performs better due to more
relevant samples around while P-DRO has a better performance when $\alpha$ is
small due to the relatively small misspecification error then. When we have
100 or more samples, we can combine both powers and transform them into better
decisions.

(a) $\alpha=0.05$

(b) $\alpha=0.1$

Figure 4: Real-world income prediction from ACS Income Data, where the bar
chart shows the averaged coefficient of prescriptiveness and the vertical line
at the top of each bar shows the corresponding standard error.

###  7.3 Portfolio Allocation

We apply our framework to a mean-CVaR portfolio allocation problem. We set
$\rho=0.05$ in the objective (8) across the two following studies.

####  7.3.1 A Synthetic Study

We simulate the data generating process according to Definition 4.1. Consider
$d_{x}=10,d_{y}=d_{z}=5$ and for each problem instance, we simulate
$f(x)=0.3Bx+2\sin\|x\|_{2}\mathbf{1}_{d_{y}}$ and
$\eta\overset{d}{=}N(0,2I_{d_{y}\times d_{y}})$ with
$B\in\mathbb{R}^{d_{y}\times d_{x}}$ and each entry $B_{ij}\sim U(-2,2)$
independently. That is, we introduce a bounded misspecified linear term to
capture the misspecification error compared with P-DRO. Denote $\bm{u}_{1}$
(or $\bm{u}_{2}$ respectively) as the vector in $R^{d_{x}}$ where each
component is independently sampled from $U(0,2)$ (or $U(2,6)$ respectively)
and $\bm{\Sigma}=(AA^{\top}+A^{\top}A)/4$ as a random positive definite matrix
with $A\in\mathbb{R}^{d_{x}\times d_{x}}$ and each entry $A_{ij}\sim U(-1,1)$
independently. We design the following three scenarios to characterize the
magnitude of covariate shifts: (1) _No Covariate Shift_ :
$\mu_{x}\overset{d}{=}\nu_{x}\overset{d}{=}N(\bm{u}_{1},\bm{\Sigma})$; (2)
_Fully Covariate Shift_ :
$\mu_{x}\overset{d}{=}N(\bm{u}_{1},\bm{\Sigma}),\nu_{x}\overset{d}{=}N(\bm{u}_{2},\bm{\Sigma})$;
(3) _Mixed Covariate Shift_ :
$\mu_{x}\overset{d}{=}N(\bm{u}_{1},\bm{\Sigma}),\nu_{x}\overset{d}{=}\frac{1}{2}N(\bm{u}_{1},\bm{\Sigma})+\frac{1}{2}N(\bm{u}_{2},\bm{\Sigma})$.

##### Results.

We randomly generate 10 instances, where the model performance cost is
approximated as the average of 20 different conditional expected cost value
with each calculated as $\mathbb{E}_{\mathbb{P}_{Y|x}}[c(z,y)]$ for the
corresponding $(x,z)$. We vary two different ambiguity sizes $\varepsilon=2,5$
in Figure 5. And in that plot, across each scenario, we set the base point of
y-axis as the 2W-DRO; for other method except 2W-DRO, the y-axis value
represents the relative performance loss when switching from 2W-DRO to that
corresponding method. We find that across three cases above, 2W-DRO incurs the
smallest cost. This can be seen from Figure 5 that switching to other three
methods always results in a positive performance loss for 2W-DRO. Moreover,
the performance difference is significantly larger under scenario (2), i.e.
with large covariate shifts. The computational surrogate version, IW-DRO, also
performs relatively well with a small $\varepsilon$. However, with a large
$\varepsilon$, this IW-DRO does not perform as well as NP-DRO, potentially due
to the increased pessimism introduced by a larger ambiguity set $\mathcal{A}$.
Instead, 2W-DRO maintains the relative performance advantage over the other
methods.

(a) $\varepsilon=2$

(b) $\varepsilon=5$

Figure 5: Simulation cost difference of portfolio optimization with different
$\varepsilon$ when $n=100$ compared with 2W-DRO, where the bar chart shows the
averaged coefficient of prescriptiveness and the vertical line at the top of
each bar shows the corresponding standard error.

####  7.3.2 A Real-world Study

We conduct a real-world empirical study of a mean-CVaR portfolio optimization
task with four real-world dataset of asset returns and covariates from the
data library of Kenneth French’s
website222http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html.
Specifically, we take the most recent (i.e. last month) observed three-factor
model (Fama and French 1992) as covariates with $d_{x}=3$ in each dataset and
the associated monthly return data, each with a different number of assets
$d_{y}=6,10,25,30$ in the four real datasets.

##### Other Benchmarks.

Besides two related benchmarks in Section 7.1 above, we add different robust
contextual optimization approaches used in CVaR Portfolio Optimization in
literature. This includes: (1) Residual-DRO models from Kannan et al. (2023)
(ER-DRO-NW) with the reference measure
$\hat{\mu}_{n,x}=\frac{1}{n}\sum_{i=1}^{n}\delta_{\\{\hat{f}(x)+y_{i}-\hat{f}(x_{i})\\}}$,
where
$\hat{f}(x)=\sum_{i=1}^{n}\frac{K((x-x_{i})/h_{n})}{\sum_{j=1}^{n}K((x-x_{j})/h_{n})}y_{i}$;
(2) Joint-covariate-uncertainty DRO model (OTCMV) from Proposition 3.4 in
Nguyen et al. (2021); (3) Robust contextual optimization model (NWR) from
Proposition 4 in Srivastava et al. (2021); (4) Equally Weighted (EW) Model,
i.e. we set $z_{i}=1/d_{z}$ for all $i\in[d_{z}]$, which is popular due to its
robust empirical performance in DeMiguel et al. (2007). And for other DRO
models except 2W-DRO, we select the ambiguity size $\varepsilon$ in each model
through cross-validation. A full list of all models and their hyperparameter
configurations are shown in Appendix. For the hyperparameter in 2W-DRO, we set
their
$(\varepsilon_{1},\varepsilon_{2})=(\hat{\varepsilon}_{NP},\mathcal{W}_{1}(\mu_{n,x}^{NP},\mu_{n,x}^{P})-\hat{\varepsilon}_{NP})$
if $\kappa_{x}>0.5$ and
$(\mathcal{W}_{1}(\mu_{n,x}^{NP},\mu_{n,x}^{P})-\hat{\varepsilon}_{P},\hat{\varepsilon}_{P})$
otherwises. Here $\hat{\varepsilon}_{NP}$ and $\hat{\varepsilon}_{P}$ are the
best hyperparameters tuned from NP-DRO and P-DRO at each round.

##### Results.

We apply the “rolling-sample” approach following DeMiguel et al. (2007) on the
monthly data from July 1963 to December 2018 ($T=666$) with the estimation
window to be $M=60$ months333Rolling-sample: For data spanning $T$ months,
given month $t$, we use the dataset from months $t$ to $t+M-1$ as the input to
solve the portfolio optimization problem. Then we use the optimal weight to
compute the returns in month $t+M$. We repeat this procedure by adding the
next and dropping the earliest month until $t=T-M$. This gives us $T-M$
monthly out-of-sample returns.. We report the experimental results of
different models with the constructed returns $\\{r_{i}\\}_{i\in[T-M]}$
through three measures: Empirical CVaR (E-CVaR), Sharpe Ratio (SR), and
Certainty-equivalent Return (CEQ). E-CVaR is computed by the empirical
condition value-at-risk of the constructed returns $\\{r_{i}\\}_{i\in[T-M]}$.
SR and CEQ are computed with the formula in DeMiguel et al. (2007):

$\displaystyle\text{SR}\coloneqq\widehat{\text{mean}}(\\{r_{i}\\}_{i\in[T-M]})/\widehat{\text{std}}(\\{r_{i}\\}_{i\in[T-M]}),\quad\text{CEQ}\coloneqq\widehat{\text{mean}}(\\{r_{i}\\}_{i\in[T-M]})-(\widehat{\text{std}}(\\{r_{i}\\}_{i\in[T-M]}))^{2}.$

In Table 1, we report the model performance of methods listed above. We see
that the two models 2W-DRO and IW-DRO outperform other benchmarks in terms of
a large SR and CEQ and a small E-CVaR across the four datasets. Especially,
2W-DRO achieves the best performance in the four datasets in terms of Sharpe
Ratio and CEQ. Both of them perform well in terms of the empirical CVaR across
the four datasets. All of these suggest that the constructed portfolios from
2W-DRO and IW-DRO balance their utility and risks better in portfolio
constructions compared with all others benchmarks.

Table 1: Comparison of different policies in the portfolio allocation problem, with boldfaced values highlighting the best method (in the corresponding column) for the given metric. Results of other methods can be found in LABEL:tab:full-exp in Appendix. | EW | ER-DRO-NW | OTCMV | NWR | NP-DRO | P-DRO | IW-DRO | 2W-DRO
---|---|---|---|---|---|---|---|---
6 Portfolios Formed on Size and Book-to-Market (2 x 3)
SR | 0.2045 | 0.2264 | 0.2251 | 0.2331 | 0.2269 | 0.2285 | 0.2284 | 0.2399
E-CVaR | 0.1089 | 0.1027 | 0.0972 | 0.0951 | 0.1014 | 0.1010 | 0.0983 | 0.0901
CEQ | 0.0076 | 0.0083 | 0.0082 | 0.0086 | 0.0083 | 0.0085 | 0.0083 | 0.0088
10 Industry Portfolios
SR | 0.2203 | 0.2597 | 0.2365 | 0.2290 | 0.2545 | 0.2597 | 0.2577 | 0.2697
E-CVaR | 0.0917 | 0.0813 | 0.0844 | 0.0864 | 0.0814 | 0.0860 | 0.0822 | 0.0793
CEQ | 0.0075 | 0.0084 | 0.0076 | 0.0077 | 0.0082 | 0.0085 | 0.0084 | 0.0090
25 Portfolios Formed on Book-to-Market and Investment (5 x 5)
SR | 0.2029 | 0.2317 | 0.2417 | 0.2362 | 0.2318 | 0.2407 | 0.2418 | 0.2430
E-CVaR | 0.1121 | 0.1002 | 0.0977 | 0.0933 | 0.1005 | 0.0979 | 0.0994 | 0.1015
CEQ | 0.0078 | 0.0085 | 0.0087 | 0.0083 | 0.0084 | 0.0090 | 0.0089 | 0.0092
30 Industry Portfolios
SR | 0.1981 | 0.2297 | 0.2528 | 0.2310 | 0.2476 | 0.2565 | 0.2681 | 0.2673
E-CVaR | 0.1042 | 0.0850 | 0.0826 | 0.0820 | 0.0860 | 0.0889 | 0.0831 | 0.0856
CEQ | 0.0072 | 0.0074 | 0.0084 | 0.0074 | 0.0082 | 0.0090 | 0.0090 | 0.0095

To show that the improved performance is due to the existence of covariate
shift, next we illustrate the phenomenon that the covariate shift environment
occurs across time in the real-world datasets in Figure 6, where the red line
denotes the evolution of $\sum_{i\in[n]}K((x-x_{i})/h_{n})$, which
approximates the number of samples near the current covariate $x$, and the
blue line denotes the evolution of $\kappa_{x}$ (defined in (23)), which
controls the model specification in both IW-DRO and 2W-DRO. Recall that the
main difference between our proposed models and existing contextual-DRO models
is whether the models take the covariate shift into account. From Figure 6, we
find that severe covariate shift indeed exists across time, e.g. the global
financial crisis around 2008. Given the significant performance improvement
between our proposed model and existing DRO benchmarks, it is evident that
incorporating covariate shift into the design of the ambiguity set (i.e.
2W-DRO, IW-DRO) leads to an improved performance. And we see their superior
results are due to their ability to identify and adapt to the covariate shift
patterns in the history data.

Figure 6: Illustrations of Covariate Shift across time for 10 Industry
Portfolios set $h_{n}=0.1n^{-\frac{1}{d_{x}+d_{y}}}$ and
$K(x)=\exp(-\|x\|^{2})$. The green dashed line denotes the threshold of 0.5.

##  8 Conclusion and Future Research

We study the contextual stochastic optimization problem under covariate
shifts, a phenomenon that usually occurs in practice, and propose a new
Wasserstein distributionally robust optimization framework with the
intersection of robustness. Our model has implications both from the
computational and statistical perspectives. There are several directions to
move forward. One is to study other tractable ambiguity sets that is realistic
in the face of other types of practical distribution shifts. The other is to
investigate the performance utility and computational tradeoff when we
intersect different kinds of ambiguity sets.

N.A.

## References

  * Ban and Rudin (2019) Ban GY, Rudin C (2019) The big data newsvendor: Practical insights from machine learning. _Operations Research_ 67(1):90–108.
  * Ben-Tal et al. (2013) Ben-Tal A, Den Hertog D, De Waegenaere A, Melenberg B, Rennen G (2013) Robust solutions of optimization problems affected by uncertain probabilities. _Management Science_ 59(2):341–357.
  * Bennouna and Van Parys (2022) Bennouna A, Van Parys B (2022) Holistic robust data-driven decisions. _arXiv preprint arXiv:2207.09560_ .
  * Bertsimas et al. (2018) Bertsimas D, Gupta V, Kallus N (2018) Robust sample average approximation. _Mathematical Programming_ 171(1):217–282.
  * Bertsimas and Kallus (2020) Bertsimas D, Kallus N (2020) From predictive to prescriptive analytics. _Management Science_ 66(3):1025–1044.
  * Bertsimas and Koduri (2022) Bertsimas D, Koduri N (2022) Data-driven optimization: A reproducing kernel hilbert space approach. _Operations Research_ 70(1):454–471.
  * Bertsimas et al. (2023) Bertsimas D, McCord C, Sturt B (2023) Dynamic optimization with side information. _European Journal of Operational Research_ 304(2):634–651.
  * Bertsimas and Van Parys (2021) Bertsimas D, Van Parys B (2021) Bootstrap robust prescriptive analytics. _Mathematical Programming_ 1–40.
  * Besbes et al. (2022) Besbes O, Ma W, Mouchtaki O (2022) Beyond iid: data-driven decision-making in heterogeneous environments. _Advances in Neural Information Processing Systems_ 35:23979–23991.
  * Besbes et al. (2023) Besbes O, Ma W, Mouchtaki O (2023) From contextual data to newsvendor decisions: On the actual performance of data-driven algorithms. _arXiv preprint arXiv:2302.08424_ .
  * Blanchet et al. (2019) Blanchet J, Kang Y, Murthy K (2019) Robust wasserstein profile inference and applications to machine learning. _Journal of Applied Probability_ 56(3):830–857.
  * Blanchet and Murthy (2019) Blanchet J, Murthy K (2019) Quantifying distributional model risk via optimal transport. _Mathematics of Operations Research_ 44(2):565–600.
  * Bolley et al. (2007) Bolley F, Guillin A, Villani C (2007) Quantitative concentration inequalities for empirical measures on non-compact spaces. _Probability Theory and Related Fields_ 137(3-4):541–593.
  * Boyd and Vandenberghe (2004) Boyd SP, Vandenberghe L (2004) _Convex optimization_ (Cambridge University Press).
  * Breiman (2001) Breiman L (2001) Random forests. _Machine learning_ 45:5–32.
  * Chen et al. (2019) Chen Z, Sim M, Xu H (2019) Distributionally robust optimization with infinitely constrained ambiguity sets. _Operations Research_ 67(5):1328–1344.
  * Delage and Ye (2010) Delage E, Ye Y (2010) Distributionally robust optimization under moment uncertainty with application to data-driven problems. _Operations Research_ 58(3):595–612.
  * DeMiguel et al. (2007) DeMiguel V, Garlappi L, Uppal R (2007) Optimal versus naive diversification: How inefficient is the 1/n portfolio strategy? _The Review of Financial Studies_ 22(5):1915–1953.
  * Dietterich (2000) Dietterich TG (2000) Ensemble methods in machine learning. _International workshop on multiple classifier systems_ , 1–15 (Springer).
  * Ding et al. (2021) Ding F, Hardt M, Miller J, Schmidt L (2021) Retiring adult: New datasets for fair machine learning. _Advances in neural information processing systems_ 34:6478–6490.
  * Duchi et al. (2023) Duchi J, Hashimoto T, Namkoong H (2023) Distributionally robust losses for latent covariate mixtures. _Operations Research_ 71(2):649–664.
  * Duchi et al. (2021) Duchi JC, Glynn PW, Namkoong H (2021) Statistics of robust optimization: A generalized empirical likelihood approach. _Mathematics of Operations Research_ 46(3):946–969.
  * Dudley (1969) Dudley RM (1969) The speed of mean glivenko-cantelli convergence. _The Annals of Mathematical Statistics_ 40(1):40–50.
  * El Balghiti et al. (2019) El Balghiti O, Elmachtoub AN, Grigas P, Tewari A (2019) Generalization bounds in the predict-then-optimize framework. _Advances in neural information processing systems_ 32\.
  * Elmachtoub and Grigas (2022) Elmachtoub AN, Grigas P (2022) Smart “predict, then optimize”. _Management Science_ 68(1):9–26.
  * Esfahani and Kuhn (2018) Esfahani PM, Kuhn D (2018) Data-driven distributionally robust optimization using the wasserstein metric: Performance guarantees and tractable reformulations. _Mathematical Programming_ 171(1):115–166.
  * Esteban-Pérez and Morales (2022) Esteban-Pérez A, Morales JM (2022) Distributionally robust stochastic programs with side information based on trimmings. _Mathematical Programming_ 195(1-2):1069–1105.
  * Fama and French (1992) Fama EF, French KR (1992) The cross-section of expected stock returns. _the Journal of Finance_ 47(2):427–465.
  * Fournier and Guillin (2015) Fournier N, Guillin A (2015) On the rate of convergence in wasserstein distance of the empirical measure. _Probability Theory and Related Fields_ 162(3):707–738.
  * Freund and Schapire (1997) Freund Y, Schapire RE (1997) A decision-theoretic generalization of on-line learning and an application to boosting. _Journal of computer and system sciences_ 55(1):119–139.
  * Gao et al. (2022) Gao R, Chen X, Kleywegt AJ (2022) Wasserstein distributionally robust optimization and variation regularization. _Operations Research_ .
  * Gao and Kleywegt (2023) Gao R, Kleywegt A (2023) Distributionally robust stochastic optimization with wasserstein distance. _Mathematics of Operations Research_ 48(2):603–655.
  * Gupta and Kallus (2022) Gupta V, Kallus N (2022) Data pooling in stochastic optimization. _Management Science_ 68(3):1595–1615.
  * Györfi et al. (2006) Györfi L, Kohler M, Krzyzak A, Walk H (2006) _A distribution-free theory of nonparametric regression_ (Springer Science & Business Media).
  * Hanasusanto and Kuhn (2013) Hanasusanto GA, Kuhn D (2013) Robust data-driven dynamic programming. _Advances in Neural Information Processing Systems_ 26:827–835.
  * Hannah et al. (2010) Hannah L, Powell W, Blei D (2010) Nonparametric density estimation for stochastic optimization with an observable state variable. _Advances in Neural Information Processing Systems_ 23:820–828.
  * Hastie et al. (2009) Hastie T, Tibshirani R, Friedman JH, Friedman JH (2009) _The elements of statistical learning: data mining, inference, and prediction_ , volume 2 (Springer).
  * Hsu et al. (2012) Hsu D, Kakade SM, Zhang T (2012) Random design analysis of ridge regression. _Conference on learning theory_ , 9–1 (JMLR Workshop and Conference Proceedings).
  * Hu et al. (2018) Hu W, Niu G, Sato I, Sugiyama M (2018) Does distributionally robust supervised learning give robust classifiers? _International Conference on Machine Learning_ , 2029–2037 (PMLR).
  * Hu et al. (2022) Hu Y, Kallus N, Mao X (2022) Fast rates for contextual linear optimization. _Management Science_ 68(6):4236–4245.
  * Iyengar et al. (2023) Iyengar G, Lam H, Wang T (2023) Hedging against complexity: Distributionally robust optimization with parametric approximation. _International Conference on Artificial Intelligence and Statistics_ , 9976–10011 (PMLR).
  * Jacobs et al. (1991) Jacobs RA, Jordan MI, Nowlan SJ, Hinton GE (1991) Adaptive mixtures of local experts. _Neural computation_ 3(1):79–87.
  * Johnson-Yu et al. (2023) Johnson-Yu S, Finocchiaro J, Wang K, Vorobeychik Y, Sinha A, Taneja A, Tambe M (2023) Characterizing and improving the robustness of predict-then-optimize frameworks. _International Conference on Decision and Game Theory for Security_ , 133–152 (Springer).
  * Kallus and Mao (2023) Kallus N, Mao X (2023) Stochastic optimization forests. _Management Science_ 69(4):1975–1994.
  * Kannan et al. (2023) Kannan R, Bayraksan G, Luedtke JR (2023) Residuals-based distributionally robust optimization with covariate information. _Mathematical Programming_ 1–57.
  * Kantorovich and Rubinshtein (1958) Kantorovich LV, Rubinshtein S (1958) On a space of totally additive functions. _Vestnik of the St. Petersburg University: Mathematics_ 13(7):52–59.
  * Lakshminarayanan et al. (2017) Lakshminarayanan B, Pritzel A, Blundell C (2017) Simple and scalable predictive uncertainty estimation using deep ensembles. _Advances in neural information processing systems_ 30\.
  * Lee and Raginsky (2017) Lee J, Raginsky M (2017) Minimax statistical learning with wasserstein distances. _arXiv preprint arXiv:1705.07815_ .
  * Liu et al. (2024) Liu J, Wang T, Cui P, Namkoong H (2024) On the need for a language describing distribution shifts: Illustrations on tabular datasets. _Advances in Neural Information Processing Systems_ 36\.
  * Nadaraya (1964) Nadaraya EA (1964) On estimating regression. _Theory of Probability & Its Applications_ 9(1):141–142.
  * Nguyen et al. (2021) Nguyen VA, Zhang F, Blanchet J, Delage E, Ye Y (2021) Robustifying conditional portfolio decisions via optimal transport. _arXiv preprint arXiv:2103.16451_ .
  * Pflug and Wozabal (2007) Pflug G, Wozabal D (2007) Ambiguity in portfolio selection. _Quantitative Finance_ 7(4):435–442.
  * Pun et al. (2023) Pun CS, Wang T, Yan Z (2023) Data-driven distributionally robust cvar portfolio optimization under a regime-switching ambiguity set. _Manufacturing & Service Operations Management_ 25(5):1779–1795.
  * Qi et al. (2021) Qi M, Grigas P, Shen ZJM (2021) Integrated conditional estimation-optimization. _arXiv preprint arXiv:2110.12351_ .
  * Quinonero-Candela et al. (2008) Quinonero-Candela J, Sugiyama M, Schwaighofer A, Lawrence ND (2008) _Dataset shift in machine learning_ (Mit Press).
  * Rahimian and Mehrotra (2019) Rahimian H, Mehrotra S (2019) Distributionally robust optimization: A review. _arXiv preprint arXiv:1908.05659_ .
  * Rockafellar and Uryasev (2002) Rockafellar RT, Uryasev S (2002) Conditional value-at-risk for general loss distributions. _Journal of Banking & Finance_ 26(7):1443–1471.
  * Sadana et al. (2024) Sadana U, Chenreddy A, Delage E, Forel A, Frejinger E, Vidal T (2024) A survey of contextual optimization methods for decision-making under uncertainty. _European Journal of Operational Research_ .
  * Shafieezadeh Abadeh et al. (2015) Shafieezadeh Abadeh S, Mohajerin Esfahani PM, Kuhn D (2015) Distributionally robust logistic regression. _Advances in Neural Information Processing Systems_ 28\.
  * Shapiro et al. (2014) Shapiro A, Dentcheva D, Ruszczyński A (2014) _Lectures on stochastic programming: modeling and theory_ (SIAM).
  * Shimodaira (2000) Shimodaira H (2000) Improving predictive inference under covariate shift by weighting the log-likelihood function. _Journal of statistical planning and inference_ 90(2):227–244.
  * Srivastava et al. (2021) Srivastava PR, Wang Y, Hanasusanto GA, Ho CP (2021) On data-driven prescriptive analytics with side information: A regularized nadaraya-watson approach. _arXiv preprint arXiv:2110.04855_ .
  * Sugiyama et al. (2007) Sugiyama M, Krauledat M, Müller KR (2007) Covariate shift adaptation by importance weighted cross validation. _Journal of Machine Learning Research_ 8(5).
  * Tanoumand et al. (2023) Tanoumand N, Bodur M, Naoum-Sawaya J (2023) Data-driven distributionally robust optimization: Intersecting ambiguity sets, performance analysis and tractability .
  * Taskesen et al. (2021) Taskesen B, Yue MC, Blanchet J, Kuhn D, Nguyen VA (2021) Sequential domain adaptation by synthesizing distributionally robust experts. _International Conference on Machine Learning_ , 10162–10172 (PMLR).
  * Tsybakov (2008) Tsybakov AB (2008) _Introduction to nonparametric estimation_ (Springer Science & Business Media).
  * Van Parys et al. (2021) Van Parys BP, Esfahani PM, Kuhn D (2021) From data to decisions: Distributionally robust optimization is optimal. _Management Science_ 67(6):3387–3402.
  * Villani (2009) Villani C (2009) _Optimal transport: old and new_ , volume 338 (Springer).
  * Wainwright (2019) Wainwright MJ (2019) _High-dimensional statistics: A non-asymptotic viewpoint_ , volume 48 (Cambridge university press).
  * Wang et al. (2021) Wang T, Chen N, Wang C (2021) Distributionally robust prescriptive analytics with wasserstein distance. _arXiv preprint arXiv:2106.05724_ .
  * Watson (1964) Watson GS (1964) Smooth regression analysis. _Sankhyā: The Indian Journal of Statistics, Series A_ 359–372.
  * Weed and Bach (2019) Weed J, Bach F (2019) Sharp asymptotic and finite-sample rates of convergence of empirical measures in wasserstein distance. _Bernoulli_ 25(4A):2620–2648.
  * Wilder et al. (2019) Wilder B, Dilkina B, Tambe M (2019) Melding the data-decisions pipeline: Decision-focused learning for combinatorial optimization. _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, 1658–1665.
  * Yang et al. (2022) Yang J, Zhang L, Chen N, Gao R, Hu M (2022) Decision-making with side information: A causal transport robust approach. _Optimization Online_ .
  * Zhang et al. (2022) Zhang L, Yang J, Gao R (2022) A simple and general duality proof for wasserstein distributionally robust optimization. _arXiv preprint arXiv:2205.00362_ .
  * Zhang et al. (2024) Zhang L, Yang J, Gao R (2024) Optimal robust policy for feature-based newsvendor. _Management Science_ 70(4):2315–2329.
  * Zhuang et al. (2020) Zhuang F, Qi Z, Duan K, Xi D, Zhu Y, Zhu H, Xiong H, He Q (2020) A comprehensive survey on transfer learning. _Proceedings of the IEEE_ 109(1):43–76.

Proofs of Statements and Additional Experiments

If not clearly stated, we would use $C_{i}$ to represent constants independent
of $n,x$ and $t$ below in our proof.

##  9 Proofs in Section 3

###  9.1 Proof of Proposition 3.1.

######  Definition 9.1 (Geodesic Metric Space)

Let $(X,d)$ be a metric space. A constant speed geodesic between two points
$x_{0},x_{i}\in X$ is a continuous curve $x:[0,1]\to X$ such that for every
$s,t\in[0,1]$, $d(x_{s},x_{t})=|s-t|d(x_{0},x_{1})$. Then the metric $(X,d)$
is a geodesic metric space.

From the discussion in Section 7 from Villani (2009), we know that the space
of probability measures equipped with any Wasserstein distance
$\mathcal{W}_{p}$ is a geodesic metric space (From Section 7 in Villani
(2009), then the result follows the same structure by Proposition 5.3 in
Taskesen et al. (2021). $\hfill\square$

###  9.2 Proof of Theorem 3.2.

######  Definition 9.2 (Interchangeability Principle, from Zhang et al.
(2022))

We say a measurable function
$h:\mathcal{Y}_{1}\times\mathcal{Y}_{2}\to\mathbb{R}$ and $\mathbb{P}$
satisfies the interchangeability principle if the function
$g(\cdot)=\sup_{y\in\mathcal{Y}_{2}}h(\cdot,y)$ is measurable under
$\mathbb{P}$ and:

$\mathbb{E}_{\hat{Y}}\left[\sup_{y\in\mathcal{Y}_{2}}h(\hat{Y},y)\right]=\sup_{\gamma\in\Pi_{\mathbb{P}}}\mathbb{E}_{(\hat{Y},Y)\sim\gamma}[h(\hat{Y},Y)],$

where $\Pi_{\mathbb{P}}$ is the set of probability distributions supported on
$\mathcal{Y}_{1}\times\mathcal{Y}_{2}$ with the first marginal being
$\mathbb{P}$.

The above interchangeability condition holds for
$(\hat{Y}_{1},\ldots,\hat{Y}_{K})\sim\Gamma$ follows the discrete distribution
and our cost function.

We provide a stronger result with a generalized
$\mathcal{W}_{d}(\mu_{1},\mu_{2})=\inf_{\pi\in\Pi(\mu_{1},\mu_{2})}\mathbb{E}_{\pi}[d(X,Y)]$
as long as $d(\cdot,\cdot)$ is a distance metric. Then when
$d(X,Y)=\|X-Y\|_{1}^{p}$, the special case of one Wasserstein ball has the
following equivalence structure:

$\\{\mu:\mathcal{W}_{p}(\mu_{1},\mu_{2})\leq\varepsilon\\}=\\{\mu:\mathcal{W}_{d}(\mu_{1},\mu_{2})\leq\varepsilon^{p}\\}.$

######  Theorem 9.3 (Strong Duality of 2W-DRO)

If $\mathcal{A}=\\{\mu|\mathcal{W}_{d}(\mu,\mu_{i})\leq\varepsilon_{i},\forall
i\in[M]\\}$ is nonempty and the interchangeability condition in Definition 9.2
holds for the function
$\tilde{c}_{z}(\bm{\hat{}}Y,y)=c(z,Y)-\sum_{i}\lambda_{i}d(\hat{Y}_{i},y)$ for
any $\bm{\lambda}\geq 0,z\in\mathcal{Z}$ and $\mathbb{P}$, we have:

$\sup_{\mu\in\mathcal{A}}\mathbb{E}_{Y\sim\mu}[c(z,Y)]=\inf_{\bm{\lambda}\geq
0}\left\\{\sum_{i\in[M]}\lambda_{i}\varepsilon_{i}+\sup_{\mathbb{P}}\mathbb{E}_{\mathbb{P}}\left[\sup_{y\in\mathcal{Y}}\left\\{c(z,y)-\sum_{i=1}^{M}\lambda_{i}d(\hat{Y}_{i},y)\right\\}\right]\right\\},$
(28)

where the supremium over $\mathbb{P}$ is the random
$(\hat{Y}_{1},\ldots,\hat{Y}_{K})\sim\mathbb{P}$ is a multivariate joint
distribution with its $i$-th marginal distribution $\mathbb{P}|_{i}=\mu_{i}$,
i.e. $\hat{Y}_{i}\sim\mu_{i},\forall i\in[M]$.

We first take similar routines as Theorem 1 in Zhang et al. (2022), where we
take the Legendre transformation of
$c_{P,z}(\bm{\varepsilon})=\sup_{\mu\in\mathcal{A}}\mathbb{E}_{\mu}[c(z;Y)]$
while highlighting the dependence of $c_{P}(z)$ on the ambiguity size
$\bm{\varepsilon}=(\varepsilon_{1},\ldots,\varepsilon_{K})^{\top}$. That is:

$\displaystyle
c_{P,z}^{*}(\bm{\lambda}):=\sup_{\bm{\varepsilon}\geq\bm{0}}\\{c_{P,z}(\bm{\varepsilon})-\bm{\lambda}^{\top}\bm{\varepsilon}\\}$
$\displaystyle=\sup_{\bm{\varepsilon}\geq\bm{0}}\sup_{\mu}\left\\{\mathbb{E}_{\mu}[c(z;Y)]-\bm{\lambda}^{\top}\bm{\varepsilon}:\mathcal{W}_{d}(\mu,\mu_{i})\leq\varepsilon_{i},\forall
i\in[K]\right\\}$
$\displaystyle=\sup_{\mu}\sup_{\bm{\varepsilon}\geq\bm{0}}\left\\{\mathbb{E}_{\mu}[c(z;Y)]-\bm{\lambda}^{\top}\bm{\varepsilon}:\mathcal{W}_{d}(\mu,\mu_{i})\leq\varepsilon_{i},\forall
i\in[K]\right\\}$
$\displaystyle=\sup_{\mu}\left\\{\mathbb{E}_{Y\sim\mu}[c(z;Y)]-\sum_{i=1}^{K}\lambda_{i}\mathcal{W}_{d}(\mu,\mu_{i})\right\\}$
$\displaystyle=\sup_{\mu}\left\\{\mathbb{E}_{Y\sim\mu}[c(z;Y)]-\sum_{i=1}^{K}\lambda_{i}\inf_{\gamma_{i}\in\Pi(\mathcal{Y},\mathcal{Y})}\mathbb{E}_{(\hat{Y}_{i},Y)\sim\gamma_{i}}[d(\hat{Y}_{i},Y)]\right\\},$
$\displaystyle=\sup_{\mu,\gamma_{1},\ldots,\gamma_{K}}\mathbb{E}_{Y\sim\mu,(\hat{Y}_{i},Y)\sim\gamma_{i},\forall
i\in[K]}\left[c(z;Y)-\sum_{i\in[K]}\lambda_{i}d(\hat{Y}_{i},Y)\right]$
$\displaystyle=\sup_{(\hat{Y}_{1},\ldots,\hat{Y}_{K})\sim\mathbb{P},\gamma\in\Pi_{\mathbb{P}}}\mathbb{E}_{(\hat{Y}_{1},\ldots,\hat{Y}_{K},Y)\sim\gamma}\left[c(z;Y)-\sum_{i\in[K]}\lambda_{i}d(\hat{Y}_{i},Y)\right]$
$\displaystyle=\sup_{(\hat{Y}_{1},\ldots,\hat{Y}_{K})\sim\mathbb{P}}\sup_{\gamma\in\Pi_{\mathbb{P}}}\mathbb{E}_{(\hat{Y}_{1},\ldots,\hat{Y}_{K},Y)\sim\gamma}\left[c(z;Y)-\sum_{i\in[K]}\lambda_{i}d(\hat{Y}_{i},Y)\right],$

where the third equality follows from the minimization of $\bm{\varepsilon}$
over $\bm{\varepsilon}\geq 0$. And the fourth equality follows from the
definition of Wasserstein distance. The sixth equality follows by equivalence
of the superimum over the joint distributions between $(\mathbb{P},\gamma)$
and $(\mu,\gamma_{1},\ldots,\gamma_{K})$. On one hand, for any possible
coupling $(\hat{Y}_{1},\ldots,\hat{Y}_{K})\sim\mathbb{P}$ and
$\gamma\in\Pi_{\mathbb{P}}$, we can find one corresponding set of
$\\{\gamma_{i}\\}_{i\in[K]}$ and set the distribution $\mu:=\Pi|_{K+1}$ to
construct a feasible solution for $\mu,\gamma_{1},\ldots,\gamma_{K}$.
Therefore,

$\sup_{(\hat{Y}_{1},\ldots,\hat{Y}_{K})\sim\mathbb{P},\gamma\in\Pi_{\mathbb{P}}}\mathbb{E}_{(\hat{Y}_{1},\ldots,\hat{Y}_{K},Y)\sim\gamma}\left[\ldots\right]\leq\sup_{\mu,\gamma_{1},\ldots,\gamma_{K}}\mathbb{E}_{Y\sim\mu,(\hat{Y}_{i},Y)\sim\gamma_{i},\forall
i\in[K]}\left[\ldots\right].$

On the other hand, for any $\mu,\gamma_{1},\ldots,\gamma_{K}$ satisfying
$(\hat{Y}_{i},Y)\sim\gamma_{i},\forall i\in[K]$, we can find a (but not
unique) joint distribution $\gamma\in\Pi_{\mathbb{P}}$ with
$(\hat{Y}_{1},\ldots,\hat{Y}_{K})\sim\mathbb{P}$. Therefore,

$\sup_{(\hat{Y}_{1},\ldots,\hat{Y}_{K})\sim\mathbb{P},\gamma\in\Pi_{\mathbb{P}}}\mathbb{E}_{(\hat{Y}_{1},\ldots,\hat{Y}_{K},Y)\sim\gamma}\left[\ldots\right]\geq\sup_{\mu,\gamma_{1},\ldots,\gamma_{K}}\mathbb{E}_{Y\sim\mu,(\hat{Y}_{i},Y)\sim\gamma_{i},\forall
i\in[K]}\left[\ldots\right].$

Then following the interchangeability principle by Definition 9.2 for the
inner term
$\sup_{\gamma\in\Pi_{\mathbb{P}}}\mathbb{E}_{(\hat{Y}_{1},\ldots,\hat{Y}_{K},Y)\sim\gamma}\left[c(z;Y)-\sum_{i\in[K]}\lambda_{i}d(\hat{Y}_{i},Y)\right]$,
we have:

$c_{P,z}^{*}(\bm{\lambda})=\sup_{(\hat{Y}_{1},\ldots,\hat{Y}_{K})\sim\mathbb{P}}\mathbb{E}_{\mathbb{P}}\left[\sup_{y\in\mathcal{Y}}\left\\{c(z;y)-\sum_{i=1}^{K}\lambda_{i}d(\hat{Y}_{i},y)\right\\}\right].$

Then we apply Legendre transform back to obtain:

$c_{P,z}(\bm{\varepsilon})=\min_{\bm{\lambda}\geq\bm{0}}\\{\bm{\lambda}^{\top}\bm{\varepsilon}+c_{P,z}^{*}(\bm{\lambda})\\}=\min_{\bm{\lambda}\geq\bm{0}}\left\\{\bm{\lambda}^{\top}\bm{\varepsilon}+\sup_{\Gamma}\mathbb{E}_{\Gamma}\left[\sup_{y\in\mathcal{Y}}\left\\{c(z,y)-\sum_{i=1}^{K}\lambda_{i}d(\hat{Y}_{i},y)\right\\}\right]\right\\}.\hfill\square$

###  9.3 Proof of Theorem 3.4.

Since in this case $\mu_{1},\mu_{2}$ are both discrete distributions, the
joint distribution $\mathbb{P}$ is discrete. Therefore, we introduce variable
$\\{\gamma_{i,j}\\}_{i,j\in[n]}$ to encode the discrete weight on each point
of the joint distribution $\mathbb{P}$. Then the primal objective can be
reformulated as:

$\displaystyle c_{P}(z)=\inf_{\lambda_{1},\lambda_{2}\geq 0,s_{i,j}}{}$
$\displaystyle\lambda_{1}\varepsilon_{1}+\lambda_{2}\varepsilon_{2}+\sup_{\bm{\gamma}}\sum_{i,j\in[n]}\gamma_{i,j}s_{i,j}$
(29) s.t.
$\displaystyle\sup_{y\in\mathcal{Y}}\left\\{c_{z}(y)-\lambda_{1}\|y_{1,i}-y\|-\lambda_{2}\|y_{2,j}-y\|\right\\}\leq
s_{i,j},\forall i,j\in[n];$
$\displaystyle\sum_{j\in[n]}\gamma_{i,j}=w_{1,i},\forall i\in[n];$
$\displaystyle\sum_{i\in[n]}\gamma_{i,j}=w_{2,j},\forall j\in[n].$

We first notice the inner supremum structure of a transportation problem if we
fix $s_{i,j}$. That is the following optimization problem:

$\sup_{\bm{\gamma}}\left\\{\sum_{i,j\in[n]}\gamma_{i,j}s_{i,j},\quad\text{s.t.}\sum_{j\in[n]}\gamma_{i,j}=w_{1,i},\forall
i\in[n];\sum_{i\in[n]}\gamma_{i,j}=w_{2,j},\forall j\in[n].\right\\}$ (30)

By transforming the objective
$\sup_{\bm{\gamma}}\sum_{i,j\in[n]}\gamma_{i,j}s_{i,j}$ into
$-\inf_{\bm{\gamma}}\sum_{i,j\in[n]}\gamma_{i,j}(-s_{i,j})$, we are solving
the primal version of the transportation problem from $\mu_{1}$ to $\mu_{2}$
with the unit transportation cost $-s_{i,j}$ for each $(i,j)$-pair. The
objective of (29) is then reformulated as $\inf_{\lambda_{1},\lambda_{2}\geq
0,s_{i,j}}\left\\{\lambda_{1}\varepsilon_{1}+\lambda_{2}\varepsilon_{2}-\inf_{\bm{\gamma}}\sum_{i,j\in[n]}\gamma_{i,j}(-s_{i,j})\right\\}$.
To remove the min-max structure (i.e. the inner infimum), we take the duality
of the transportation problem of (30) as:

$\sup_{\bm{u},\bm{v}}\left\\{\sum_{i\in[n]}w_{1,i}u_{i}+\sum_{j\in[n]}w_{2,j}v_{j},\quad\text{s.t.}-(u_{i}+v_{j})\geq
s_{i,j},\forall i,j\in[n].\right\\}$

Then, we remove the inner min-max structure to reformulate the objective as
$\inf_{\lambda_{1},\lambda_{2}\geq
0,\bm{u},\bm{v},s_{i,j},\alpha_{i,k},\beta_{i,k}}\\{\lambda_{1}\varepsilon_{1}+\lambda_{2}\varepsilon_{2}-\sum_{i=1}^{n}w_{1,i}u_{i}-\sum_{j=1}^{n}w_{2,j}v_{j}\\}$.
Then we only need to handle the infinite constraints of:

$\sup_{y\in\mathcal{Y}}\left\\{c_{z}(y)-\lambda_{1}\|y_{1,i}-y\|-\lambda_{2}\|y_{2,j}-y\|\right\\}\leq
s_{i,j}.$ (31)

Like the single Wasserstein case in Esfahani and Kuhn (2018), we take the
convex conjugate function for each $c_{z,k}$ due to Assumption 3.3. For each
$k\in[K]$, the left-hand side of (31) is equivalent as:

$\displaystyle\sup_{y\in\mathcal{Y}}\left\\{c_{z,k}(y)-\lambda_{1}\|y_{1,i}-y\|-\lambda_{2}\|y_{2,j}-y\|\right\\}$
$\displaystyle=\min_{\|\alpha_{i,k}\|_{*}\leq\lambda_{1},\|\beta_{j,k}\|_{*}\leq\lambda_{2}}\sup_{y\in\mathcal{Y}}\left\\{c_{z,k}(y)-\alpha_{i,k}^{\top}y-\beta_{j,k}^{\top}y+\alpha_{i,k}^{\top}y_{1,i}+\beta_{j,k}^{\top}y_{i,j}\right\\}$
$\displaystyle=\min_{\|\alpha_{i,k}\|_{*}\leq\lambda_{1},\|\beta_{j,k}\|_{*}\leq\lambda_{2}}\sup_{y\in\mathcal{Y}}\left\\{c_{z,k}(y)-(\alpha_{i,k}+\beta_{j,k})^{\top}y\right\\}+\min_{\|\alpha_{i,k}\|_{*}\leq\lambda_{1},\|\beta_{j,k}\|_{*}\leq\lambda_{2}}\left\\{\alpha_{i,k}^{\top}y_{1,i}+\beta_{j,k}^{\top}y_{2,j}\right\\}.$

Through a convex conjugate of $c_{z,k}$, we obtain the following constraints
for (31):

$s_{i,j}\geq[-c_{z,k}+\mathbf{1}_{\mathcal{Y}}]^{*}(\alpha_{i,k}+\beta_{j,k})-\alpha_{i,k}^{\top}y_{1,i}-\beta_{j,k}^{\top}y_{2,j},\forall
i,j\in[n],k\in[K].$

Finally plugging the things together, the primal objective we obtain can be
reformulated as:

$\displaystyle c_{P}(z)=\inf_{\lambda_{1},\lambda_{2}\geq
0,\bm{u},\bm{v},s_{i,j},\alpha_{i,k},\beta_{i,k}}{}$
$\displaystyle\lambda_{1}\varepsilon_{1}+\lambda_{2}\varepsilon_{2}-\sum_{i=1}^{n}w_{1,i}u_{i}-\sum_{j=1}^{n}w_{2,j}v_{j}$
(32) s.t. $\displaystyle s_{i,j}\leq-(u_{i}+v_{j}),\forall i,j\in[n];$
$\displaystyle
s_{i,j}\geq[-c_{z,k}+\mathbf{1}_{\mathcal{Y}}]^{*}(\alpha_{i,k}+\beta_{j,k})-\alpha_{i,k}^{\top}y_{1,i}-\beta_{j,k}^{\top}y_{2,j},\forall
i,j\in[n],k\in[K];$
$\displaystyle\|\alpha_{i,k}\|_{*}\leq\lambda_{1},\|\beta_{j,k}\|_{*}\leq\lambda_{2},\forall
i,j\in[n],k\in[K].$

$\hfill\square$

Introducing the auxiliary variable
$s_{i}\geq\sup_{y\in\mathcal{Y}}[c(z,y)-\lambda|y-y_{i}|]$, then the objective
function is linear in $\lambda,\\{s_{i}\\}_{i\in[n]}$ but has infinite
constraints. We deal with $s_{i}$ as follows:

$\displaystyle s_{i}$
$\displaystyle\geq\sup_{y\in\mathcal{Y}}[c(z,y)-\lambda|y-y_{i}|]$ (33)
$\displaystyle=\sup_{y\in\mathcal{Y}}[c_{k}(z,y)-\max_{|p_{ik}|\leq\lambda}p_{ik}(y-y_{i})],\forall
i\in[n],k\in[K].$
$\displaystyle=\min_{|p_{ik}|\leq\lambda}\max_{y\in\mathcal{Y}}[c_{k}(z,y)-p_{ik}(y-y_{i})],\forall
i\in[n],k\in[K],$

where the second equality holds due to the Cauchy inequality and the fact that
$c(z,Y)=\max_{k\in[K]}c_{k}(z,Y)$. The third equality holds because of Sion’s
minimiax theorem for convex functions.

From the definition of convex conjugate, it is clear that:

$\displaystyle\sup_{y\in\mathcal{Y}}[c_{k}(z,y)-p_{ik}(y-y_{i})]$
$\displaystyle=[-c_{z,k}+\mathbb{I}_{\cdot\in\mathcal{Y}}]^{*}(p_{ik})+p_{ik}y_{i}$
$\displaystyle=\inf_{v_{ik}}\\{[-c_{z,k}]^{*}(p_{ik}-v_{ik})+\sigma_{\mathcal{Y}}^{*}(v_{ik})\\}+p_{ik}y_{i},$

The first equality holds because of conjugate of involution property in
Chapter 3 in Boyd and Vandenberghe (2004). The second equality holds because
$c_{k}(z,Y)$ is upper semi-continuous. Therefore, we obtain the result in
Proposition 3.4. $\hfill\square$

###  9.4 Proof of Corollary 3.6

Recall the uncertainty parameter
$y\in\mathcal{Y}(\subseteq\mathbb{R}^{d_{y}}):=\\{y:Ey\leq f\\}$ and
$c_{z,k}(y)=a_{z,k}^{\top}y+b_{z,k}$, with $a_{z,k},b_{z,k}$ being linear in
$z$. Firstly, from Theorem 3.4 and compare the result in Corollary 3.6, we
only need to reformulate the following sets of constraints:

$s_{i,j}\geq[-c_{z,k}+\sigma_{\mathcal{Y}}]^{*}(\alpha_{i,k}+\beta_{j,k})-\alpha_{i,k}^{\top}y_{1,i}-\beta_{j,k}^{\top}y_{2,j},\forall
i,j\in[n],k\in[K];$

into the following set of constraints:

$\displaystyle s_{i,j}$ $\displaystyle\geq
b_{z,k}-(\alpha_{i,k}^{\top}y_{1,i}+\beta_{j,k}^{\top}y_{2,j})+\gamma_{i,j,k}^{\top}f,\forall
i,j\in[n],k\in[K];$ $\displaystyle\alpha_{i,k}+\beta_{j,k}$
$\displaystyle=E^{\top}\gamma_{i,j,k}-a_{z,k},\forall i,j\in[n],k\in[K].$

We focus on the term
$[-c_{z,k}+\sigma_{\mathcal{Y}}]^{*}(\alpha_{i,k}+\beta_{j,k})$ in Theorem
3.4. Specifically, for the two components of the left-hand side, by strong
duality and the definition of conjugate function, we have:

$\sigma_{\mathcal{Y}}(v_{ik})=\sup_{y}\\{v_{ik}^{\top}y,~{}\text{s.t.}~{}Ey\leq
f\\}=\inf_{\gamma_{ik}\geq
0}\\{\gamma_{ik}^{\top}f,~{}\text{s.t.}~{}E^{\top}\gamma_{ik}=v_{ik}\\}.$
$[-c_{z,k}]^{*}(\alpha_{i,k}+\beta_{j,k})=[-a_{z.k}]^{*}(\alpha_{i,k}+\beta_{j,k})=\begin{cases}b_{z,k},&\alpha_{i,k}+\beta_{j,k}=-a_{z,k},\\\
\infty,&\alpha_{i,k}+\beta_{j,k}\not=-a_{z,k}\end{cases}$

Through the property of convex conjugate, we plug those above into (32) and
obtain the result. $\hfill\square$

##  10 Proofs for Section 4.1

We prove the measure concentration result depending on kernels with different
levels when the response variable $y$ is univariate. The multivariate case
follows similarly by replacing the dependence in terms of $t$ from
$O(\exp(-t^{2}))$ to $O(\exp(-t^{d_{y}})$. And all the other proofs would
follow directly. Before delving into details of measure concentration results,
we first represent some technical lemmas.

###  10.1 Technical Tools

We show in the following that the neighborhood
$N_{n}^{r}=O_{p}(c_{x}nh_{n}^{d_{x}})$ when there is no case of covariate
shift, which also plays an important role under distribution shift. We first
restate the following multiplicative Chernoff bound:

######  Lemma 10.1 (Multiplicative Chernoff Bound)

Suppose $X_{1},\ldots,X_{n}$ are independent and identically distributed
random variables taking values in $[0,B]$. Denote
$\bar{X}=\frac{1}{n}\sum_{i=1}^{n}X_{i}$ and $\mu=\mathbb{E}[X]$. Then for
$0<\delta<1$, we have:

$\mathbb{P}(\bar{X}\leq(1-\delta)\mu)\leq\exp(-\delta^{2}\mu
n/2),\quad\mathbb{P}(|\bar{X}-\mu|\geq\delta\mu)\leq 2\exp(-\delta^{2}\mu
n/3),\quad\forall\delta\in(0,1).$

######  Lemma 10.2 (Nature and Impossibilities of Local Estimators)

Under $0<\mu_{\min}\leq\mu_{x}(\cdot)\leq\mu_{\max}<\infty$, we have:

$\displaystyle\mathbb{E}[N_{n}^{r}]$
$\displaystyle\in[C_{1}\mu_{\min}nh_{n}^{d_{x}},\max\\{C_{2}\mu_{\max}nh_{n}^{d_{x}},n\\}],$
$\displaystyle\lim_{n\to\infty}\frac{\mathbb{E}[N_{n}^{r}]}{\mu_{x}nh_{n}^{d_{x}}}$
$\displaystyle=C\in(0,\infty),$
$\displaystyle\mathbb{P}(C_{3}nh_{n}^{d_{x}}\leq N_{n}^{r}\leq
C_{4}nh_{n}^{d_{x}})$ $\displaystyle\leq\exp(-C_{5}\mu_{\min}nh_{n}^{d_{x}}),$

where $C_{1},C_{2},C_{3},C_{4},C_{5}$ are independent with $n,x$.

###### Proof 10.3

Proof of Lemma 10.2.In this lemma we would build some results for the numbers
of nearby samples, which would be illustrated in the following non-parametric
estimator results. We would provide a lower probabilistic bound for
$N_{n}^{r}$. Note that the volume of the intersection of the ball
$\\{x:\|x-x_{i}\|_{2}\leq rh_{n}\\}$ and $[0,1]^{d_{x}}$ is of the order
$\min\left\\{2^{-{d_{x}}}C_{d_{x}}\pi
r^{d_{x}}h_{n}^{d_{x}},2^{-d_{x}}C_{d_{x}}\pi\right\\}$, where
$C_{d}\pi^{\frac{d}{2}}R^{d}$ is the volume of a $d$-dimensional ball444In
fact, $C_{d}=\frac{1}{\Gamma(\frac{d}{2}+1)}$.. Here the worst case is taken
when $x$ is in the corner of $[0,1]^{d_{x}}$. Since the other quantity is
independent of $t$, $x$ and $n$, we simply use $C_{3}h_{n}^{d_{x}}$. By
Assumption 4, we have:

$\mu_{\min}Vol(\\{x:\|x-x_{i}\|_{2}\leq
rh_{n}\\})\leq\mathbb{P}(x_{i}:\|x-x_{i}\|\leq
rh_{n})\leq\mu_{\max}Vol(\\{x:\|x-x_{i}\|_{2}\leq rh_{n}\\}).$

Furthermore, we have

$\lim_{n\to\infty}\frac{\mathbb{P}(x_{i}\in[0,1]^{d_{x}}:\|x-x_{i}\|\leq
rh_{n})}{\mu_{x}rh_{n}^{d_{x}}}=1.$

Since $N_{n}^{r}$ is the sum of $n$ i.i.d. Bernoulli random variables with
success probability $\mathbb{P}(\|x-x_{i}\|_{2}\leq rh_{n})$, we have:
$\mathbb{E}[N_{n}^{r}]=n\mathbb{P}(x_{i}:\|x-x_{i}\|\leq
rh_{n})\in[C_{1}\mu_{\min}nh_{n}^{d_{x}},C_{2}\mu_{\max}nh_{n}^{d_{x}}]$ for
some $C_{1},C_{2}$ independent with $n,x$. And
$\lim_{n\to\infty}\mathbb{E}[N_{n}^{r}]=C\mu_{x}nh_{n}^{d_{x}}$, where $C$ is
independent with $n,x$.

And for the concentration,
$\text{Var}(N_{n}^{r})=\mathbb{E}[N_{n}^{r}](1-\mathbb{E}[N_{n}^{r}]/n)$. From
Lemma 10.1, we would obtain:

$\displaystyle\mathbb{P}(N_{n}^{r}\leq C_{4}nh_{n}^{d_{x}}/2)$
$\displaystyle\leq\mathbb{P}\left(\frac{1}{n}(N_{n}^{r}-\mathbb{E}[N_{n}^{r}])\leq-\frac{1}{2n}\mathbb{E}[N_{n}^{r}]\right)\leq\exp\left(-C_{6}n\mu_{\min}h_{n}^{d_{x}}\right).$
(34)

On the other side, we would obtain one lower bound argument of $N_{n}^{r}$.
$\hfill\square$

The following result is from Theorem 5.10 in Kantorovich and Rubinshtein
(1958):

######  Lemma 10.4 (Dual of Wasserstein Distance)

Dual of type-$p$ $(1\leq p<\infty)$ Wasserstein distance between two
distributions $\mu$ and $\nu$ supported on $\Xi$ is computed by:

$\mathcal{W}_{p}^{p}(\mu,\nu)=\sup_{u,v\in
L^{1}(\Xi)}\left\\{\int_{\Xi}u(\xi)\mu(d\xi)+\int_{\Xi}v(y)\nu(dy):u(\xi)+v(y)\leq\|\xi-y\|_{p}^{p},\forall\xi,y\in\Xi\right\\}.$
(35)

The key step in proving the result of boxed kernel is from the following
steps, which helps to reduce the boxed kernel to the case of naive kernel.

###### Lemma 10.5

Suppose $\mu_{n,x}=\sum_{i=1}^{n}w_{i}\delta_{f(x)+\eta_{i}}$ with $w_{i}\geq
0,\forall i\in[n];\sum_{i\in[n]}w_{i}=1$,
$\mu_{n,x}^{avg}=\sum_{i=1}^{n}\frac{1}{n}\delta_{f(x)+\eta_{i}}$, and
$\mu_{Y|x}\sim f(x)+\eta$, where $\\{\eta_{i}\\}_{i\in[n]}$ are i.i.d.
sampled, then we have:

$(n\min_{i\in[n]}w_{i})^{\frac{1}{p}}\mathcal{W}_{p}(\mu_{Y|x},\mu_{n,x}^{avg})\leq\mathcal{W}_{p}(\mu_{Y|x},\mu_{n,x})\leq(n\max_{i\in[n]}w_{i})^{\frac{1}{p}}\mathcal{W}_{p}(\mu_{Y|x},\mu_{n,x}^{avg})$