Wasserstein Distributionally Robust Optimization: Theory and Applications in Machine Learning

Daniel Kuhn, Peyman Mohajerin Esfahani, Viet Anh Nguyen, Soroosh Shafieezadeh-Abadeh

Introduction

and the optimal risk is defined as the risk of the least risky admissible loss function, that is,

If the distribution $\mathds{P}$ is unknown, we lack an important input parameter for the risk evaluation problem (1) and the decision problem (2). In this case, the unknown true distribution $\mathds{P}$ could be replaced with a nominal distribution $\widehat{\mathds{P}}_{N}$ estimated from the $N$ training samples. Note that unlike $\mathds{P}$ , the nominal distribution $\widehat{\mathds{P}}_{N}$ is accessible as it is constructed from observable quantities. Therefore, the nominal risk evaluation and decision problems (that is, problems (1) and (2) with $\widehat{\mathds{P}}_{N}$ instead of $\mathds{P}$ ) are at least in principle solvable. The following example showcases common methods for constructing the nominal distribution $\widehat{\mathds{P}}_{N}$ .

In the remainder we will primarily work with the following non-parametric and parametric models for the nominal distribution.

In the absence of any structural information, it is convenient to set $\widehat{\mathds{P}}_{N}$ to the discrete empirical distribution, that is, the uniform distribution on the $N$ training samples,

where $\delta_{\widehat{\xi}_{i}}$ denotes the Dirac point mass at the $i^{\rm th}$ training sample $\widehat{\xi}_{i}$ .

where only the mean vector $\widehat{\mu}$ and the covariance matrix $\widehat{\Sigma}$ depend on the training samples and are constructed via maximum likelihood estimation.

As a function of the training data, the nominal distribution $\widehat{\mathds{P}}_{N}$ constitutes itself a random object, which is governed by the distribution $\mathds{P}^{N}$ of the $N$ independent training samples. $\square$

Even if the most sophisticated statistical tools are deployed, the nominal distribution $\widehat{\mathds{P}}_{N}$ will invariably differ from the unknown true distribution $\mathds{P}$ that generated the training samples. Moreover, if $\widehat{\mathds{P}}_{N}$ is used instead of $\mathds{P}$ , the solutions of the risk evaluation problem (1) and the decision problem (2) are likely to inherit any estimation errors in $\widehat{\mathds{P}}_{N}$ . In the context of financial portfolio theory it has even been observed that estimation errors in the input parameters of an optimization problem are often amplified by the optimization . To make things worse, one can generally show that even if the distributional input parameters of a decision problem are unbiased, the optimization results tend to be optimistically biased. Thus, implementing the optimal decisions leads to disappointment in out-of-sample tests. In decision analysis this phenomenon is sometimes termed the optimizer’s curse , and in stochastic optimization it is referred to as the optimization bias .

One can show that the Wasserstein distance is a metric, that is, it is nonnegative, symmetric and subadditive, and it vanishes only if $\mathds{Q}=\mathds{Q}^{\prime}$ [107, p. 94]. One can further show that $W_{p}(\mathds{Q},\mathds{Q}^{\prime})$ is finite whenever both $\mathds{Q}$ and $\mathds{Q}^{\prime}$ have finite $p^{\rm th}$ -order moments [107, p. 95].

The optimization problem in (5) constitutes an infinite-dimensional linear program over the transportation plan $\pi$ . This linear program admits a strong dual, which in turn provides an alternative characterization of the Wasserstein distance.

For any $p\in[1,\infty)$ , the $p^{\rm th}$ power of the type- $p$ Wasserstein distance between $\mathds{Q}$ and $\mathds{Q}^{\prime}$ admits the dual representation

For a proof of Theorem 1.4 see [107, § 5]. The dual problem can be interpreted as the profit maximization problem of a third party that reallocates the dirt from $\mathds{Q}$ to $\mathds{Q}^{\prime}$ on behalf of the problem owner by buying dirt at the origin $\xi$ at unit price $\phi(\xi)$ and selling dirt at the destination $\xi^{\prime}$ at unit price $\psi(\xi^{\prime})$ . The constraints ensure that the problem owner prefers to use the services of the third party for every origin-destination pair $(\xi,\xi^{\prime})$ instead of reallocating the dirt independently at her own transportation cost $\|\xi-\xi^{\prime}\|^{p}$ . The optimal price functions $\phi^{\star}$ and $\psi^{\star}$ —if they exist—are called Kantorovich potentials [107, p. 99].

The Lipschitz modulus can be viewed as the slope of the steepest line segment connecting any two points on the graph of $\phi$ . The following result simplifies Theorem 1.4 for $p=1$ .

The type-1 Wasserstein distance between $\mathds{Q}$ and $\mathds{Q}^{\prime}$ admits the dual representation

Kantorovich and Rubinstein originally established this result for compactly supported distributions. A modern proof for arbitrary distributions can be found in [107, Remark 6.5]. Theorem 1.5 asserts that the type-1 Wasserstein distance between $\mathds{Q}$ and $\mathds{Q}^{\prime}$ equals the difference between the expected values of a test function $\phi$ under $\mathds{Q}$ and $\mathds{Q}^{\prime}$ , respectively, maximized across all Lipschitz-continuous test functions with Lipschitz modulus of at most 1.

The above reasoning suggests that in order to approximate the (optimal) risk well, one should construct an estimator $\widehat{\mathds{P}}_{N}$ that has a small Wasserstein distance to the unknown true distribution $\mathds{P}$ with high confidence. Unfortunately, however, estimators are subject to fundamental performance limitations and cannot be improved beyond a certain level.

Depending on the available structural information on $\mathds{P}$ , the nominal distributions portrayed in Example 1.1, which will be used throughout this tutorial, are essentially optimal within certain estimator families.

Elliptical distributions: Assume that $\mathds{P}$ is known to be an elliptical distribution with a known density generator $g$ but unknown mean vector $\mu$ and covariance matrix $\Sigma$ . In this case, the problem of finding an estimator $\widehat{\mathds{P}}_{N}$ for the distribution $\mathds{P}$ reduces to finding an estimator $\widehat{\theta}_{N}$ for the vector $\theta=(\mu,\Sigma)$ of unknown distribution parameters. Under mild regularity conditions, the Cramér-Rao inequality guarantees that the covariance matrix of $\sqrt{N}\cdot\widehat{\theta}_{N}$ exceeds the inverse Fisher information matrix in a positive semidefinite sense for any unbiased estimator $\widehat{\theta}_{N}$ . As the maximum likelihood estimator $\widehat{\theta}_{N}^{\,\rm ML}$ is asymptotically unbiased and efficient, i.e., the mean of $\widehat{\theta}_{N}^{\,\rm ML}$ converges to $\theta$ and the variance of $\sqrt{N}\cdot\widehat{\theta}_{N}^{\,\rm ML}$ converges to the inverse Fisher information matrix as $N$ grows, it is asymptotically optimal among all conceivable unbiased estimators.

We emphasize that, by mobilising more powerful results from statistics, the above optimality guarantees could be extended to even larger families of estimators. $\Box$

Using the proposed ambiguity set, we define the worst-case risk as

Problem (7) constitutes a distributionally robust optimization problem. It seeks decisions that have minimum risk under the most adverse distributions in the ambiguity set. Intuitively, problem (7) can thus be viewed as a zero-sum game, where the decision-maker first selects an admissible loss function with the goal to minimize the risk, in response to which some fictitious adversary or ‘nature’ selects a distribution from within the ambiguity set with the goal to maximize the risk. The hope is that by minimizing the worst-case risk, we actually push down the risk under all distributions in the ambiguity set—in particular under the unknown true distribution $\mathds{P}$ , which is contained in the ambiguity set if $\varepsilon$ is large enough. Thus, there is reason to hope that the solutions of distributionally robust optimization problems with carefully calibrated ambiguity sets display low out-of-sample risk.

The distributionally robust risk evaluation and decision problems (6) and (7) are attractive for a multitude of diverse reasons.

Fidelity: Distributionally robust models are more ‘honest’ than their nominal counterparts as they acknowledge the presence of distributional uncertainty. They also benefit from information about the type and magnitude of the estimation errors, which is conveniently encoded in the geometry and size of the ambiguity set.

Managing expectations: Due to the optimizer’s curse, the solutions of nominal decision problems equipped with noisy estimators display an optimistic in-sample risk, which cannot be realized out of sample; see Example 1.2. In contrast, the solutions of distributionally robust decision problems are guaranteed to display an out-of-sample risk that falls below the worst-case optimal risk whenever the ambiguity set contains the unknown true distribution. Thus, nominal decision problems over-promise and under-deliver, while distributionally robust decision problems under-promise and over-deliver.

Computational tractability: The distributionally robust problems (6) and (7) can often be reformulated as (or tightly approximated by) finite convex programs that are solvable in polynomial time. Section 2 will showcase some key tractability results.

Performance guarantees: For judiciously calibrated ambiguity sets, one can prove that the worst-case optimal risk for any fixed sample size $N$ provides an upper confidence bound on the out-of-sample risk attained by the optimizers of (7) (finite sample guarantee) and that the optimizers of (7) converge almost surely to an optimizer of (2) as $N$ tends to infinity (asymptotic guarantee); see Section 3.

Regularization by robustification: The optimizer’s curse is reminiscent of overfitting phenomena that plague most statistical learning models. One can show that distributionally robust learning models equipped with a Wasserstein ambiguity set are often equivalent to regularized learning models that minimize the sum of a nominal objective and a norm term that penalizes hypothesis complexity. Similarly, one can show that some distributionally robust maximum likelihood estimation models produce shrinkage estimators. Thus, Wasserstein distributional robustness offers new probabilistic interpretations for popular regularization techniques. The empirical success of regularization methods in statistics fuels hope that Wasserstein distributionally robust models can effectively combat the optimizer’s curse across many application areas. Connections between robustification and regularization will be explored in Section 4.

Anticipating black swans: If uncertainty is modeled by the empirical distribution, then the nominal decision problem evaluates the admissible loss functions only at the training samples. However, possible future uncertainty realizations that differ from all training samples but could have devastating consequences (‘black swans’) are ignored. If the empirical distribution may be perturbed within a Wasserstein ball with a positive radius, on the other hand, then (possibly small amounts of) probability mass can be moved anywhere in the support set $\Xi$ . Thus, the Wasserstein distributionally robust decision problem faithfully anticipates the possibility of black swans. We emphasize that all distributions in a Kullback-Leibler divergence ball must be absolutely continuous with respect to the nominal distribution, which implies that the corresponding distributionally robust decision problems ignore the possibility of black swans.

Optimality principle: Data-driven optimization aims to use the training data directly to construct an estimator for the objective of problem (2) (a predictor) and a decision that minimizes this predictor (a prescriptor) without the detour of constructing an estimator for $\mathds{P}$ . It has been shown that optimal predictors and the corresponding prescriptors can be constructed by solving a meta-optimization model that minimizes the in-sample risk of the predictor-prescriptor pairs subject to constraints guaranteeing that the in-sample risk is actually attainable out of sample. It has been shown that this meta-optimization problem admits a unique solution: the best predictor-prescriptor pair is obtained by solving a distributionally robust optimization problem over all distributions in some neighborhood of the empirical distribution [77, Theorem 7]. Thus, if one aims to transform training data to decisions, it is in some precise sense optimal to do this by solving a data-driven distributionally robust optimization problem.

Distributionally robust optimization models with Wasserstein ambiguity sets were introduced in . Reformulations of these models as nonconvex optimization problems as well as initial attempts to solve these problems via algorithms from global optimization are reported in and [83, § 7.1]. In the next section we will review convex reformulations and approximations that were discovered in and significantly generalized in .

Computation

The aim of this section is to show that the worst-case risk evaluation problem (6) and the distributionally robust decision problem (7) are computationally tractable in many situations of practical interest. Note first that checking whether a fixed distribution $\mathds{Q}$ is feasible in (6) requires computing the Wasserstein distance $W_{p}(\mathds{Q},\widehat{\mathds{P}}_{N})$ . It is therefore instructive to study the complexity of evaluating Wasserstein distances between arbitrary distributions.

Computing the Wasserstein distance between two discrete distributions amounts to solving a tractable linear program that is susceptible to the network simplex algorithm as well as dual ascent methods or specialized auction algorithms , etc. The set of feasible transportation plans is termed the transportation polytope and displays many useful theoretical properties, which are surveyed in . The need to evaluate Wasserstein distances between increasingly fine-grained histograms has recently motivated efficient approximation schemes. When augmented with an entropic regularization term, for instance, the finite-dimensional transportation problem can be solved quickly by using Sinkhorn’s algorithm . Variants of this approach use Tikhonov regularizers , Bregman divergences or Tsallis entropies instead of the entropic regularization term. A survey of algorithms for the finite-dimensional transportation problem is provided in .

As soon as at least one of the two involved distributions ceases to be discrete, the Wasserstein distance can no longer be evaluated in polynomial time. Even in the simplest imaginable scenario where one distribution is uniform on a hypercube and the other distribution is discrete with two atoms, computing the Wasserstein distance becomes intractable .

Computing the type- $p$ Wasserstein distance between two distributions $\mathds{Q}$ and $\mathds{Q}^{\prime}$ is #P-hard even if $\|\cdot\|$ is the Euclidean norm, $\mathds{Q}$ is the uniform distribution on the standard hypercube $^{m}$ , and $\mathds{Q}^{\prime}$ is a discrete distribution supported on only two points.

If $p=2$ and $\|\cdot\|$ is the Euclidean norm, then the Wasserstein distance admits an analytical lower bound that depends only on the distributions’ first- and second-order moments. This bound is available for any pair of distributions even if their exact Wasserstein distance cannot be computed efficiently. Moreover, the bound is exact for elliptical distributions.

The bound is exact if $\mathds{Q}$ and $\mathds{Q}^{\prime}$ are elliptical distributions with the same density generator.

The inequality (8) may be loose if $\mathds{Q}$ and $\mathds{Q}^{\prime}$ are elliptical distributions with different density generators. Maybe unexpectedly, however, the Wasserstein distance between two elliptical distributions with the same density generator $g$ is actually independent of $g$ . In its general form, Theorem 2.2 is due to Gelbrich . The exact formula for the type-2 Wasserstein distance between normal distributions has been discovered earlier in .

As any Wasserstein ball with a strictly positive radius contains non-discrete distributions (the nominal distribution can be smeared out even if the transportation budget is small), it is perhaps surprising that the worst-case risk evaluation problem (6) may be tractable at all. Indeed, Theorem 2.1 indicates that checking feasibility is already hard in general. We will see below, however, that the extremal distributions determining the worst-case risk are often structurally equivalent to the nominal distribution. Thus, there is hope that problems (6) and (7) become tractable if we choose a nominal distribution with a particularly simple structure (e.g., a discrete or an elliptical distribution).

In the remainder of this section, we will first review tractable bounds on the worst-case risk and present a strong duality result that paves the way towards exact tractable reformulations (Section 2.1). Next, we will delineate efficient methods to compute the worst-case risk as well as the underlying worst-case distributions in situations when the nominal distribution is discrete (Section 2.2) or elliptical (Section 2.3).

Before attempting to derive exact tractable reformulations for the worst-case risk (6), we focus on the simpler task of establishing efficiently computable upper and lower bounds. To derive a pessimistic upper bound, we note that the transportation cost $\|\xi-\xi^{\prime}\|^{p}$ is a convex function of the random variable $\|\xi-\xi^{\prime}\|$ for any $p\geq 1$ . Jensen’s inequality thus implies

where the equality follows from the definition of the worst-case risk, while the second inequality is a direct consequence of the Kantorovich-Rubinstein theorem (see Theorem 1.5). We summarize the above reasoning in the following theorem.

An optimistic lower bound on the worst-case risk can be obtained by replacing the Wasserstein ball in (6) with a smaller ambiguity set. If the distributions in the restricted Wasserstein ball admit a finite parameterization, then the lower bounding problem coincides with a finite optimization problem. Depending on the parameterization, this problem may even be convex. If $\widehat{\mathds{P}}_{N}$ is the empirical distribution, for example, one may restrict the original Wasserstein ball to a subset that contains only perturbed empirical distributions of the form

In the next sections we will describe specific settings in which (6) and (10) are tractable.

2 Tractability Results for Empirical Nominal Distributions

Assume now that the Wasserstein ambiguity set is centered at the empirical distribution defined in (3). In this case, under a mild convexity assumption, the worst-case risk (6) can be exactly expressed as the optimal value of a finite convex optimization problem.

In the limit when $p$ tends to 1 and $q$ to $\infty$ , the function $\varphi(q)$ decays as $1/q$ , while $\|u_{ij}/\gamma\|_{*}^{q}$ grows exponentially whenever $\|u_{ij}\|_{*}>\gamma$ . Thus, we have

For $p=1$ , the constraints of the finite convex program (2.6) are thus equivalent to

In the opposite limit when $p$ tends to $\infty$ and $q$ to 1, the function $\varphi(q)$ converges to 1, and therefore it is easy to see that the constraints of problem (2.6) simplify to

Thus, the variable $\gamma$ disappears from the constraints. As the objective function coefficient of $\gamma$ is nonnegative, this implies that $\gamma=0$ at optimality. $\Box$

Under the conditions of Theorem 2.6 we have

For a proof of Theorem 2.8 we refer to [60, Theorem 4.4] and . We emphasize that problem (15) is always solvable because it has a compact feasible set and an upper semicontinuous objective function, and thus the use of the maximization operator is justified.

For $p=1$ , the last constraint of (15) simplifies to

To analyze the limit when $p$ tends to $\infty$ , we divide the last constraint of (15) by $\varepsilon^{p}$ and observe that $\|\theta_{ij}/(\varepsilon\,\alpha_{ij})\|^{p}$ grows exponentially with $p$ if $\|\theta_{ij}/\alpha_{ij}\|>\varepsilon$ . Otherwise, $\|\theta_{ij}/(\varepsilon\,\alpha_{ij})\|^{p}$ remains bounded by 1 for all $p$ . For $p=\infty$ , the last constraint of (15) is therefore equivalent to the requirement that $\|\theta_{ij}\|\leq\varepsilon\,\alpha_{ij}$ for all $i\in[N]$ and $j\in[J]$ .

An intimate connection between distributionally robust optimization with type- $\infty$ Wasserstein balls and classical robust optimization has first been discovered in . $\Box$

Even though problem (15) is guaranteed to have an optimal solution, the worst-case risk (6) may not be attained by any distribution if $p=1$ . An instance of problem (6) that fails to be solvable is constructed in Example 2.10 below, which replicates [60, Example 2].

If $\nu_{\infty}=\emptyset$ , one can show that

is an extremal distribution that solves (6). For $p>1$ , the last constraint in (15) ensures that $\theta_{ij}=0$ whenever $\alpha_{ij}=0$ because otherwise $\alpha_{ij}\|\theta_{ij}/\alpha_{ij}\|^{p}$ evaluates to $\infty$ . This implies that the set $\nu_{\infty}$ is empty. Thus, for $p>1$ , the worst-case risk (6) of a piecewise concave loss function is always attained by the discrete distribution $\mathds{Q}^{\star}$ constructed above.

If $\nu_{\infty}\neq\emptyset$ , which is only possible in the special case $p=1$ , the distributions

are feasible and asymptotically optimal in (6) as $n\geq|\nu_{\infty}|$ tends to infinity. Intuitively, these distributions send some atoms with decaying probabilities to infinity along specific recession directions $\theta_{ij}^{\star}$ , $(i,j)\in\nu_{\infty}$ , of the support set. Note that moving an atom to infinity is possible even when only a finite (type-1) transportation budget is available provided that the probability mass transported is inversely proportional to the transportation distance.

For $p>1$ , atoms can also migrate to infinity at a finite transportation cost provided that their probabilities are inversely proportional to the $p^{\rm th}$ power of the transportation distance. As piecewise concave loss functions grow at most linearly, however, the decay in probability always outweighs the increase in loss. This reasoning provides an intuitive explanation for our insight that $\nu_{\infty}=\emptyset$ and that the supremum in (6) is always attained for $p>1$ .

One can show that the supremum of the worst-case risk evaluation problem (6) is never attained under the conditions of Theorem 2.11, that is, any asymptotically optimal sequence of distributions must push some (decreasing amount of) probability mass to infinity. As in the case of a piecewise concave loss function, such a sequence can be constructed explicitly. To do so, choose a maximizer $z^{\star}$ of problem (16), which is generally intractable as pointed out in Remark 2.12. Moreover, select $i_{0}\in[N]$ and $\xi^{\star}\in\arg\max_{\|\xi\|\leq 1}\xi^{\top}z^{\star}$ . Then, the distributions

can be shown to be feasible and asymptotically optimal in (6) as $n\geq 1$ tends to infinity.

Note that substituting the SDP (17) into the distributionally robust decision problem (7) yields a tractable SDP if the set $\mathcal{L}$ of admissible loss functions is defined through SDP constraints in $Q$ and $q$ . In order to construct an extremal distribution that solves problem (6) for a fixed convex quadratic loss function, it is useful to derive the dual of the SDP (17).

Suppose that all conditions of Theorem 2.13 hold. If $\lambda_{\rm max}(Q)$ denotes the largest eigenvalue of $Q$ , then

Problem (18) represents a quadratically constrained quadratic program (QCQP) with a compact feasible set and is therefore solvable. As $Q$ is not necessarily negative semidefinite, problem (18) is generally nonconvex. This is perhaps puzzling because (18) is obtained by ‘massaging’ the dual of (17) and because dual optimization problems are convex by construction. The apparent contradiction is resolved by noting that nonconvex QCQPs of the form (18) with a single constraint are equivalent to convex SDPs by virtue of the celebrated $\mathcal{S}$ -procedure [14, Appendix B.1].

are feasible and asymptotically optimal in (6) as $n\geq 1$ tends to infinity.

If the worst-case risk over a Wasserstein ball centered at the empirical distribution is attained, then there always exists an extremal distribution with $N+1$ atoms that can be characterized in quasi-closed form [38, Corollary 2]. In practice, however, it is often convenient to ignore this minimal representability and to search over candidate distributions with more than $N+1$ atoms, e.g., by solving a finite convex optimization problem such as (2.8). For generic nominal distributions, necessary and sufficient conditions for the existence of an extremal distribution are detailed in [38, Corollary 1].

3 Tractability Results for Elliptical Nominal Distributions

We first define an uncertainty set in the space of mean vectors and covariance matrices.

The uncertainty set $\mathcal{U}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ can conveniently be used in classical robust optimization. Indeed, a robust constraint that requires a concave function $h(\mu,\Sigma)$ to be nonpositive for all $(\mu,\Sigma)\in\mathcal{U}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ can be reformulated as a convex constraint that involves the conjugate of $-h(\mu,\Sigma)$ and the support function of the uncertainty set $\mathcal{U}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ [3, Theorem 2], that is,

This constraint is computationally tractable for many commonly used constraint functions because the support function of $\mathcal{U}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ is SDP-representable .

Unlike the mean vector $\mu=\mathds{E}^{\mathds{Q}}[\xi]$ and the second-order moment matrix $M=\mathds{E}^{\mathds{Q}}[\xi\xi^{\top}]$ , both of which constitute linear functions of the underlying distribution $\mathds{Q}$ , the covariance matrix $\Sigma=M-\mu\mu^{\top}$ is nonlinear in $\mathds{Q}$ . The condition $(\mu,\Sigma)\in\mathcal{U}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ thus appears to be nonconvex in $\mathds{Q}$ . To gain a clearer understanding, it is instructive to introduce the uncertainty set $\mathcal{V}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ for $(\mu,M)$ induced by the uncertainty set $\mathcal{U}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ for $(\mu,\Sigma)$ , that is,

Maybe surprisingly, even though it is defined as the pre-image of a convex set under a nonlinear transformation, one can prove that $\mathcal{V}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ is convex. This implies, counterintuitively, that the condition $(\mu,\Sigma)\in\mathcal{U}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ is actually convex in $\mathds{Q}$ because it is equivalent to the requirement $(\mu,M)\in\mathcal{V}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ and because the moments $(\mu,M)$ are linear in $\mathds{Q}$ .

Thanks to its convexity, the uncertainty set $\mathcal{V}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ can again conveniently be used in classical robust optimization. Indeed, a robust constraint that requires a concave function $h(\mu,M)$ to be nonpositive for all $(\mu,M)\in\mathcal{V}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ can be reformulated as a simple convex constraint involving the conjugate of $-h(\mu,M)$ and the support function of $\mathcal{V}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ . This constraint is computationally tractable for many commonly used constraint functions because the support function of $\mathcal{V}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ is SDP-representable .

A useful ambiguity set in the space of probability distributions is the Gelbrich hull, which is constructed as the pre-image of $\mathcal{U}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ under the mean-covariance projection.

Thus, the Gelbrich hull can be expressed as the pre-image of the convex set $\mathcal{V}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ under a linear transformation, which shows that is is actually convex. We emphasize that convexity is not apparent from Definition 2.19, which introduces the Gelbrich hull as the pre-image of a convex set under a nonlinear transformation.

If we define $\mathcal{P}(\Xi,\mu,\Sigma)$ as the Chebyshev ambiguity set that contains all distributions on $\Xi$ with mean vector $\mu$ and covariance matrix $\Sigma$ , then the Gelbrich hull can also be expressed as

Theorem 2.20 immediately implies that the (optimal) Gelbrich risk provides an upper bound on the (optimal) worst-case risk whenever $p\geq 2$ .

Note that (26b) follows immediately from the definition of the uncertainty set $\mathcal{V}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ and the formula for the covariance matrix in terms of the mean vector and the second-order moment matrix. The inner problems in (26a) and (26b) both represent the same distributionally robust optimization problem over a Chebyshev ambiguity set but with different parameterizations. This problem can be viewed as an infinite-dimensional linear program over all probability distributions $\mathds{Q}$ that satisfy the linear equality constraints $\mathds{E}^{\mathds{Q}}[\xi]=\mu$ and $\mathds{E}^{\mathds{Q}}[\xi\xi^{\top}]=M$ . Therefore, the optimal value of the inner maximization problem is concave in the right hand side parameters $\mu$ and $M$ but generally nonconcave in the alternative parameters $\mu$ and $\Sigma$ . The outer problem in (26a) hedges against ambiguity in the mean vector and the covariance matrix, while the one in (26b) hedges against ambiguity in the first- and second-order moments. The formulation (26a) is conceptually appealing because of its connection to the Wasserstein distance and because it is more natural to characterize a distribution in terms of its mean vector and covariance matrix. The formulation (26b), on the other hand, is computationally attractive because it expresses the outer problem as a convex program that maximizes a manifestly concave function over the convex set $\mathcal{V}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ .

In order to construct an extremal distribution for the Gelbrich risk evaluation problem (24), it is again expedient to derive the dual of the SDP (27).

If all conditions of Theorem 2.23 hold, then we have

Note that problem (28) has a continuous objective function as well as a compact feasible set and is therefore solvable. Any optimal solution $(\mu^{\star},\Sigma^{\star},\{\alpha_{j}^{\star},\theta_{j}^{\star},\Theta_{j}^{\star}\}_{j})$ can in principle be used to construct an extremal distribution $\mathds{Q}^{\star}$ that attains the supremum in the Gelbrich risk evaluation problem (24). Specifically, for any $j\in[J]$ let $\mathds{Q}_{j}^{\star}$ be any distribution supported on

If $\alpha^{\star}_{j}>0$ , we impose the additional requirement that $\mathds{Q}^{\star}_{j}$ has mean value $\theta_{j}^{\star}/\alpha_{j}^{\star}$ and second-order moment matrix $\Theta_{j}^{\star}/\alpha_{j}^{\star}$ . Such a distribution is indeed guaranteed to exist. One can then show that the mixture distribution $\mathds{Q}^{\star}=\sum_{j\in[J]}\alpha^{\star}_{j}\cdot\mathds{Q}^{\star}_{j}$ is optimal in (24). By construction, this distribution $\mathds{Q}^{\star}$ has mean vector $\mu^{\star}$ and covariance matrix $\Sigma^{\star}$ . We emphasize that problem (28) can be reformulated as a tractable SDP by applying the variable substitution $M\leftarrow\Sigma+\mu\mu^{\top}$ , replacing the constraint $(\mu,\Sigma)\in\mathcal{U}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ with $(\mu,M)\in\mathcal{V}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ and recalling from Lemma 2.18 that the uncertainty set $\mathcal{V}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ is SDP-representable. Thus, problem (28) can be solved in polynomial time. Even though the mixture components $\mathds{Q}^{\star}_{j}$ , $j\in[J]$ , are guaranteed to exist, however, one can prove that it is NP-hard to construct them. In other words, even though it is easy to solve (28) and even though any solution of (28) gives rise to a solution $\mathds{Q}^{\star}$ of the Gelbrich risk evaluation problem (24), constructing $\mathds{Q}^{\star}$ remains hard.

While exactly computable in polynomial time, the Gelbrich risk of a piecewise quadratic loss function may only provide a loose upper bound on the worst-case risk under the Wasserstein ambiguity set, which is often the actual quantity of interest. One can prove, however, that the Gelbrich risk (24) coincides with the worst-case risk (6) with respect to a type-2 Wasserstein ball if the loss function is quadratic and the nominal distribution is elliptical.

where the equality holds because $\mathds{Q}^{\star}$ is optimal in the Gelbrich risk evaluation problem (24), while the two inequalities follow from the feasibility of $\mathds{Q}^{\star}$ in the worst-case risk evaluation problem (6) and Corollary 2.21, respectively. Thus, all inequalities in the above expression are exact, which implies that $\mathds{Q}^{\star}$ is actually optimal in (6).

Next, we show how $\mathds{Q}^{\star}$ can be constructed from the optimality conditions of the SDP (29).

If all conditions of Theorem 2.25 hold, $\widehat{\Sigma}\succ 0$ and there exists $\gamma^{\star}\geq 0$ with $\gamma^{\star}I\succ Q$ that solves the nonlinear algebraic equation

then the Gelbrich risk (24) is attained by any distribution with mean vector

Moreover, if $\widehat{\mathds{P}}_{N}=\mathcal{E}_{g}(\widehat{\mu},\widehat{\Sigma})$ is elliptical, $p=2$ and $\|\cdot\|=\|\cdot\|_{2}$ is the Euclidean norm, then the elliptical distribution $\mathds{Q}^{\star}=\mathcal{E}_{g}(\mu^{\star},\Sigma^{\star})$ attains the worst-case risk in (6).

One can show that if $Q\succeq 0$ , then $\gamma^{\star}$ exists and $\Sigma^{\star}\succeq\lambda_{\min}(\widehat{\Sigma})I$ . To give an intuition for Theorem 2.26, note that the SDP (29) can be converted to an equivalent nonlinear program (NLP) in the single decision variable $\gamma$ by using Schur complements to show that

at optimality. The resulting NLP minimizes a strictly convex objective function that explodes as $\gamma$ drops to $\lambda_{\max}(Q)$ or as $\gamma$ tends to infinity. Equation (30) represents its first-order optimality condition, whose unique solution $\gamma^{\star}$ can be computed efficiently to any precision via bisection or the Newton-Raphson method. Using (30), one can then show that any distribution with mean vector $\mu^{\star}$ and covariance matrix $\Sigma^{\star}$ as defined in (31a) and (31b), respectively, is indeed feasible and optimal in (24).

It is instructive to contrast Theorem 2.25 with Theorem 2.13, both of which provide exact tractable SDP reformulations for the problem of evaluating the worst-case risk of a quadratic loss function with respect to a type-2 Wasserstein ball. We highlight that the SDP (29) derived in Theorem 2.25 for elliptical nominal distributions accommodates only two linear matrix inequalities, while the SDP (17) derived in Theorem 2.13 for empirical nominal distributions involves $N$ linear matrix inequalities and may thus be considerably harder to solve.

Performance Guarantees

We now argue that for judiciously calibrated Wasserstein ambiguity sets, the worst-case risk (6) associated with a finite sample size $N$ provides an upper confidence bound on the true risk (1) for all admissible loss functions (finite sample guarantee) and that the worst-case optimal risk (7) converges almost surely to the true optimal risk (2) as $N$ tends to infinity (asymptotic guarantee). Intuitively, the finite sample guarantee ensures that the out-of-sample risk will fall short of the worst-case risk with high confidence when we implement an optimizer of the distributionally robust decision probelm (7), while the asymptotic guarantee formalizes the simple intuition that more data enables us to make better decisions.

Concentration inequalities for the nominal distribution $\widehat{\mathds{P}}_{N}$ and its moments can be used to derive finite sample and asymptotic guarantees. If $\widehat{\mathds{P}}_{N}$ is the empirical distribution, for instance, one can prove that $\widehat{\mathds{P}}_{N}$ converges exponentially fast to the data-generating distribution $\mathds{P}$ , in probability with respect to the Wasserstein distance, as $N$ tends to infinity.

The concentration inequality portrayed in Theorem 3.1 gives rise to the following finite sample guarantees [60, Theorem 3.5].

Assume that all conditions of Theorem 3.1 hold and $\varepsilon_{p,N}(\eta)$ is defined as in (3.1). Then, for all $\eta\in(0,1)$ and $\varepsilon\geq\varepsilon_{N}(\eta)$ we have

Requiring the Wasserstein ball to cover $\mathds{P}$ with high confidence is only a sufficient but not a necessary condition for the finite sample guarantees (32a) and (32b). Indeed, these guarantees can be sustained even if the Wasserstein radius is reduced below $\varepsilon_{p,N}(\eta)$ , which is essentially the smallest radius for which the Wasserstein ball represents a $(1-\eta)$ -confidence set for $\mathds{P}$ . The minimal Wasserstein radius that preserves the finite sample guarantees (32a) and (32b) often decays significantly faster than $\mathcal{O}(N^{-\frac{p}{m}})$ without suffering from a curse of dimensionality. If $p=1$ , the data-generating distribution is absolutely continuous with respect to the Lebesgue measure and the set $\mathcal{L}$ of admissible loss functions admits a smooth parameterization, for example, one can show that a Wasserstein radius of the order $\mathcal{O}(\sqrt{\log m/N})$ maintains finite sample guarantees akin to (32a) and (32b), which is consistent with recent findings in the compressed sensing and high-dimensional statistics literature [10, Theorem 1]. $\Box$

As the number $N$ of training samples grows, one can simultaneously reduce the Wasserstein radius $\varepsilon$ and the significance level $\eta$ without sacrificing the finite sample guarantees (32a) and (32b), which allows us to prove asymptotic consistency [60, Theorem 3.6].

Next, we describe a concentration inequality for the sample mean and the sample covariance matrix that has ramifications for the Gelbrich risk minimization problem (25).

Suppose that the unknown true distribution $\mathds{P}$ has mean vector $\mu$ and covariance matrix $\Sigma$ and that there are $\alpha>2$ and $A>0$ such that $\mathds{E}^{\mathds{P}}[\exp(\|\xi\|_{2}^{\alpha})]\leq A$ . Then, there is $c>1$ that depends on $\mathds{P}$ only through $\mu$ , $\Sigma$ , $\alpha$ , $A$ , and $m$ such that for any $\eta\in(0,1]$ the sample mean $\widehat{\mu}$ and the sample covariance matrix $\widehat{\Sigma}$ satisfy the concentration inequality $\mathds{P}^{N}[(\mu,\Sigma)\in\mathcal{U}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})]\geq 1-\eta$ whenever $\varepsilon$ exceeds \be ε_N(η) = log(c /η)N . \ee

Theorem 3.5 asserts that the uncertainty set $\mathcal{U}_{\varepsilon}(\widehat{\mu},\widehat{\Sigma})$ with radius $\varepsilon\geq\varepsilon_{N}(\eta)$ represents a $(1-\eta)$ -confidence set for the mean vector and covariance matrix of the unknown data-generating distribution $\mathds{P}$ . The critical radius $\varepsilon_{N}(\eta)$ of this confidence set decays as $\mathcal{O}(N^{-\frac{1}{2}})$ and is therefore—unlike the critical radius (3.1)—not subject to a curse of dimensionality.

Theorem 3.5 strengthens [85, Theorem 2.3], which leverages a generalized central limit theorem to show that the type-2 Wasserstein distance between two normal distributions with true and empirical moments, respectively, decays asymptotically as $\mathcal{O}(N^{-\frac{1}{2}})$ . A generalization of this result to elliptical distributions is discussed in [85, Remark 2.4].

Assume that all conditions of Theorem 3.5 hold and $\varepsilon_{N}(\eta)$ is defined as in (3.5). Then, for all $\eta\in(0,1)$ and $\varepsilon\geq\varepsilon_{N}(\eta)$ we have

Theorem 3.6 asserts that the Gelbrich risk (24) offers an upper confidence bound on the true risk (1) under the unknown data-generating distribution uniformly across all loss functions. It also asserts that the optimal value of the Gelbrich risk optimization problem (25) provides an upper confidence bound on the out-of-sample performance of its optimizers.

Distributionally Robust Optimization in Machine Learning

We now demonstrate that the theory of data-driven distributionally robust optimization with Wasserstein ambiguity sets has interesting ramifications for statistical learning and motivates new approaches for addressing fundamental learning tasks such as classification (Section 4.1), regression (Section 4.2), maximum likelihood estimation (Section 4.3) or minimum mean square error estimation (Section 4.4). We conclude with an overview of other applications of distributionally robust optimization in machine learning (Section 4.5).

The distributionally robust classification problem (34) encapsulates two interesting special cases. First, if the Wasserstein radius is set to $\varepsilon=0$ , then (34) collapses to the standard empirical risk minimization problem that minimizes the average prediction error across the training samples. Moreover, if the parameter $\kappa$ appearing in the definition of the norm tends to infinity, then (34) reduces to a classical regularized empirical risk minimization problem.

Recall that the norm $\|\xi\|=\|x\|+\frac{\kappa}{2}|y|$ on the input-output space encodes the transportation cost in the definition of the Wasserstein distance. Thus, $\kappa$ can be viewed as the cost of switching an output from $+1$ to $-1$ or vice versa. If $\kappa=\infty$ , then all distributions in the Wasserstein ball are obtained by perturbing the empirical distribution along the input space because perturbations along the output space would be infinitely expensive. By setting $\kappa=\infty$ , one thus postulates that there is only input uncertainty but no output uncertainty.

Proposition 4.1 gives commonly used regularization techniques a robustness interpretation, which applies under the premise that there is no output uncertainty. It identifies the regularization weight with the Wasserstein radius $\varepsilon$ and the regularization function with the dual of the norm that determines the transportation cost along the input space.

2 Distributionally Robust Regression

Assume now that the norm on the input-output space satisfies $\|\xi\|=\|x\|+\frac{\kappa}{2}|y|$ , where $\|x\|$ is an arbitrary norm on the input space, while $\kappa>0$ quantifies the relative importance of outputs versus inputs. In the absence of output uncertainty (that is, for $\kappa=\infty$ ), there is again an intimate relation between robustification and regularization.

Moreover, if $L(z)$ is the square error and $p=2$ , then problem (35) reduces to

Proposition 4.2 asserts that if there is no output uncertainty, then the distributionally robust regression problem (35) reduces to a regularized empirical risk minimization problem, where the regularization function is given by the dual of the norm on the input space. For Lipschitz continuous univariate loss functions $L(z)$ and for $p=1$ , one simply minimizes the sum of the empirical risk and the regularization term weighted by the product of Wasserstein radius and the Lipschitz modulus of $L(z)$ . Note that the Huber loss, the $\delta$ -insensitive loss and the pinball loss are all Lipschitz continuous with Lipschitz moduli $\delta$ , 1 and $\max\{\delta,1-\delta\}$ , respectively. For the squared loss we need to set $p=2$ because the type-2 Wasserstein ball is the largest Wasserstein ball for which the worst-case expected loss is finite. In this case, one minimizes a combination of the square root of the empirical loss and the regularization term. If one measures distances in the input space using the $\infty$ -norm, then this convex program reduces to the so-called generalized LASSO (Least Absolute Shrinkage and Selection Operator) estimation problem. For further details on distributionally robust regression see .

3 Distributionally Robust Maximum Likelihood Estimation

While $\Sigma$ serves as an input for many problems in engineering, science or economics, it is often the precision matrix $\Sigma^{-1}$ that appears in their solutions. For example, in mean-variance portfolio analysis the portfolio variance to be minimized depends on the covariance matrix of the asset returns, while the optimal portfolio weights depend on the precision matrix. Similarly, linear discriminant analysis uses the covariance matrix of the features as an input and outputs a maximum likelihood classifier that depends on the precision matrix. Moreover, the optimal fingerprint method for climate change detection requires the covariance matrix of the internal climate variability as an input and outputs a climate change signal depending on the precision matrix. Thus, it is often more important to know the precision matrix than the covariance matrix. To ensure that the precision matrix is well defined, we will henceforth assume that $\Sigma\succ 0$ . Unfortunately, the sample covariance matrix is rank-deficient in the big-data regime when the dimension of $\xi$ exceeds the sample size ( $m>N$ ) even if $\Sigma$ has full rank. In this case, one cannot invert $\widehat{\Sigma}$ to obtain a meaningful precision matrix estimator.

From now on we will assume that the unknown true distribution $\mathds{P}$ of $\xi$ is normal. Thus, the problem of maximizing the log-likelihood of the training samples reduces to the following convex program over all candidate mean vectors $\mu$ and precision matrices $X$ [14, § 7.1]. \be inf_μ∈R^m, X ∈S_+^m { -logdetX + 1N ∑_i=1^N (^ξ_i - μ)^⊤X (^ξ_i - μ) } \eeUnfortunately, this maximum likelihood estimation (MLE) problem is unbounded for $N\leq m$ and (almost surely) solved by $\mu^{\star}=\widehat{\mu}$ and $X^{\star}=\widehat{\Sigma}^{-1}$ for $N>m$ . Thus, we fail again to find an estimator in the big-data regime and simply recover the sample mean and the sample covariance matrix in the small-data regime. To overcome this deficiency, we robustify the MLE problem against all distributions within a type-2 Wasserstein ball centered at the normal nominal distribution $\widehat{\mathds{P}}_{N}=\mathcal{N}(\widehat{\mu},\widehat{\Sigma})$ , that is, we solve the robust MLE problem \be inf_μ∈R^m, X∈S_+^m {-logdetX + sup_Q ∈B_ε, 2(^P_N) E^Q [(ξ-μ)^⊤X (ξ-μ) ]}. \eeIf $\varepsilon=0$ , then the robust MLE problem (4.3) reduces to the nominal MLE problem (4.3) because—by the definition of the sample mean and the sample covariance matrix—the (normal) nominal distribution has the same first- and second-order moments as the (discrete) empirical distribution and because the loss function in the expectation is quadratic in $\xi$ . One can show via Theorem 2.25 that (4.3) is equivalent to a convex SDP with a determinant term in the objective function. Provided that the Wasserstein radius $\varepsilon$ is strictly positive, this SDP is solvable even in the big-data regime when $m>N$ . Thus, it yields a valid precision matrix estimator even if the sample covariance matrix is rank-deficient. Moreover, as SDPs are tractable, the optimal estimator can be computed in polynomial time. In fact, the SDP at hand is highly symmetric and can therefore even be solved in closed form [68, Theorem 3.1].

Theorem 4.3 asserts that the robust MLE estimator $\mu^{\star}$ for the mean vector coincides with the sample mean $\widehat{\mu}$ . More interestingly, it further asserts that the robust MLE estimator $X^{\star}$ for the precision matrix has the same eigenvectors $v_{i}$ as the sample covariance matrix $\widehat{\Sigma}$ , while its eigenvalues $x^{\star}_{i}$ are obtained by applying the nonlinear transformation (36) to the corresponding eigenvalues $\lambda_{i}$ of $\widehat{\Sigma}$ . This transformation involves a single unknown parameter $\gamma^{\star}$ , which is the unique positive solution of the algebraic equation (36). As $X^{\star}$ is obtained by transforming the eigenvalues of the sample covariance matrix, it can be interpreted as a nonlinear shrinkage estimator. We thus refer to it as the Wasserstein shrinkage estimator.

As $X^{\star}$ and $\widehat{\Sigma}$ share the same eigenvectors, $X^{\star}$ is rotation-equivariant, that is, the estimator applied to the rotated data $R\,\widehat{\xi}_{i}$ , $i\in[N]$ , coincides with the rotated estimator $RX^{\star}R$ of the original data for every possible rotation matrix $R$ . Moreover, as all eigenvalues of $X^{\star}$ are strictly positive, the estimator is always invertible. Finally, Theorem 4.3 indicates that $X^{\star}$ can be computed highly efficiently by computing the spectral decomposition of $\widehat{\Sigma}$ and by solving the scalar algebraic equation (36), which can be accomplished by bisection.

One can show that the Wasserstein shrinkage estimator displays numerous desirable properties [68, Proposition 3.5]. First, its eigenvalues $x_{i}^{\star}$ decrease with $\varepsilon$ and eventually converge to 0. This makes intuitive sense as for large values of $\varepsilon$ nothing is known about $\xi$ , and thus the safest bet is that all of its components have high variance and low precision. Moreover, one can show that the order of the eigenvalues $x_{i}^{\star}$ matches the order of the inverse sample eigenvalues $1/\lambda_{i}$ irrespective of $\varepsilon>0$ , which is expected in the absence of any structural information. Finally, one can show that the condition number of $X^{\star}$ decreases monotonically to 1 as $\varepsilon$ grows. Thus, the condition number of $X^{\star}$ improves with the level of ambiguity.

A statistical theory that shows how to optimally choose $\varepsilon$ is developed in . Surprisingly, the Wasserstein radius that attains the lowest possible out-of-sample loss scales as $\varepsilon\propto 1/N$ instead of the canonical inverse square-root scaling, which may be expected for this problem.

So far we have assumed that there is no structural information about the distribution of $\xi$ besides normality. In some practical situation, however, the precision matrix $X$ may have a known sparsity pattern. Indeed, one can show that an element $X_{ij}$ of the precision matrix vanishes if and only if the random variables $\xi_{i}$ and $\xi_{j}$ are conditionally independent given all other components of $\xi$ . Conditional independencies of this type naturally arise, for example, in the analysis of spatio-temporal data. In the presence of sparsity information, the robust MLE problem is still equivalent to a tractable SDP. Even though it loses its analytical solvability, one can devise a tailored sequential quadratic approximation algorithm with rigorous convergence guarantees to solve the problem numerically, see [68, § 4].

4 Distributionally Robust Minimum Mean Square Error Estimation

Note that (37) constitutes an infinite-dimensional functional optimization problem and thus appears to be hard. However, by establishing a minimax theorem for (37) and exploiting the properties of elliptical distributions, one can show that the outer infimum in (37) is attained by an affine estimator. Combining this structural insight with Theorem 2.25 allows us to prove that the estimation problem (37) is in fact equivalent to a convex program .

If $\widehat{\Sigma}\succ 0$ , then the estimation problem (37) is equivalent to the nonlinear convex SDP

It is possible to eliminate all nonlinearities in (38) by using Schur complements and to reformulate the nonlinear convex SDP as a standard linear SDP, which is formally tractable. However, larger problem instances quickly exceed the capabilities of general-purpose solvers. Instead, there is merit in addressing the nonlinear SDP (38) directly with a customized first-order Frank-Wolfe algorithm, which starts at $S^{(0)}=\widehat{\Sigma}$ and constructs iterates

where $\gamma^{\star}$ is the unique solution with $\gamma^{\star}I\succ\nabla f(S^{(k)})$ of the algebraic equation

which can be solved via bisection [98, Theorem 3.2]. For a judiciously chosen step-size rule, the Frank-Wolfe algorithm also offers rigorous convergence guarantees [98, Theorem 3.3].

In some applications one has additional structural information about the relation between the signal $x$ and the observation $y$ (e.g., the measurement noise may be known to be independent of the signal, or the observation may be governed by a linear measurement model, etc.). Such structural information can be used to restrict the Wasserstein ambiguity set in (37), thereby reducing the conservativeness of the distributionally robust MMSE estimator .

5 Other Applications in Machine Learning

Ideas from distributionally robust optimization also permeate several other areas of statistics and machine learning. For example, a distributionally robust optimization model involving two Wasserstein balls centered at two distinct empirical distributions can be used to develop a computationally tractable convex approximation for the minimax robust hypothesis testing problem that aims to minimize the maximum of the worst-case type-I and type-II errors of a prescribed hypothesis test . Another example is data-driven inverse optimization, where one observes random signals as well as optimal solutions of an optimization problem parameterized by these signals. The aim is to predict the solution corresponding to a new unseen signal from $N$ independent historical observations without any knowledge of the optimization problem’s objective function. This problem can be framed as a structural regression problem that minimizes the worst-case expected prediction loss with respect to a Wasserstein ambiguity set over a space of candidate objective functions . Data-driven inverse optimization lends itself, for example, to learning the purchasing behavior of consumers, the production costs of electricity generators, the route choice preferences of passengers in a multimodal transportation system or the hidden optimality principles governing a biological system. As a third example, distributionally robust optimization models with Wasserstein ambiguity sets can be used to efficiently compute the worst-case misclassification probability of a given classifier, which amounts to evaluating the worst-case expectation of the (nonconvex) zero-one loss . Using similar techniques, one can also efficiently compute the worst-case probability of an undesirable event described by the conjunction or disjunction of several linear inequalities for the random vector $\xi$ . If the undesirable event can be influenced so as drive its worst-case probability below a prescribed tolerance, we face a distributionally robust chance constraint. Even though distributionally robust chance constrained programs with Wasserstein ambiguity sets around the empirical distribution are intractable in general, they are sometimes equivalent to mixed-integer linear programs that can be solved with off-the-shelf software . In contrast, distributionally robust chance constrained programs with moment ambiguity sets can often be reformulated as (or tightly approximated by) tractable conic programs .

To conclude, we highlight two opportunities for tailoring a distributionally robust decision problem with a Wasserstein ambiguity set around the empirical distribution to a given training dataset. Recall first that finite sample guarantees hold whenever $\varepsilon$ is large enough for the Wasserstein ball to contain the unknown data-generating distribution with high confidence $1-\beta$ . Recall also that the distributionally robust decision problem can often be reformulated as a tractable convex program whose size scales with the sample size $N$ . If the computational burden is unmanageable for the given sample size, we can select $K\ll N$ , approximate $\widehat{\mathds{P}}_{N}$ with the closest $K$ -point distribution $\mathds{Q}^{\star}_{K}$ in Wasserstein distance and replace the original Wasserstein ball of radius $\varepsilon$ around $\widehat{\mathds{P}}_{N}$ with a new inflated Wasserstein ball of radius $\varepsilon+W_{p}(\widehat{\mathds{P}}_{N},\mathds{Q}^{\star}_{K})$ around $\mathds{Q}^{\star}_{K}$ . By construction, the inflated Wasserstein ball contains the data-generating distribution with the same confidence $1-\beta$ . But the size of the corresponding decision problem is only proportional to $K$ . This approach provides a systematic method for reducing the computational burden without sacrificing robustness guarantees (but at the expense of increasing the model’s level of conservatism). The approximation of a rich $N$ -point distribution with a sparse $K$ -point distribution is referred to as scenario reduction in the stochastic programming literature. While the exact computation of $\mathds{Q}^{\star}_{K}$ is hard, there exist efficient approximation algorithms for scenario reduction .

An important input for any distributionally robust optimization model with a Wasserstein ambiguity set is the norm that determines the transportation cost in the definition of the Wasserstein distance. The flexibility to choose this norm could be exploited to improve the out-of-sample performance of the model’s optimizers. A method for learning the best Mahalanobis norm from the training data is described in . It is shown that this metric learning framework encompasses adaptive regularization as a special case.

Acknowledgments. This research was funded by the SNSF grant BSCGI0_157733.