A Data-Driven Distributionally Robust Bound on the Expected Optimal Value of Uncertain Mixed 0-1 Linear Programming
Guanglin Xu Samuel Burer

The Young Substellar Companion ROXs 12 B: 
 
Near-Infrared Spectrum, System Architecture, and Spin-Orbit Misalignment[Some 
of the data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership 
among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. 
The Observatory was made possible by the generous financial support of the W.M. Keck Foundation.]
    Lucas A. Cieza
    December 30, 2023
========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





This paper studies the expected optimal value of a mixed 0-1 programming 
problem with uncertain objective coefficients following a joint distribution. We assume that
the true distribution is not known exactly, but a set of independent samples can be
observed. Using the Wasserstein metric, we construct an ambiguity set centered at the
empirical distribution from the observed samples and containing the true distribution with 
a high statistical guarantee. The problem
of interest is to investigate the bound on the expected optimal value over the Wasserstein
ambiguity set. Under standard assumptions, we reformulate the problem into a copositive
program, which naturally leads to a tractable semidefinite-based approximation. We compare 
our approach with a moment-based approach from the literature on three applications.
Numerical results illustrate the effectiveness of our approach. 




Keywords: Distributionally robust optimization; Wasserstein metric; 
copositive programming; semidefinite programming








§ INTRODUCTION
 


We consider the following uncertain mixed 0-1
linear programming problem:

    v(ξ) := max{ (Fξ)^Tx : [     Ax = b, x ≥ 0; x_j ∈{0,1} ∀  j ∈ ]}

where A ∈^m × n, F ∈^n × k,
and b ∈^m are the problem data,  x ∈^n_+ is the 
vector of decision variables, ⊆{ 1, …, n }
is an index set of binary variables, and the objective coefficients are linear in
the random vector ξ∈^k via F.
Problem (<ref>) entails two extreme
classes of programs: if = ∅, then
(<ref>) represents the regular linear program
with uncertain objective coefficients;
if = { 1, …, n }, then (<ref>)
represents the regular binary program with uncertain
coefficients. In general, problem (<ref>) is NP-hard
<cit.>. 

The optimal value v(ξ) is a random variable as ξ
is a random vector. We assume 
that ξ follows a multivariate distribution
 supported on a nonempty set Ξ⊆^k,
which is, in particular, defined as a slice of a closed, convex, 
full-dimensional cone Ξ⊆_+ ×^k-1:

    Ξ := {ξ∈Ξ: e_1^Tξ =  ξ_1 = 1 },

where e_1 is the first standard basis vector in ^k. In
words, Ξ is the homogenization of Ξ. We
choose this homogenized version for notational convenience.
Note that it, in fact, enables us to model affine effects of the uncertain
parameters in (<ref>). 

The expected optimal value of (<ref>), 
denoted by v_, is defined as

    v_ := _ [v(ξ)] = ∫_Ξ v(ξ)   d ( ξ).

The problem of computing v_ has been extensively studied in
the literature. Hagstrom <cit.>
showed that computing v_ for the longest path problem over a
directed acyclic graph is #𝒫-complete even if the arc lengths
are each independently distributed and restricted to taking two possible
values. Aldous <cit.> studied a linear assignment
problem with random cost coefficients following either an independent
uniform distribution on [0,1] or an exponential distribution with
parameter 1 and proved that the asymptotic value of v_
approaches π^2/6 as the number of assignments goes
to infinity. For additional studies, see <cit.>.

In practice, it is difficult or impossible to know 
completely, and computing v_ is thus not well defined 
in this situation. An alternative is to construct an ambiguity set, denoted
by , that contains a family of distributions supported on Ξ
and consistent with any known properties of .
Ideally, the ambiguity set will possess some statistical guarantee, 
e.g., the probability that ∈ will be at least 1 - β, where
β is the significance level. In analogy with v_, 
we define v_



for any ∈. Then, we are interested in computing 
the maximum expected optimal value v_
over the ambiguity set :

    v_^+ :=  sup_∈ v_.

Note that, when the probability of ∈ is at least
1 - β, the probability of v_≤ v_^+ is
at least 1 - β. 

There are three main issues (somehow conflicting) regarding 
the computation of v_^+. First, one would like an ambiguity 
set  with a high statistical guarantee to contain the true distribution . 
In this way, the computed v_^+ will be an upper bound 
on v_ with a high confidence level. (We will introduce several 
approaches in the following paragraph.) Second, one would like
v_^+ to be tight in the sense that it is as close to
v_ as possible. Generally, if  enforces more information 
about , then v_^+ will be closer to v_. Finally, the third 
concern is the complexity of the resulting optimization problem, 
i.e., whether the problem can be solved in polynomial time.



Bertsimas et al. <cit.> constructed moment ambiguity sets
using the first two marginal moments of each ξ_i. Denote
the first and second of each uncertain parameter by μ_i and σ_i
respectively. They computed v_^+ over all 
joint distributions sharing the same first two marginal moments and proved 
polynomial-time computability if the corresponding deterministic 
problem is solvable in polynomial time. However, the computed bound
may not be tight with respect to v_ since the marginal-moment 
model does not capture the dependence of the random variables. 
In a closely related direction, Natarajan et al. <cit.>
proposed an ambiguity set that was constructed from the
known marginal distributions of each random variable ξ_i,
and they computed v_^+ by solving a concave
maximization problem. As an extension to the marginal moment-based
approach, Natarajan et al. <cit.> proposed
a cross-moment model that was based on an ambiguity
set constructed using both marginal and cross moments.
Compared to the marginal-moment approach, the cross-moment approach 
has tighter upper bounds as the model captures the dependence of the random variables.
However, computing the bound requires solving a completely
positive program, which itself can only be approximated in general. 
Thus, the authors proposed semidefinite programming (SDP) 
relaxations to approximate v_^+.  




Moment-based ambiguity sets are also used prominently
in a parallel vein of research, called distributionally robust
optimization (DRO); see <cit.>. 
The popularity of the moment-based approach is mainly due to the fact that
it often leads to tractable optimization problems and relatively
simple models. Its weakness, however, is that moment-based sets
are not guaranteed to converge to the true distribution  when 
the sample size increases to infinity, even though the
estimations of the first and second moments are themselves guaranteed
to converge. 

As an attractive alternative to moment-based ambiguity sets, 
distance-based ambiguity sets haven been proposed in recent years. 
This approach defines  as a ball in the space of probability distributions equipped with a distance
measure, and the center of the ball is typically the empirical distribution
derived from a series of independent realizations of the random vector
ξ. The key ingredient of this approach is the distance function.
Classical distance functions include the Kullback-Leibler
divergence <cit.>,
the ϕ-divergence <cit.>,
the Prohorov metric <cit.>, empirical Burg-entropy
divergence balls <cit.>, and the 
Wasserstein metric <cit.>.

In this paper, we apply the Wasserstein metric to construct a data-driven
ambiguity set  centered at the empirical distribution _N
derived from N independent observations of ξ. This approach has
several benefits. The conservativeness of the ambiguity set
can be controlled by tuning a single parameter, the radius of the
Wasserstein ball; we will discuss this parameter in detail
in Section <ref>. Also, under mild conditions on ,
the Wasserstein ambiguity provides a natural confidence set for . 
Specifically, the Wasserstein ball around the empirical distribution on 
N independent identical samples contains  with confidence 1-β 
if its radius exceeds an explicit threshold ϵ_N(β) that can be computed 
via a closed form equation <cit.>. 
We then formulate v_^+ in (<ref>) over the constructed 
Wasserstein ambiguity set. That is, we model the maximum value of 
v_ over the ambiguity set  constructed by the Wasserstein 
metric. In Section <ref>, we reformulate problem (<ref>) 
into a copositive problem under some standard assumptions. As the copositive 
reformulation is computationally intractable, we apply a standard approach 
based on semidefinite programming techniques to approximate v_^+ 
from above. In Section <ref>, we numerically verify our approach 
on three applications from the literature. In particular, we compare our 
approach with the moment-based approach proposed in 
<cit.>. We have several important 
observations from the experimental results. First, we find that the gaps 
between the bound from our semidefinite programs 
and the true expected optimal value becomes narrower as the sample size 
increases. However, the moment-based bound remains the same regardless
of the increase in the sample size. Second, we observe that our bound 
converges to the true expected optimal value on the first two applications
where the underlying deterministic problems are linear programs. Although
our bound on the third application is not able to converge to the true
expected optimal value, it is tighter than the moment-based bound after
the sample size increases to a certain level.  We conclude our research 
and discuss some future directions in Section <ref>. 

We point out some similarities of our paper to a recent technical
report by Hanasusanto and Kuhn <cit.>. In their
report, they proposed a Wasserstein-metric ambiguity set for a two-stage
DRO problem. In particular, they applied copositive programming techniques
to reformulate the second-stage worst-case value function, which is essentially
a max-min optimization problem, while we use copositive techniques to reformulate
a max-max optimization problem; see (<ref>). Furthermore, 
they directly used a hierarchy schema to approximate the copositive cones,
while we derive natural SDP approximations based on the copositive 
reformulation. Note that their hierarchy of approximations lead to SDP 
approximations as well. Finally, they developed an approach to derive an 
empirical Wasserstein radius, which is in spirit similar to our approach 
in this paper. 




 §.§ Notation, terminology, and basic techniques

 
We denote by ^n the n-dimensional Euclidean space
and by _+^n the nonnegative orthant in ^n.
For a scalar p ≥ 1, the p-norm of z ∈^n is
defined z_p := (∑_i=1^n |z_i|^p)^1/p, e.g., z_1 =
∑_i=1^n |z_i|. We will drop the subscript for the 2-norm, i.e.,
z := z_2. For v,w ∈^n, the inner product of v and
w is denoted by v^Tw := ∑_i=1^n v_iw_i. For the 
specific dimensions k and n of the problem in this paper, 
we denote by e_i the i-th standard basis vector in ^k, 
and similarly, denote by f_j the j-th standard basis vector 
in ^n. We will also define g_1 := e_1  0∈^k+n. 
We denote by δ_ξ the Dirac distribution concentrating unit mass
at ξ∈^k. For any N ∈ℕ, 
we define [N] := {1, …, N}. 








Let ^m × n denote the space of real m × n
matrices, and A∙ B := trace(A^TB) denote the 
trace of the inner product of two matrices A, B ∈^m × n. 
We denote by ^n the space of n × n symmetric matrices, 
and for X ∈^n, X ≽ 0 represents that X is positive 
semidefinite. In addition, we denote by (X) the vector containing 
the diagonal entries of X, and denote by (v) the diagonal 
matrix with vector v along its diagonal. I ∈^n denotes 
the identity matrix. 

Finally, letting K⊆^n be a closed,
convex cone, and K^* be its dual cone, we give
a brief introduction to copositive programming with
respect to the cone K. The copositive cone with respect
to K is defined as

    ( K) := { M ∈^n: x^T M x ≥ 0  ∀ x ∈ K},

and its dual cone, the completely positive cone with respect to
K, is given as

    ( K) := { X ∈^n: X = ∑_i x^i (x^i)^T,  x^i ∈ K},

where the summation over i is finite but its cardinality is
unspecified. The term copositive programming refers to linear
optimization over ( K) or, via duality, linear optimization
over ( K). In fact, these problems are sometimes called
generalized copositive programming or set-semidefinite
optimization <cit.> in contrast
with the standard case K = _+^n. In this paper, we work with
generalized copositive programming, although we use the shorter phrase
for convenience.






§ A WASSERSTEIN-BASED AMBIGUITY SET
 

In this section, we define the Wasserstein metric and discuss 
a standard method to construct a Wasserstein-based ambiguity set.
Using this ambiguity set, we fully specify problem (<ref>). 

Denote by Θ_N := {ξ̂^1, …, ξ̂^N} 
the set of N independent samples of ξ 
governed by . The uniform empirical 
distribution based on Θ_N is 
_N := 1N ∑_i=1^N δ_ξ̂^i
where δ_ζ is the Dirac distribution concentrating
unit mass at ζ∈^k.


Let ^2(Ξ) be the set of all probability distributions
 that are supported on Ξ and that satisfy _[ξ 
- ξ'^2] = ∫_Ξξ - ξ'^2 d (ξ) < ∞ where 
ξ' ∈Ξ is some reference point, e.g., ξ' = ξ̂^i for some i ∈ [N]. 




The 2-Wasserstein distance between any , ' ∈ M^2(Ξ)
is  

    W^2(, ') := inf{ (∫_Ξ^2ξ - ξ'^2  Π (dξ, dξ')  )^1/2 : 
    [ Π is a joint distribution of ξ and ξ';   with marginals  and ', respectively ]}.

 


The Wasserstein distance is essentially the minimum cost of redistributing mass
from  to '. It is also called the “earth mover's distance" in
the community of computer science; see <cit.>.
In fact, the Wasserstein distance between two discrete distributions with 
a finite number of positive masses corresponds to a transportation planning 
problem in finite dimensions. 


Example <ref> illustrates the Wasserstein distance between two
discrete distributions.
 
Consider two discrete distributions: := ∑_i=1^Mq_i δ_ξ_i
and ' := ∑_j=1^M' q'_jδ_ξ'_j where q_i ≥ 0   i = 1, …, M, 
q'_j ≥ 0  j = 1, …, M', and ∑_i=1^M q_i = ∑_j=1^M'q'_j = 1.
Define c_ij = ξ_i - ξ'_j^2  ∀  i = 1, …, M   j = 1, …, M'.  
Then, the 2-Wasserstein distance between  and ' equals the square root
of the optimal value of the following transportation planning problem:

    [                                   min_π               ∑_i=1^M ∑_j=1^M' c_ijπ_ij;                                              ∑_j=1^M'π_ij = q_i  ∀  i = 1, …, M;                                             ∑_i=1^Mπ_ij = q'_j  ∀  j = 1, …, M';                                         π_ij≥ 0  ∀  i = 1, …, M,  j = 1, …, M', ]

where π is the joint distribution of ξ and ξ' with marginals of  and '
and π is the matrix variable in this optimization problem. 



With this setting, our ambiguity set contains a family of 
distributions that are close to _N with respect to the 
Wasserstein metric. In particular, we define our ambiguity set  
as a 2-Wasserstein ball of radius ϵ that is centered at the 
uniform empirical distribution _N:

    (_N, ϵ) := {∈^2(Ξ) :  W^2(, _N) ≤ϵ}.




The reader is referred to <cit.> for 
the general case of ^r(Ξ) and W^r(, ') for any
r ≥ 1. We use the 2-Wasserstein distance in this paper for  
two reasons. First, the Euclidean distance is one of 
the most popular distances considered in the relevant literature;
see <cit.>. 
Second, we will find that problem (<ref>) with an ambiguity 
set based on the 2-Wasserstein distance can be reformulated into 
a copositive program; see Section <ref>. 




Then, we replace the generic ambiguity set  with the Wasserstein 
ball (_N, ϵ) in problem (<ref>) to compute
a data-driven upper bound: 



    [                         v_(_N, ϵ)^+ =                        sup_Π,  ∈^2(Ξ)                      ∫_Ξ  v(ξ)  d (ξ);                                                                                    ∫_Ξ^2ξ - ξ'^2  Π(dξ, dξ') ≤ϵ^2;                                                                             Π is a joint distribution of ξ and ξ';                                                                             with marginals  and _N, respectively. ]


We next close this subsection by making some remarks.
First, the Wasserstein ball radius in problem (<ref>)
controls the conservatism of the optimal value. A larger radius 
is more likely to contain the true distribution and thus a more likely 
valid upper bound on v_, but even if it is valid, it could be a 
weaker upper bound. Therefore, it is crucial to choose an 
appropriate radius for the Wasserstein ball. Second, the 
Kullback-Leibler divergence ball is also considered in recent research; 
see <cit.>. However, in the case of
our discrete empirical distribution, the Kullback-Leibler divergence ball is a 
singleton containing only the empirical distribution itself, with probability 
one. Third, the ambiguity sets constructed by goodness-of-fit 
tests in <cit.>
also possess statistical guarantees, however, they often lead to
complicated and intractable optimization problems for the case of 
high-dimensional uncertain parameters. 
 











 §.§ An empirical Wasserstein radius

The papers <cit.> 
present a theoretical radius ϵ_N(β) for 
datasets of size N, which guarantees a desired confidence 
level 1 - β for ∈(_N, ϵ_N(β)) 
under the following standard assumption on :
[Light-tailed distribution] There exists an
exponent a > 1 such that 

    _ [exp(ξ^a)] = ∫_Ξexp(ξ^a)d (ξ) <  ∞.


Note that ϵ_N(β) depends on N 
and β. However, ϵ_N(β) is known to be 
conservative in practice; see <cit.> for 
example. In other words, (_N, ϵ_N(β)) 
might contain significantly more irrelevant distributions so 
that the computed v_(_N,  ϵ_N(β))^+ 
is significantly larger than v_. So, we propose a 
procedure to derive an empirical radius that provides 
a desired confidence level 1 - β but is much smaller 
than ϵ_N(β). Our approach is based on 
the data set Θ_N. In particular, we 
apply a procedure, similar to cross validation in spirit, 
that computes an empirical confidence level (between 
0 and 1) for a given radius ϵ; see details 
in the next paragraphs. Our procedure guarantees that 
a larger radius leads to a higher confidence level. 
Therefore, by iteratively testing different ϵ, 
we can find a radius with a desired confidence level based 
on the data set Θ_N. Although the derived 
ϵ(Θ_N,  β) depends on the 
data set Θ_N, our experimental results in 
Section <ref> indicate that it can be used 
for other datasets of  the same sample size. We will 
show the numerical evidence in Section <ref>. 
Our approach is also similar in spirit to the one used 
in <cit.>.

Our procedure requires an oracle to compute (or approximate) 
v_^+. Later in Section <ref>, we will propose 
a specific approximation; see (<ref>). Assume also that, 
in addition to the dataset Θ_N, we predetermine a 
set ℰ containing a large, 
yet finite, number of candidate radii ϵ. We randomly divide 
Θ_N into training and validation datasets K times.  
We enforce the same dataset size denoted by N_T on each of the 
K training datasets. 







Next, for each ϵ∈ℰ, we derive an empirical 
probability based on the following procedure: (i) we use each of 
the K training datasets to approximate 
v_(_N_T,  ϵ)^+ with a value called 
v_WB(ϵ) by calling the oracle, where _N_T
represents the empirical distribution from the training set; (ii) we then use 
the corresponding K validation datasets to simulate the expected 
optimal values denoted by v_SB [This process 
is to solve a linear program or integer program corresponding to 
each sample in the validation dataset and then to take the average 
of the optimal values.]; and (iii) we finally compute 
the percentage of the K instances where v_WB(ϵ) 
≥ v_SB. Let us call this empirical probability as the 
empirical confidence level. Thus, the empirical confidence level 
can roughly approximate the confidence level that the underlying 
distribution is contained in the Wasserstein-based ambiguity 
set with the radius ϵ. Note that the percentage 
computed is non-decreasing in ϵ and equal to 1 for 
some large ϵ_0. Therefore, the set containing all the 
empirical confidence levels is essentially an empirical cumulative distribution. 
Then, given a desired confidence level, we can 
choose a corresponding empirical radius ϵ∈ℰ. 
The numerical results in Section <ref> indicate that
our choices of ϵ indeed return the desired confidence levels.  
We specify the above procedure in Algorithm <ref>.





§ PROBLEM REFORMULATION AND TRACTABLE BOUND
 

In this section, we propose a copositive programming reformulation
for problem (<ref>) under some mild assumptions. 
As copositive programs are computationally intractable, we then propose
semidefinite-based relaxations for the purposes of computation. 



Let us first define the feasible set for x ∈^n  in (<ref>) 
as follows:  

    := { x ∈^n :  [     Ax = b, x ≥ 0; x_j ∈{0,1} ∀  j ∈ ]}.

We now introduce the following standard assumptions: 

 
The set ⊆^n is nonempty and bounded. 


 
Ax = b,  x ≥ 0  0 ≤ x_j ≤ 1  ∀  j ∈. 

Assumption <ref> can be easily enforced. 
For example (see also <cit.>, <cit.>), 
if = ∅, then the assumption is redundant; if problem (<ref>)
is derived from the network flow problems, for instance
the longest path problem on a directed acyclic graph, then
Assumption <ref> is implied from the network
flow constraints; if  is a nonempty set and the 
assumption is not implied by the constraints, we can add
constraints x_j + s_j = 1,  s_j ≥ 0  ∀  j ∈.

 
The support set Ξ⊆^k is convex, closed, 
and computationally tractable.


For example, Ξ􏰐 could be represented using a polynomial 
number of linear, second-order-cone, and semidefinite inequalities. 
In particular, the set Ξ possesses a polynomial-time separation 
oracle <cit.>.

  
Ξ is bounded. 


By Assumption <ref>, we know that
v(ξ) is finite and attainable for any ξ∈Ξ.



Note that under Assumptions 
<ref>-<ref>, 
v_^+ is finite and attainable and thus we can 
replace sup with max in (<ref>) under 
these conditions. Assumption <ref> 
could be merged with Assumption <ref>, but 
it is stated separately to highlight its role in proving the 
exactness of the copositive programming reformulation 
below. 




 §.§ A copositive reformulation
 





We reformulate problem (<ref>) via conic programming duality
theory and probability theory. We introduce a useful result from the literature
as follows.

 
v_(_N, ϵ)^+ equals the optimal value of

    [                                    sup          1/N∑_i=1^N ∫_Ξ  v(ξ)  d _i(ξ);                                        1/N∑_i=1^N ∫_Ξξ - ξ̂_i^2  d _i(ξ) ≤ϵ^2;                                                    _i ∈ M^2(Ξ)    ∀  i ∈ [N], ]

where _i represents the distribution of ξ conditional on ξ' = ξ̂^i for all i ∈ [N].  


As _i represents the distribution of ξ conditional on 
ξ' = ξ̂^i, the joint probability Π in problem (<ref>) 
can be decomposed as Π = 1/N∑_i ∈ [N]_i 
by the law of total probability. Thus, the optimal value of (<ref>)
coincides with v_(_N, ϵ)^+, which completes the proof.
 

We next provide a copositive programming reformulation
for problem (<ref>). As the first step, we use
a standard duality argument to write the dual of (<ref>)
(see also <cit.>):

    v_(_N, ϵ)^+ =    sup__i ∈^2(Ξ)inf_λ≥ 01/N∑_i=1^N ∫_Ξ v(ξ)  d _i(ξ) + λ (ϵ^2 - 1/N∑_i=1^N ∫_Ξξ - ξ̂_i^2  d _i(ξ)  ) 
    ≤   inf_λ≥ 0sup__i ∈^2(Ξ)λ ϵ^2 + 1/N∑_i=1^N ∫_Ξ (v(ξ) - λξ - ξ̂_i^2)  d _i(ξ) 
    
       =    inf_λ≥ 0λ ϵ^2 + 1/N∑_i=1^N sup_ξ∈Ξ (v(ξ) - λξ - ξ̂_i^2),

where (<ref>) follows from the max-min inequality, 
while equation (<ref>) follows from the
fact that M^2(Ξ) contains all the Dirac distributions supported on Ξ. 


















By Assumption <ref>, v(ξ) is finite for all ξ∈Ξ.
Then, the inequality in (<ref>) becomes an equality for any
ϵ > 0 due to a straightforward generalization of a strong duality 
result for moment problems in Proposition 3.4 in <cit.>; see
also Theorem 1 in <cit.> and Lemma 7 in 
<cit.>. 
By introducing auxiliary variables s_i, the minimization problem in (<ref>)
is equivalent to 

    [                                    v_(_N, ϵ)^+ =                                     inf_λ,   s_i                           λ ϵ^2 + 1/N∑_i=1^N s_i;                                                                                                   sup_ξ∈Ξ (v(ξ) - λξ - ξ̂_i^2) ≤ s_i    ∀  i ∈ [N];                                                                                                                                              λ≥ 0. ]

For each i ∈ [N], consider the following maximization problem
corresponding to the left-hand side of the constraints in (<ref>):

    [                             h^i(λ) :=                                   sup (Fξ)^Tx - λ(ξ^Tξ - 2ξ̂_i^Tξ + ξ̂_i^2);                                                                                                    Ax = b,  x ≥ 0;                                                                                              x_j ∈{0, 1}   ∀  j ∈;                                                                                                e_1^Tξ = 1,  ξ∈Ξ̂, ]

which is a mixed 0-1 bilinear program. Under Assumption <ref>, 
it holds also that the optimal value of (<ref>) equals the optimal value of
an associated copositive program <cit.>, 
which we now describe. 

Define

    z := ξ x∈^k + n,   
        E := [ -be_1^T       A ]∈^m × (k + n),


    H^i(λ) := [ -λ (I - ξ̂_i e_1^T - e_1ξ̂_i^T + ξ̂_i^2e_1e_1^T)                                          1/2 F^T;                                            1/2 F                                                0 ]∈^k+n,

and for any j ∈, define

    Q_j :=  0  f_j0  f_j^T - 1/20  f_je_1  0^T 
                  -  1/2e_1  00  f_j^T ∈^k+n.
 
where f_j denotes the j-th standard basis vector in ^n. 

Because both  and Ξ are bounded by Assumptions <ref>
and <ref>, there exists a scalar r > 0 such that the constraint
z^T z = ξ^T ξ + x^T x ≤ r is redundant for (<ref>). 
Furthermore, it is well-known that we can use the following
quadratic constraints to represent the binary variables
in the description of : 

    x_j^2 - x_j = 0    ⇔   Q_j∙ zz^T = 0.
 
After adding the redundant constraint and representing the
binary variables, we homogenize problem (<ref>)
as follows:

    [                      max            H^i(λ) ∙ zz^T;                               Ez = 0,  g_1^Tz = 1;                                      I ∙ zz^T ≤ r;                          Q_j ∙ zz^T = 0    ∀  j ∈;                                        z ∈Ξ×_+^n, ]

where g_1 = e_1  0∈^k+n
and e_1 denotes the standard basis 
vector in ^k.  
The copositive representation is thus

    [                   max            H^i(λ) ∙ Z;                                 (EZE^T) = 0;                           g_1 g_1^T ∙ Z = 1;                                   I ∙ Z ≤ r;                       Q_j ∙ Z = 0    ∀  j ∈;                                Z ∈(Ξ×_+^n). ]

Letting u^i ∈^m, ρ^i ∈_+, α^i ∈, and
v^i ∈^|| be the respective dual multipliers of
(EZE^T) = 0, I ∙ Z ≤ r, g_1 g_1^T ∙ Z = 1,
and Q_j ∙ Z = 0, standard conic duality
theory implies the dual of (<ref>) is

    [                                                 min_α^i, ρ^i, u^i, v^i                                                            α^i + r ρ^i;                                                                        α^i g_1 g_1^T - H^i(λ) + E^T (u^i)E + ∑_j ∈ v^i_jQ_j + ρ^i I ∈(Ξ×_+^n);                                                                                                                                      ρ^i ≥ 0. ]

Holding all other dual variables fixed, for ρ^i > 0 large,
the matrix variable in (<ref>) is strictly
copositive—in fact, positive definite—which establishes that
Slater's condition is satisfied, thus ensuring strong duality:
the optimal value of (<ref>) equals
the optimal value of (<ref>). Therefore,
we can reformulate problem (<ref>)
as follows:

    [                                                                      v_(_N, ϵ)^+ =                                                                                min                                                     λϵ^2 + 1/N∑_i=1^N (α^i + rρ^i);                                                                                                                                                                       α^i g_1 g_1^T - H^i(λ) + E^T (u^i)E + ∑_j ∈ v^i_jQ_j + ρ^i I ∈(Ξ×_+^n)  ∀  i ∈ [N];                                                                                                                                                                                                                                   ρ^i ≥ 0     ∀  i ∈ [N];                                                                                                                                                                                                                                                    λ≥ 0. ]



Note that if Assumption <ref> fails, the constraint
I ∙ Z ≤ r should be excluded from (<ref>)
and thus the terms rρ^i and ρ^i I in the objective
function and the constraint, respectively, should be excluded 
in (<ref>) as well. As such, strong duality 
between (<ref>) and (<ref>) cannot be
established in this case. However, the modified (<ref>) still
provides an upper bound on h^i(λ).  Accordingly,
the modified problem (<ref>) still provides an
upper bound on v_(_N, ϵ)^+.



 §.§ A semidefinite-based relaxation
 

As problem (<ref>) is difficult to solve in general, we 
propose a tractable approximation based on semidefinite programming
techniques. In particular, we propose an inner approximation
of (Ξ×_+^n) in
(<ref>) so that the resulting problem has an optimal
value that is an upper bound on v_^+. Now, define

    (Ξ×_+^n) :=
        {
            S + M
             : [ S_11∈(Ξ), Rows(S_21) ∈Ξ^*;           S_22≥ 0,  M ≽ 0 ]},

where (Ξ) is an inner approximation of (Ξ), i.e., 
(Ξ) ⊆(Ξ). 
Immediately, we have a relationship between (Ξ×_+^n)
and (Ξ×_+^n):
 
(Ξ×_+^n) ⊆(Ξ×_+^n).
 


Let arbitrary p  q∈Ξ×_+^n be given. We need
to show

    p  q^T
        (
            S + M
        )
        p  q = p  q^T S p  q +
        p  q^T M p  q≥ 0.


    p  q^T
        (
            S + M
        )
        p  q =     p  q^T S p  q +
        p  q^T M p  q
    
       =     p^T S_11p + 2q^TS_21p + q^TS_22q +  p  q^T M p  q
    ≥    0

The first term is nonnegative because p ∈Ξ and
S_11∈(Ξ) ⊆(Ξ); the second
term is nonnegative because p ∈Ξ, q ≥ 0, and
Rows(S_21) ∈Ξ^*; the third term is
nonnegative because q ≥ 0 and S_22≥ 0; the
last term is nonnegative because M ≽ 0. 

When Ξ= {ξ∈^k: Pξ≥ 0} is a polyhedral cone based
on some matrix P ∈^p × k, a typical inner approximation (Ξ)
of (Ξ) is given by

    (Ξ) := { S_11 = P^T Y P : Y ≥ 0 },

where Y ∈^p is a symmetric matrix variable. This corresponds
to the RLT approach of <cit.>.
When Ξ= {ξ∈^k : (ξ_2, …, ξ_k)^T≤ξ_1 } is
the second-order cone,
it is known <cit.> that

    (Ξ) = { S_11 = τ J + M_11 : τ≥ 0,  M_11≽ 0 },

where J = (1, -1, …, -1). Because of this simple structure,
it often makes sense to take (Ξ) = (Ξ) in practice.

Now consider the following problem by replacing (Ξ×_+^n) with
(Ξ×_+^n) in (<ref>). 

    [                                                                        v̅_(, ϵ)^+=                                                                                min                                                     λϵ^2 + 1/N∑_i=1^N (α^i + rρ^i);                                                                                                                                                                       α^i g_1 g_1^T - H^i(λ) + E^T (u^i)E + ∑_j ∈ v^i_jQ_j + ρ^i I ∈(Ξ×_+^n)  ∀  i ∈ [N];                                                                                                                                                                                                                                   ρ^i ≥ 0     ∀  i ∈ [N];                                                                                                                                                                                                                                                    λ≥ 0. ]


Obviously, we have the following result:

v_(_N, ϵ)^+ ≤v̅_(, ϵ)^+. 






§ NUMERICAL EXPERIMENTS
 

In this section, we validate our proposed Wasserstein-ball
approach (WB) on three applications. We will compare WB
with the moment-based approach (MB) proposed in 
<cit.> where the exact values of the 
first two moments of the distributions are known. 
In practice, the moments of the distribution are often not known 
exactly. To this end, Delage and Ye <cit.> proposed
a data-driven approach to handle this case. However, in this paper, 
we assume that the moments are known exactly for MB. Actually, this 
choice favors MB, but the goal of our experiments is to demonstrate that 
our approach provides a valid upper bound that gets closer to v_ as
the size of the data set increases, while the MB provides an upper bound,
which does not improve with the size of the data set.       


All computations are conducted with Mosek version 8.0.0.28 beta <cit.> on 
an Intel Core i3 2.93 GHz Windows computer with 4GB of RAM and are implemented 
using the modeling language YALMIP <cit.> in MATLAB  (R2014a) 
version 8.3.0.532. In order to demonstrate the effectiveness of WB, we also implement
a Monte Carlo simulation-based approach (SB) which requires a sufficiently large
number of randomly generated samples. For the project management problem in
Section <ref>, a linear program is solved for each sample of the Monte Carlo 
simulation, while for the knapsack problem in Section <ref>, an integer
program is solved for each sample. We employ CPLEX 12.4 to solve these 
linear programs and integer programs.  

 


 §.§ Statistical sensitivity analysis of highest-order statistic
  

The problem of finding the maximum value from a set
ζ = (ζ_1, …, ζ_n) of n numbers can be formulated 
as the optimization problem:

    max{ζ^Tx : e^T x = 1,  x ≥ 0}.

For example, suppose ζ_1 = max{ζ_1, …, ζ_n }, 
then the optimal solution to (<ref>) is x_1^* = 1, x_2^* =
⋯ = x_n^* = 0. For the statistical sensitivity analysis problem,
we consider a random vector ζ following a joint distribution . 
In the situation where the true distribution is not known exactly,
our focus is to investigate the upper bound on the expected
maximum value over an ambiguity set containing distributions
that possess partial shared information. 

We consider an instance with n = 3 and the true distribution
 of ζ is assumed to be jointly lognormal with first and 
second moments given by μ_log∈^3 and 
Σ_log∈^3, respectively. 

In our experiments, we use the following procedure to randomly
generate μ_log and Σ_log. We first sample 
μ∈^3 from a uniform distribution [0, 2]^3. Then, we 
randomly generate a matrix Σ∈^3 as follows:
we set the vector of standard deviations to σ = 14 e ∈^3,
sample a random correlation matrix C ∈^3 using the MATLAB
command `gallery(`randcorr',3)', and set Σ = (σ)C(σ)+μμ^T. 
We set μ and Σ as the first and second moments respectively
of the corresponding normal distribution of . Then μ_log 
and Σ_log can be computed based on the following 
formulae <cit.>:

    [                           (μ_log)_i = e^μ_i + 0.5Σ_ii,; (Σ_log)_ij = e^μ_i + μ_j + 0.5(Σ_ii + Σ_jj)(e^Σ_ij-1). ]



We can cast this problem into our framework by setting
m = 1, k = n+1, ξ = (1, ζ_1, …, ζ_n),
F = (0, I), and = ∅. Obviously, Assumptions
<ref> and <ref> are satisfied. 
Assumption <ref> is vacuous. Although Assumption 
<ref> does
not hold, problem (<ref>) can still provide
a valid upper bound on the expected optimal value as 
discussed in Section <ref>.



  §.§.§ The deviation of empirical Wasserstein radii


In this experiment, we consider a particular underlying
distribution  that is generated by the procedure mentioned
above. Also, we consider eight cases for the size
of the dataset: N ∈{10, 20, 40, 80, 160, 320, 640,
1280 }. For each case, we randomly generate a dataset Θ_N 
containing N independent samples
from  and use the procedure in Section <ref>
to determine a desired radius from a pre-specified set 
ℰ = {0.001, 0.005, 0.01, 0.05, 0.1, 0.2, 
0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 2.0}[
From preliminary experiments, the largest element 2.0 in 
set E returned 1 as the empirical confidence level for all the 
experiments we conducted. Thus, we believe it is sufficient to have
2.0 as the largest element here.]. 
In particular, we set K = 100 in Algorithm <ref>.  
Figure <ref> shows the trend of the reliabilities over different
Wasserstein radii for N ∈{ 20, 80, 320, 1280 }. Clearly, 
smaller Wasserstein radii tend to have lower empirical confidence 
levels. Furthermore, as the sample size increases, the empirical confidence
level increases as well for the same Wasserstein radius. The result of this
experiment indicates that we can practically choose a Wasserstein
radius with a desired statistical guarantee for each case of N. 
We remark that the derived radii can be used for datasets of
the same sizes generated from different distributions of the same family
(lognormal distributions in this application).  




























  §.§.§ Instances with the same underlying distribution
 

Our next experiment is to focus on a particular joint lognormal
distribution . We consider eight cases: N ∈{ 10, 20, 40, 80, 160, 320,
640, 1280 }. For each case, we test 100 trials and in each
trial we randomly generate N independent samples from  
and choose the Wasserstein radius with an empirical confidence level of 0.90.
We compare our approach with MB where the first
two moments are directly given as μ_log and Σ_log. 
We also randomly generate 100000 independent
samples from  to simulate the true expected optimal value.

We demonstrate experimental results in Figure  <ref>. 
Note that the solid black line represents the simulated value of
the true expected optimal value, while the dashed black line 
represents the upper bound calculated by the moment-based 
approach. Furthermore, we solve an instance of (<ref>) 
for each of the 100 trails in each case of N. We use the blue, red,
and green lines to respectively represent the 80^th quantile,
the median, and the 20^th quantile of the values from the 
100 trials in each case. Figure <ref> shows that our approach provides
weaker bounds on the expected optimal value for smaller sample
sizes. However, as the size of samples increases, our approach
provides stronger bounds and the bounds get relatively close to the
simulated value. In contrast, the value from MB remains the same 
regardless of sample size. 





  §.§.§ Instances with different underlying distributions


In this experiment, we consider eight cases N ∈{ 10, 20, 40, 
80, 160, 320, 640, 1280 }. In each case, we randomly 
generate 100 trials. For each trial in each case we generate N
samples from a random lognormal distribution whose first and second 
moments are generated by using the procedure at the beginning of
this section. For each trial in each case, we solve an instance of 
(<ref>) with a Wasserstein radius corresponding to 
an empirical confidence level of 0.90. We also simulate the true expected 
optimal values by randomly generating 100000 samples from the
true distributions.
 
For each trial in each case, we denote the optimal value from (<ref>) 
by v̅_WB^+ and the simulated value by v_SB. 
Then, we calculate the relative gap between WB and SB as 

    gap(WB) := v̅_WB^+ - v_SBv_SB.
 
We take the average of the relative gaps over the 100 trials for each case. 
Then, for each trial in each case, we solve MB with the first two
moments computed by (<ref>). Denote the optimal value
from MB by v̅_MB^+. Similarly, we calculate the relative gap 
between MB  and SB as

    gap(MB) := v̅_MB^+ - v_SBv_SB.
 
We then take the average of the relative gaps over the 100 trials in each case. 
Figure <ref> illustrates the average relative gaps from both WB and
MB over the eight cases. Clearly, the upper bound from WB approaches the simulated
value along with the increase of the size of samples, while the average relative gap 
between the bound from MB and the simulated value does not. 
 

Table <ref> shows 
the percentage of the 100 trials where the optimal values
from WB are greater than or equal to the corresponding simulated optimal values in the eight cases.
The result demonstrates that the derived empirical Wasserstein radii indeed
provide desired statistical guarantees in practice. 












 §.§ Project management problem
 

In this application, we consider a project management problem,
which can be formulated as a longest-path problem on a directed
acyclic graph. The arcs denote activities and the nodes denote the
completions of a set of activities. Arc lengths denote the time
to complete the activities. Thus, the longest path from the starting
node s to the ending node t gives the time needed to compete
the whole project. Let ζ_ij be the length (time) of arc (activity)
from node i to node j. The problem can be solved as a linear
program due to the network flow structure as follows:


    [                                                                                                                                           max                                                                                                                          ∑_(i,j) ∈ Aζ_ij x_ij;                                                                                                                                               ∑_i : (i,j) ∈ A x_ij - ∑_j : (i,j) ∈ A x_ji = 
    {[              1,     if i = s; 0,     if i ∈ N, and i ≠ s,t;             -1,     if i = t ].;                                                                                                                                                                                                                                                                     x_ij≥ 0,  ∀  (i,j) ∈ A, ]

where A denotes the set containing all the arcs, N
denotes the set containing all nodes on the network, and x_ij denotes
the number of units of flow sent from node i to node j through arc 
(i,j) ∈ A. For the stochastic project
management problem, the activity times are random. In such cases,
due to the resource allocation and management constraints, the project manager
would like to quantify the worst-case expected completion time of the project,
which is corresponding to the worst-case longest path of the network. 

We consider an instance with a network structure shown in Figure <ref>.
 
This network consists of 7 arcs and 6 nodes. There
are 3 paths from the starting node to the ending node on the network.
In the experiments of this example, we consider truncated joint normal
distributions. We use the following procedure to generate a truncated joint
normal distribution : denoting | A| by the cardinality of set A, 
we generate ζ≥ 0 from a jointly normal
distribution with first and second moments given by μ∈^| A|
and Σ∈^| A|, respectively. Specifically, we sample
μ from a uniform distribution [0, 5]^| A| while the matrix
Σ is generated randomly using the following procedure: we set the
vector of standard deviations to σ = e, sample a random correlation
matrix C ∈^| A| using the MATLAB command
`gallery(`randcorr',| A|)', and set Σ = (σ)C(σ)+μμ^T. 
Skipping the details, we can cast the network flow problem
into our framework. It is straightforward to check that Assumptions 
<ref>, <ref>, and <ref>
are satisfied and Assumption <ref> is vacuous. 



  §.§.§ Instances with the same underlying distribution


The first experiment of this example focuses on a particular underlying distribution
. We consider seven cases: N ∈{10, 20, 40, 80, 160, 320,
640}. For each case, we run 100 trials and in each trial we randomly
generate a dataset Θ_N containing N independent samples from . 
We use the procedure in Section <ref> to
compute an empirical confidence level set for each case. Then, we use computed
empirical confidence level sets to derive empirical Wasserstein radii for the following computations.  
For each trial in each case, we solve an instance of (<ref>)
with a Wasserstein radius corresponding to an empirical confidence level of 0.90. 
We compare WB with MB where the first two moments are approximated
by using the sample mean and variance from 100000 samples.
The computed moments are close to their theoretical counterparts as the 
sample size is considerably large. We also simulate the expected optimal 
value over the 100000 samples. Figure <ref>
shows that WB provides weaker bounds on the expected optimal value for 
smaller sample sizes. However, as the size of samples increases, WB provides
stronger bounds and the bounds get relatively close to the
simulated value. In contrast, the bounds from MB remains the same 
regardless of the change of sample sizes. 




 



  §.§.§ Instances with different underlying distributions


In this experiment, we consider seven cases: N ∈{ 10, 20, 40, 80, 
160, 320, 640}. For each case, we randomly generate 100 trials in 
which N independent samples are drawn from a randomly generated 
truncated joint normal distribution. Then, for each trial in each case, we solve 
and instance of (<ref>) with a Wasserstein 
radius corresponding to a 0.90 empirical confidence level. 
We solve MB where the first two moments are approximated 
by computing the sample mean and variance of 100000 samples. 
We also simulate the expected optimal value over the 100000 samples 
for each trial in each case. We compute the relative gap between the WB and
SB as well as the relative gap between MB and SB. 
Then, for each case, we take the average of the relative
gaps from both WB and MB over the 100 trials. 
Figure <ref> illustrates the average relative gaps 
over the seven cases. 
Clearly, the upper bound from WB approaches to the simulated
value along with the increase in the size of samples, while the gap between
the bound from MB and the simulated value remains relatively the same
as the sample size increases. Table <ref> shows

the percentage of the 100 trials where the optimal values
from WB are greater than or equal to the corresponding
simulated optimal values over the seven cases.  






 §.§ Knapsack problem
 

A standard knapsack problem is defined as follows: given a set of items
and each with a weight and a value, the problem is to determine the number 
of items to include in a knapsack such that the total weight is less than or equal to
a given capacity limit and the total value is maximized; see the detail in 
<cit.>. Let w_i and ζ_i be the weight and value 
of item i  (i = 1, …, n) respectively. Let W be the maximum weight 
capacity of the knapsack. Then, the knapsack problem can be formulated 
as an integer program:

    v(ζ) := max{∑_i = 1^n ζ_ix_i : ∑_i=1^n w_i x_i ≤ W,  x_i ∈{0, 1}},

where x_i represents the number of item i to include in the knapsack. 
Assume that the values of the items are random and follow an unknown 
joint distribution. Assume also that we can collect a dataset containing N 
samples with each corresponding to an observation of the n item values. 
In such cases, we would like to compute a data-driven distributionally robust 
upper bound on the expected maximum value of the knapsack. We can 
approximate the upper bound by solving problem (<ref>).

We consider an instance with n=4, w = (5,  4,  6,  3 )^T,
and W = 10. The true distribution  of ζ is assumed to be jointly lognormal
with first and second moments given by μ_log∈^4 and 
Σ_log∈^4, respectively. Similar to the procedure described in 
Section <ref>, we sample μ∈^4 from a uniform distribution [0, 2]^4. 
Then, we randomly generate a matrix Σ∈^4 as follows:
we set the vector of standard deviations to σ = 14 e ∈^4,
sample a random correlation matrix C ∈^4 using the MATLAB
command `gallery(`randcorr',4)', and set Σ = (σ)C(σ)+μμ^T. 
Then μ_log and Σ_log can be computed
based on (<ref>). We can easily cast this problem into
our framework. For simplicity, we skip the details. It is also
straightforward to check that the conditions in Assumptions <ref>,
<ref>, and <ref> are satisfied. Although 
Assumption <ref> is not satisfied, we still can solve
(<ref>) to obtain a valid upper bound on the expected optimal
value of the knapsack problem. 



  §.§.§ Instances with the same underlying distribution


In the first experiment, we focus on a particular underlying distribution
 and consider eight cases: N ∈{10, 20, 40, 
80, 160, 320, 640, 1280 }. For each case, we run 100 
trials and in each trial we randomly generate a dataset Θ_N 
containing N independent samples from . Similarly, we derive
empirical Wasserstein radii for each case. Then, for each trial in each
case, we solve an instance of (<ref>) with an empirical Wasserstein radius 
corresponding to an empirical confidence level of 0.90. We compare our approach 
with MB where the first two moments are computed by (<ref>). 
We simulate the expected optimal value over 100000 samples. 
Note that we solve an integer program for each sample in the simulation. 
Figure <ref> shows that WB provides weaker bounds on the 
expected optimal value for smaller sample sizes. However, as the size of 
samples increases, WB provides stronger bounds and the bounds get 
relatively close to the simulated value. In contrast, the bounds from MB 
remains the same regardless of the change of sample sizes. We remark
that the upper bound computed by WB may not converge to the true
expected optimal value as the sample size increases to infinity; see the trend 
shown in Figure <ref>. This is due to the fact that problem 
(<ref>) is a relaxation of problem (<ref>) and the
fact the relaxation is not tight in this example. 




 



  §.§.§ Instances with different underlying distributions


This experiment considers eight cases: N ∈{ 10, 20, 40, 80, 160, 
320, 640, 1280 }. For each case, we randomly generate 
100 trials with each trial is drawn from a randomly generated 
joint lognormal distribution. Then, for each trial in
each case, we solve an instance of (<ref>) with a Wasserstein 
radius corresponding to a 0.90 empirical confidence level. 
We simulate the expected optimal value over 100000 samples 
for each trial in each case. Next, we compute the relative gap 
between the values of WB and SB as well as the relative gap between
the values of MB and SB. Then, for each case, we take the average of 
the relative gaps from both WB and MB over the 100 trials. 
Figure <ref> illustrates the average relative gaps over 
the seven cases. Clearly, the WB gap becomes narrower as the 
sample size increases, while the MB gap remains relatively the same. 

Table <ref> shows

the percentage of the 100 trials where the optimal values
from WB are greater than or equal to the corresponding
simulated optimal values over the eight cases.  







§ CONCLUDING REMARKS
 

In this paper, we have studied the expected optimal value of a mixed 0-1
programming problem with uncertain objective coefficients following a
joint distribution whose information is not known exactly but a set 
of independent samples can be collected. Using the samples, we
have constructed a Wasserstein-based ambiguity set that contains
the true distribution with a desired confidence level. We proposed
an approach to compute the upper bound on the expected optimal
value. Then under mild assumption, the problem was reformulated
to a copositive program, which leads to a semidefinite-based relaxation. 
We have validated the effectiveness of our approach over three
applications.



§ ACKNOWLEDGEMENTS

We would like to thank Kurt Anstreicher, Qihang Lin, Luis Zuluaga,
and Tianbao Yang for many useful suggestions, which helped us
improve the results of the paper. 









plain














