Scenario Reduction Revisited: Fundamental Limits and Guarantees

TL;DR

This paper investigates the fundamental limits of scenario reduction using Wasserstein distance, identifying distributions most resistant to reduction, and introduces new approximation algorithms with guarantees for both discrete and continuous cases.

Contribution

It provides the first polynomial-time constant-factor approximation algorithms and exact exponential-time algorithms for scenario reduction under Wasserstein distance.

Findings

Identifies distributions with maximum Wasserstein distance to their reduced forms.

Provides sharp bounds on the benefits of continuous over discrete scenario reduction.

Introduces the first polynomial-time approximation algorithms and exact algorithms for scenario reduction.

Abstract

The goal of scenario reduction is to approximate a given discrete distribution with another discrete distribution that has fewer atoms. We distinguish continuous scenario reduction, where the new atoms may be chosen freely, and discrete scenario reduction, where the new atoms must be chosen from among the existing ones. Using the Wasserstein distance as measure of proximity between distributions, we identify those $n$-point distributions on the unit ball that are least susceptible to scenario reduction, i.e., that have maximum Wasserstein distance to their closest $m$-point distributions for some prescribed $m<n$. We also provide sharp bounds on the added benefit of continuous over discrete scenario reduction. Finally, to our best knowledge, we propose the first polynomial-time constant-factor approximations for both discrete and continuous scenario reduction as well as the first exact exponential-time algorithms for continuous scenario reduction.

Figures

Tables

Table 1. Table 1 : Runtimes of different methods for discrete scenario reduction.

	DPCV	LOC-1	LOC-2	MILP
Average (secs)	2.16	2.60	3.59	1,349.71
Worst-case (secs)	8.21	9.89	13.69	36,120.99

Equations

$$\displaystyle d_{l}(\mathbb{P},\mathbb{Q})=\left[\min_{\mathbf{\Pi}\in\mathbb{R}_{+}^{n\times m}}\left\{\sum_{i\in I}\sum_{j\in J}\pi_{ij}\|\bm{\xi}_{i}-\bm{\zeta}_{j}\|^{l}:\begin{array}[]{l}\sum_{j\in J}\pi_{ij}=p_{i}\ \forall i\in I\\[5.69054pt] \sum_{i\in I}\pi_{ij}=q_{j}\ \forall j\in J\end{array}\right\}\right]^{1/l},$$

$$C_l({\hat{\mathbb{P}}}_n,m)=\underset{\mathbb{Q}}{\min }\left\{d_l({\hat{\mathbb{P}}}_n,\mathbb{Q}):\mathbb{Q}\in \mathcal{P}({\mathbb{R}}^d,m)\right\},$$

$$D_l({\hat{\mathbb{P}}}_n,m)=\underset{\mathbb{Q}}{\min }\left\{d_l({\hat{\mathbb{P}}}_n,\mathbb{Q}):\mathbb{Q}\in \mathcal{P}(\text{supp}({\hat{\mathbb{P}}}_n),m)\right\},$$

$${\overline{C}}_l(n,m)=\underset{{\hat{\mathbb{P}}}_n\in {\mathcal{P}}_{\mathrm{E}}({\mathbb{R}}^d,n)}{\max }\left\{C_l({\hat{\mathbb{P}}}_n,m):\Vert \mathbf{\xi }{\Vert }_2\le 1\forall \mathbf{\xi }\in \text{supp}({\hat{\mathbb{P}}}_n)\right\},$$

$$\mathfrak{P}(I,m)=\big{\{}\{I_{1},\ldots,I_{m}\}\;:\;\emptyset\neq I_{1},\ldots,I_{m}\subseteq I,\;\;\cup_{j}I_{j}=I,\;\;I_{i}\cap I_{j}=\emptyset\;\;\forall i\neq j\big{\}},$$

$$C_l({\hat{\mathbb{P}}}_n,m)=\underset{\{I_j\}\in \mathfrak{P}(I,m)}{\min }{\left[\frac{1}{n}\sum _{j\in J}\underset{{\mathbf{\zeta }}_j\in {\mathbb{R}}^d}{\min }\sum _{i\in I_j}\Vert {\mathbf{\xi }}_i-{\mathbf{\zeta }}_j{\Vert }^l\right]}^{1/l}.$$

$$\underset{\mathbb{Q}\in \mathcal{P}(\mathrm{\Xi },\infty )}{\min }{d}_l({\hat{\mathbb{P}}}_n,\mathbb{Q})={\left[\frac{1}{n}\sum _{i\in I}\underset{\mathbf{\zeta }\in \mathrm{\Xi }}{\min }\Vert {\mathbf{\xi }}_i-\mathbf{\zeta }{\Vert }^l\right]}^{1/l},$$

$$C_l({\hat{\mathbb{P}}}_n,m)=\underset{\{{\mathbf{\zeta }}_j\}\subseteq {\mathbb{R}}^d}{\min }{\left[\frac{1}{n}\sum _{i\in I}\underset{j\in J}{\min }\Vert {\mathbf{\xi }}_i-{\mathbf{\zeta }}_j{\Vert }^l\right]}^{1/l}.$$

$${\mathbf{\zeta }}_j^{\star }=\text{gmed}(I_j)\in {\mathrm{argmin}}_{{\mathbf{\zeta }}_j\in {\mathbb{R}}^d}\sum _{i\in I_j}\Vert {\mathbf{\xi }}_i-{\mathbf{\zeta }}_j\Vert ,$$

$$\raisebox{7.11317pt}{$\overline{C}_{2}(n,m)\;=\;$}\begin{array}[]{c@{\quad}l}\displaystyle\max_{\{\bm{\xi}_{i}\}\,\subseteq\,\mathbb{R}^{d}}&\displaystyle\min_{\{I_{j}\}\in\mathfrak{P}(I,m)}\left[\frac{1}{n}\sum_{j\in J}\sum_{i\in I_{j}}\left\|\bm{\xi}_{i}-\text{mean}(I_{j})\right\|_{2}^{2}\right]^{1/2}\\[17.07164pt] \displaystyle\text{s.t.}&\displaystyle\|\bm{\xi}_{i}\|_{2}\leq 1\quad\forall i\in I.\end{array}$$

$$\raisebox{22.05087pt}{$\overline{C}_{2}^{2}(n,m)\;=\;$}\begin{array}[]{c@{\quad}l@{\quad}l}\displaystyle\max_{\tau\in\mathbb{R},\;\{\bm{\xi}_{i}\}\,\subseteq\,\mathbb{R}^{d}}&\displaystyle\frac{1}{n}\tau\\[5.69054pt] \displaystyle\text{s.t.}&\displaystyle\tau\leq\sum_{j\in J}\sum_{i\in I_{j}}\left\|\bm{\xi}_{i}-\text{mean}(I_{j})\right\|_{2}^{2}&\displaystyle\forall\{I_{j}\}\in\mathfrak{P}(I,m)\\[14.22636pt] &\displaystyle\bm{\xi}_{i}{}^{\top}\bm{\xi}_{i}\leq 1&\forall i\in I.\end{array}$$

$$\mspace{-20.0mu}\begin{aligned} \big{\|}\bm{\xi}_{i}-&\text{mean}(I_{j})\big{\|}_{2}^{2}\;\;=\;\;\frac{1}{|I_{j}|^{2}}\Big{\|}|I_{j}|\bm{\xi}_{i}-\sum_{k\in I_{j}}\bm{\xi}_{k}\Big{\|}_{2}^{2}\\ &=\;\;\frac{1}{|I_{j}|^{2}}\Bigg{(}|I_{j}|^{2}\bm{\xi}_{i}{}^{\top}\bm{\xi}_{i}-2|I_{j}|\sum_{k\in I_{j}}\bm{\xi}_{i}{}^{\top}\bm{\xi}_{k}+\sum_{k\in I_{j}}\bm{\xi}_{k}{}^{\top}\bm{\xi}_{k}+\sum_{\begin{subarray}{c}k,k^{\prime}\in I_{j}\\ k\neq k^{\prime}\end{subarray}}\bm{\xi}_{k}{}^{\top}\bm{\xi}_{k^{\prime}}\Bigg{)}.\end{aligned}$$

$$\mathbf{S}=[{\mathbf{\xi }}_1,\dots ,{\mathbf{\xi }}_n{]}^{\top }[{\mathbf{\xi }}_1,\dots ,{\mathbf{\xi }}_n]\in {\mathbb{S}}^n,\mathbf{S}⪰0\text{ and }\text{rank}(\mathbf{S})\le \min \{n,d\}$$

$$\tau\leq\sum_{j\in J}\frac{1}{|I_{j}|^{2}}\sum_{i\in I_{j}}\Bigg{(}|I_{j}|^{2}s_{ii}-2|I_{j}|\sum_{k\in I_{j}}s_{ik}+\sum_{k\in I_{j}}s_{kk}+\sum_{\begin{subarray}{c}k,k^{\prime}\in I_{j}\\ k\neq k^{\prime}\end{subarray}}s_{kk^{\prime}}\Bigg{)}.$$

$$\begin{array}[]{c@{\quad}l@{\quad}l}\displaystyle\max_{\tau\in\mathbb{R},\;\mathbf{S}\in\mathbb{S}^{n}}&\displaystyle\frac{1}{n}\tau\\ \displaystyle\text{s.t.}&\displaystyle\tau\leq\sum_{j\in J}\frac{1}{|I_{j}|^{2}}\sum_{i\in I_{j}}\Bigg{(}|I_{j}|^{2}s_{ii}-2|I_{j}|\sum_{k\in I_{j}}s_{ik}+\sum_{k\in I_{j}}s_{kk}+\sum_{\begin{subarray}{c}k,k^{\prime}\in I_{j}\\ k\neq k^{\prime}\end{subarray}}s_{kk^{\prime}}\Bigg{)}\\[19.91692pt] &\displaystyle\mspace{405.0mu}\forall\{I_{j}\}\in\mathfrak{P}(I,m)\\[2.84526pt] &\displaystyle\mathbf{S}\succeq\bm{0},\;\;s_{ii}\leq 1\quad\forall i\in I,\end{array}$$

$$\begin{array}[]{c@{\quad}l@{\quad}l}\displaystyle\max_{\tau,\alpha,\beta\in\mathbb{R}}&\displaystyle\frac{1}{n}\tau\\ \displaystyle\text{s.t.}&\displaystyle\tau\leq(n-m)\alpha,\;\;\alpha+\beta\leq 1\\[2.84526pt] &\displaystyle\alpha\geq 0,\;\;\alpha+n\beta\geq 0,\end{array}$$

$$\frac{1}{|I_{j}|^{2}}\sum_{i\in I_{j}}\Bigg{(}|I_{j}|^{2}(\alpha+\beta)-2|I_{j}|\left(\alpha+|I_{j}|\beta\right)+\sum_{k\in I_{j}}(\alpha+\beta)+\sum_{\begin{subarray}{c}k,k^{\prime}\in I_{j}\\ k\neq k^{\prime}\end{subarray}}\beta\Bigg{)}=(|I_{j}|-1)\alpha,$$

$${\mathbf{\xi }}_i=y\mathbf{e}+(x-y){\mathbf{e}}_i\in {\mathbb{R}}^n\text{with}x=\sqrt{\frac{n-1}{n}}\text{ and }y=\frac{-1}{\sqrt{n(n-1)}}.$$

$$\raisebox{11.38109pt}{$\overline{C}_{1}(n,m)\;=\;$}\begin{array}[]{c@{\quad}l}\displaystyle\max_{\{\bm{\xi}_{i}\}\,\subseteq\,\mathbb{R}^{d}}&\displaystyle\min_{\{I_{j}\}\in\mathfrak{P}(I,m)}\frac{1}{n}\sum_{j\in J}\sum_{i\in I_{j}}\left\|\bm{\xi}_{i}-\text{gmed}(I_{j})\right\|_{2}\\[14.22636pt] \displaystyle\text{s.t.}&\displaystyle\|\bm{\xi}_{i}\|_{2}\leq 1\quad\forall i\in I.\end{array}$$

$${\overline{C}}_1(n,m)$$

$$\displaystyle\raisebox{7.8245pt}{$\leq$}\;\;\begin{array}[]{c@{\quad}l}\displaystyle\max_{\{\bm{\xi}_{i}\}\,\subseteq\,\mathbb{R}^{d}}&\displaystyle\min_{\{I_{j}\}\in\mathfrak{P}(I,m)}\left[\frac{1}{n}\sum_{j\in J}\sum_{i\in I_{j}}\left\|\bm{\xi}_{i}-\text{mean}(I_{j})\right\|_{2}^{2}\right]^{1/2}\\[17.07164pt] \displaystyle\text{s.t.}&\displaystyle\|\bm{\xi}_{i}\|_{2}\leq 1\quad\forall i\in I,\end{array}$$

$$\sum _{i\in I_j}{\Vert {\mathbf{\xi }}_i-\text{gmed}(I_j)\Vert }_2\le \sum _{i\in I_j}{\Vert {\mathbf{\xi }}_i-\text{mean}(I_j)\Vert }_2\forall j\in J,$$

$${\Vert {\mathbf{\xi }}_i-\text{mean}(I_j)\Vert }_2$$

$$={\left[{\left(x-\frac{x+(|I_j|-1)y}{|I_j|}\right)}^2+(|I_j|-1){\left(y-\frac{x+(|I_j|-1)y}{|I_j|}\right)}^2\right]}^{1/2}$$

$$={\left[\frac{|I_j|-1}{|I_j|}\right]}^{1/2}(x-y)={\left[\frac{n(|I_j|-1)}{(n-1)|I_j|}\right]}^{1/2}\forall i\in I_j,$$

$$C_1({\hat{\mathbb{P}}}_n,m)$$

$$=\underset{\{I_j\}\in \mathfrak{P}(I,m)}{\min }\frac{1}{\sqrt{n(n-1)}}\sum _{j\in J}\sqrt{|I_j|(|I_j|-1)}.$$

$$C_1({\hat{\mathbb{P}}}_n,m)$$

$$\sum _{j\in J}\sqrt{z_j(z_j+1)}$$

$$\displaystyle\geq\;\Bigg{[}\sum_{j\in J}z_{j}(z_{j}+1)+\sum_{\begin{subarray}{c}j,j^{\prime}\in J\\ j\neq j^{\prime}\end{subarray}}z_{j}z_{j^{\prime}}\Bigg{]}^{1/2}$$

$$\displaystyle=\;\left[\Bigg{(}\sum_{j\in J}z_{j}\Bigg{)}^{2}+\sum_{j\in J}z_{j}\right]^{1/2}\;=\;\sqrt{(n-m)(n-m+1)},$$

$${\overline{C}}_2^2(n,m)\le {\overline{C}}_1(n,m)\le {\overline{C}}_2(n,m).$$

$$C_l({\hat{\mathbb{P}}}_n,m)=r\cdot C_l\left(\frac{1}{n}\sum _{i\in I}{\delta }_{\frac{{\mathbf{\xi }}_i-\mathbf{\mu }}{r}},m\right)\le r\cdot {\overline{C}}_l(n,m)\lesssim r\cdot \sqrt{1-p},$$

Full text

1 Introduction

The vast majority of numerical solution schemes in stochastic programming rely on a discrete approximation of the true (typically continuous) probability distribution governing the uncertain problem parameters. This discrete approximation is often generated by sampling from the true distribution. Alternatively, it could be constructed directly from real historical observations of the uncertain parameters. To obtain a faithful approximation for the true distribution, however, the discrete distribution must have a large number $n$ of support points or scenarios, which may render the underlying stochastic program computationally excruciating.

An effective means to ease the computational burden is to rely on scenario reduction pioneered by Dupačová et al (2003), which aims to approximate the initial $n$-point distribution with a simpler $m$-point distribution ($m<n$) that is as close as possible to the initial distribution with respect to a probability metric; see also Heitsch and Römisch (2003). The modern stability theory of stochastic programming surveyed by Dupačová (1990) and Römisch (2003) indicates that the Wasserstein distance may serve as a natural candidate for this probability metric.

Our interest in Wasserstein distance-based scenario reduction is also fuelled by recent progress in data-driven distributionally robust optimization, where it has been shown that the worst-case expectation of an uncertain cost over all distributions in a Wasserstein ball can often be computed efficiently via convex optimization (Mohajerin Esfahani and Kuhn 2015, Zhao and Guan 2015, Gao and Kleywegt 2016). A Wasserstein ball is defined as the family of all distributions that are within a certain Wasserstein distance from a discrete reference distribution. As distributionally robust optimization problems over Wasserstein balls are harder to solve than their stochastic counterparts, we expect significant computational savings from replacing the initial $n$-point reference distribution with a new $m$-point reference distribution. The benefits of scenario reduction may be particularly striking for two-stage distributionally robust linear programs, which admit tight approximations as semidefinite programs (Hanasusanto and Kuhn 2016).

Suppose now that the initial distribution is given by $\mathbb{P}=\sum_{i\in I}p_{i}\delta_{\bm{\xi}_{i}}$, where $\bm{\xi}_{i}\in\mathbb{R}^{d}$ and $p_{i}\in[0,1]$ represent the location and probability of the $i$-th scenario of $\mathbb{P}$ for $i\in I=\{1,\ldots,n\}$. Similarly, assume that the reduced target distribution is representable as $\mathbb{Q}=\sum_{j\in J}q_{j}\delta_{\bm{\zeta}_{j}}$, where $\bm{\zeta}_{j}\in\mathbb{R}^{d}$ and $q_{j}\in[0,1]$ stand for the location and probability of the $j$-th scenario of $\mathbb{Q}$ for $j\in J=\{1,\ldots,m\}$. Then, the type-$l$ Wasserstein distance between $\mathbb{P}$ and $\mathbb{Q}$ is defined through

[TABLE]

where $l\geq 1$ and $\|\cdot\|$ denotes some norm on $\mathbb{R}^{d}$, see, e.g., Heitsch and Römisch (2007) or Pflug and Pichler (2011). The linear program in the definition of the Wasserstein distance can be viewed as a minimum-cost transportation problem, where $\pi_{ij}$ represents the amount of probability mass shipped from $\bm{\xi}_{i}$ to $\bm{\zeta}_{j}$ at unit transportation cost $\|\bm{\xi}_{i}-\bm{\zeta}_{j}\|^{l}$. Thus, $d_{l}^{l}(\mathbb{P},\mathbb{Q})$ quantifies the minimum cost of moving the initial distribution $\mathbb{P}$ to the target distribution $\mathbb{Q}$.

For any $\mathrm{\Xi}\subseteq\mathbb{R}^{d}$, we denote by $\mathcal{P}_{\mathrm{E}}(\mathrm{\Xi},n)$ the set of all uniform discrete distributions on $\mathrm{\Xi}$ with exactly $n$ distinct scenarios and by $\mathcal{P}(\mathrm{\Xi},m)$ the set of all (not necessarily uniform) discrete distributions on $\mathrm{\Xi}$ with at most $m$ scenarios. We henceforth assume that $\mathbb{P}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$. This assumption is crucial for the simplicity of the results in Sections 2 and 3, and it is almost surely satisfied whenever $\mathbb{P}$ is obtained via sampling from a continuous probability distribution. Hence, we can think of $\mathbb{P}$ as an empirical distribution. To remind us of this interpretation, we will henceforth denote the initial distribution by $\hat{\mathbb{P}}_{n}$. Note that the pairwise difference of the scenarios can always be enforced by slightly perturbing their locations, while the uniformity of their probabilities can be enforced by decomposing the scenarios into clusters of close but mutually distinct sub-scenarios with (smaller) uniform probabilities.

We are now ready to introduce the continuous scenario reduction problem

[TABLE]

where the new scenarios $\bm{\zeta}_{j}$, $j\in J$, of the target distribution $\mathbb{Q}$ may be chosen freely from within $\mathbb{R}^{d}$, as well as the discrete scenario reduction problem

[TABLE]

where the new scenarios must be chosen from within the support of the empirical distribution, which is given by the finite set $\text{supp}(\hat{\mathbb{P}}_{n})=\{\bm{\xi}_{i}:i\in I\}$. Even though the continuous scenario reduction problem offers more flexibility and is therefore guaranteed to find (weakly) better approximations to the initial empirical distribution, to our best knowledge, the existing stochastic programming literature has exclusively focused on the discrete scenario reduction problem.

Note that if the support points $\bm{\zeta}_{j}$, $j\in J$, are fixed, then both scenario reduction problems simplify to a linear program over the probabilities $q_{j}$, $j\in J$, which admits an explicit solution (Dupačová et al 2003, Theorem 2). Otherwise, however, both problems are intractable. Indeed, if $l=1$, then the discrete scenario reduction problem represents a metric $k$-median problem with $k=m$, which was shown to be $\mathcal{NP}$-hard by Kariv and Hakimi (1979). If $l=2$ and distances in $\mathbb{R}^{d}$ are measured by the 2-norm, on the other hand, then the continuous scenario reduction problem constitutes a $k$-means clustering problem with $k=m$, which is $\mathcal{NP}$-hard even if $d=2$ or $m=2$; see Mahajan et al (2009) and Aloise et al (2009).

Heitsch and Römisch (2003) have shown that the discrete scenario reduction problem admits a reformulation as a mixed-integer linear program (MILP), which can be solved to global optimality for $n\lesssim 10^{3}$ using off-the-shelf solvers. For larger instances, however, one must resort to approximation algorithms. Most large-scale discrete scenario reduction problems are nowadays solved with a greedy heuristic that was originally devised by Dupačová et al (2003) and further refined by Heitsch and Römisch (2003). For example, this heuristic is routinely used for scenario (tree) reduction in the context of power systems operations, see, e.g., Römisch and Vigerske (2010) or Morales et al (2009) and the references therein. Despite its practical success, we will show in Section 4 that this heuristic fails to provide a constant-factor approximation for the discrete scenario reduction problem.

This paper extends the theory of scenario reduction along several dimensions.

• (i)

We establish fundamental performance guarantees for continuous scenario reduction when $l\in\{1,2\}$, i.e., we show that the Wasserstein distance of the initial $n$-point distribution to its nearest $m$-point distribution is bounded by $\sqrt{\frac{n-m}{n-1}}$ across all initial distributions on the unit ball in $\mathbb{R}^{d}$. We show that for $l=2$ this worst-case performance is attained by some initial distribution, which we construct explicitly. We also provide evidence indicating that this worst-case performance reflects the norm rather than the exception in high dimensions $d$. Finally, we provide a lower bound on the worst-case performance for $l=1$.

• (ii)

We analyze the loss of optimality incurred by solving the discrete scenario reduction problem instead of its continuous counterpart. Specifically, we demonstrate that the ratio $D_{l}(\hat{\mathbb{P}}_{n},m)/C_{l}(\hat{\mathbb{P}}_{n},m)$ is bounded by $\sqrt{2}$ for $l=2$ and by 2 for $l=1$. We also show that these bounds are essentially tight.

• (iii)

We showcase the intimate relation between scenario reduction and $k$-means clustering. By leveraging existing constant-factor approximation algorithms for $k$-median clustering problems due to Arya et al (2004) and the new performance bounds from (ii), we develop the first polynomial-time constant-factor approximation algorithms for both continuous and discrete scenario reduction. We also show that these algorithms can be warmstarted using the greedy heuristic by Dupačová et al (2003) to improve practical performance.

• (iv)

We present exact mixed-integer programming reformulations for the continuous scenario reduction problem.

Continuous scenario reduction is intimately related to the optimal quantization of probability distributions, where one seeks an $m$-point distribution approximating a non-discrete initial distribution. Research efforts in this domain have mainly focused on the asymptotic behavior of the quantization problem as $m$ tends to infinity, see Graf and Luschgy (2000). The ramifications of this stream of literature for stochastic programming are discussed by Pflug and Pichler (2011). Techniques familiar from scenario reduction lend themselves also for scenario generation, where one aims to construct a scenario tree with a prescribed branching structure that approximates a given stochastic process with respect to a probability metric, see, e.g., Pflug (2001) and Hochreiter and Pflug (2007).

The rest of this paper unfolds as follows. Section 2 seeks to identify $n$-point distributions on the unit ball that are least susceptible to scenario reduction, i.e., that have maximum Wasserstein distance to their closest $m$-point distributions, and Section 3 discusses sharp bounds on the added benefit of continuous over discrete scenario reduction. Section 4 presents exact exponential-time algorithms as well as polynomial-time constant-factor approximations for scenario reduction. Section 5 reports on numerical results for a color quantization experiment. Unless otherwise specified, below we will always work with the 2-norm on $\mathbb{R}^{d}$.

Notation:

We let $\mathbb{I}$ be the identity matrix, $\mathbf{e}$ the vector of all ones and $\mathbf{e}_{i}$ the $i$-th standard basis vector of appropriate dimensions. The $ij$-th element of a matrix $\mathbf{A}$ is denoted by $a_{ij}$. For $\mathbf{A}$ and $\mathbf{B}$ in the space $\mathbb{S}^{n}$ of symmetric $n\times n$ matrices, the relation $\mathbf{A}\succeq\mathbf{B}$ means that $\mathbf{A}-\mathbf{B}$ is positive semidefinite. Generic norms are denoted by $\|\cdot\|$, while $\|\cdot\|_{p}$ stands for the $p$-norm, $p\geq 1$. For $\mathrm{\Xi}\subseteq\mathbb{R}^{d}$, we define $\mathcal{P}(\mathrm{\Xi},m)$ as the set of all probability distributions supported on at most $m$ points in $\mathrm{\Xi}$ and $\mathcal{P}_{\mathrm{E}}(\mathrm{\Xi},n)$ as the set of all uniform distributions supported on exactly $n$ distinct points in $\mathrm{\Xi}$. The support of a probability distribution $\mathbb{P}$ is denoted by $\text{supp}(\mathbb{P})$, and the Dirac distribution concentrating unit mass at $\bm{\xi}$ is denoted by $\delta_{\bm{\xi}}$.

2 Fundamental Limits of Scenario Reduction

In this section we characterize the Wasserstein distance $C_{l}(\hat{\mathbb{P}}_{n},m)$ between an $n$-point empirical distribution $\hat{\mathbb{P}}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{\bm{\xi}_{i}}$ and its continuously reduced optimal $m$-point distribution $\mathbb{Q}\in\mathcal{P}(\mathbb{R}^{d},m)$. Since the positive homogeneity of the Wasserstein distance $d_{l}$ implies that $C_{l}(\hat{\mathbb{P}}^{\prime}_{n},m)=\lambda\cdot C_{l}(\hat{\mathbb{P}}_{n},m)$ for the scaled distribution $\hat{\mathbb{P}}^{\prime}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{\lambda\bm{\xi}_{i}}$, $\lambda\in\mathbb{R}_{+}$, we restrict ourselves to empirical distributions $\hat{\mathbb{P}}_{n}$ whose scenarios satisfy $\left\lVert\bm{\xi}_{i}\right\rVert_{2}\leq 1$, $i=1,\ldots,n$. We thus want to quantify

[TABLE]

which amounts to the worst-case (i.e., largest) Wasserstein distance between any $n$-point empirical distribution $\hat{\mathbb{P}}_{n}$ over the unit ball and its optimally selected continuous $m$-point scenario reduction. By construction, this worst-case distance satisfies $\overline{C}_{l}(n,m)\geq 0$, and the lower bound is attained whenever $n=m$. One also verifies that $\overline{C}_{l}(n,m)\leq\overline{C}_{l}(n,1)\leq 1$ since the Wasserstein distance to the Dirac distribution $\delta_{\bm{0}}$ is bounded above by $1$. Our goal is to derive possibly tight upper bounds on $\overline{C}_{l}(n,m)$ for the Wasserstein distances of type $l\in\{1,2\}$.

In the following, we denote by $\mathfrak{P}(I,m)$ the family of all $m$-set partitions of the index set $I$, i.e.,

[TABLE]

and an element of this set (i.e. a specific $m$-set partition) as $\{I_{j}\}\in\mathfrak{P}(I,m)$. Our derivations will make extensive use of the following theorem.

Theorem 2.1

For any type-$l$ Wasserstein distance induced by any norm $\|\cdot\|$, the continuous scenario reduction problem can be reformulated as

[TABLE]

Problem (2) can be interpreted as a Voronoi partitioning problem that asks for a Voronoi decomposition of $\mathbb{R}^{d}$ into $m$ cells whose Voronoi centroids $\bm{\zeta}_{1},\ldots,\bm{\zeta}_{m}$ minimize the cumulative $l$-th powers of the distances to $n$ prespecified points $\bm{\xi}_{1},\ldots,\bm{\xi}_{n}$.

Proof

of Theorem 2.1   Theorem 2 of Dupačová et al (2003) implies that the smallest Wasserstein distance between the empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ and any distribution $\mathbb{Q}$ supported on a finite set $\mathrm{\Xi}\subset\mathbb{R}^{d}$ amounts to

[TABLE]

where $\mathcal{P}(\mathrm{\Xi},\infty)$ denotes the set of all probability distributions supported on the finite set $\mathrm{\Xi}$. The continuous scenario reduction problem $C_{l}(\hat{\mathbb{P}}_{n},m)$ selects the set $\mathrm{\Xi}^{\star}$ that minimizes this quantity over all sets in $\mathrm{\Xi}\subset\mathbb{R}^{d}$ with $|\mathrm{\Xi}|=m$ elements:

[TABLE]

One readily verifies that any optimal solution $\{\bm{\zeta}_{1}^{\star},\ldots,\bm{\zeta}_{m}^{\star}\}$ to problem (3) corresponds to an optimal solution $\{I_{1}^{\star},\ldots,I_{m}^{\star}\}$ to problem (2) with the same objective value if we identify the set $I_{j}^{\star}$ with all observations $\bm{\xi}_{i}$ that are closer to $\bm{\zeta}_{j}^{\star}$ than any other $\bm{\zeta}_{j^{\prime}}^{\star}$ (ties may be broken arbitrarily). Likewise, any optimal solution $\{I_{1}^{\star},\ldots,I_{m}^{\star}\}$ to problem (2) with inner minimizers $\{\bm{\zeta}_{1}^{\star},\ldots,\bm{\zeta}_{m}^{\star}\}$ translates into an optimal solution $\{\bm{\zeta}_{1}^{\star},\ldots,\bm{\zeta}_{m}^{\star}\}$ to problem (3) with the same objective value. ∎

Remark 1

(Minimizers of (2))

For $l=2$, the inner minimum corresponding to the set $I_{j}$ is attained by the mean $\bm{\zeta}^{\star}_{j}=\text{mean}(I_{j})=\frac{1}{|I_{j}|}\sum_{i\in I_{j}}\bm{\xi}_{i}$. Likewise, for $l=1$, the inner minimum corresponding to the set $I_{j}$ is attained by any geometric median

[TABLE]

which can be determined efficiently by solving a second-order cone program whenever a $p$-norm with rational $p\geq 1$ is considered (Alizadeh and Goldfarb (2003)).

The rest of this section derives tight upper bounds on $\overline{C}_{l}(n,m)$ for Wasserstein distances of type $l=2$ (Section 2.1) as well as upper and lower bounds for Wasserstein distances of type $l=1$ (Section 2.2). We summarize and discuss our findings in Section 2.3.

2.1 Fundamental Limits for the Type-2 Wasserstein Distance

We now derive a revised upper bound on $\overline{C}_{l}(n,m)$ for the type-$2$ Wasserstein distance. The result relies on auxiliary lemmas that are relegated to the appendix.

Theorem 2.2

The worst-case type-$2$ Wasserstein distance satisfies $\overline{C}_{2}(n,m)\leq\sqrt{\frac{n-m}{n-1}}$.

Note that whenever the reduced distribution satisfies $m>1$, the bound of Theorem 2.2 is strictly tighter than the naïve bound of $1$ from the previous section.

Proof

of Theorem 2.2   From Theorem 2.1 and Remark 1 we observe that

[TABLE]

Introducing the epigraphical variable $\tau$, this problem can be expressed as

[TABLE]

For each $j\in J$ and $i\in I_{j}$, the squared norm in the first constraint of (4) can be expressed in terms of the inner products between pairs of empirical observations:

[TABLE]

Introducing the Gram matrix

[TABLE]

then allows us to simplify the first constraint in (4) to

[TABLE]

Note that the second constraint in problem (4) can now be expressed as $s_{ii}\leq 1$, and hence all constraints in (4) are linear in the Gram matrix $\mathbf{S}$.

Our discussion implies that we obtain an upper bound on $\overline{C}_{2}(n,m)$ by reformulating problem (4) as a semidefinite program in terms of the Gram matrix $\mathbf{S}$

[TABLE]

where we have relaxed the rank condition in the definition of the Gram matrix (5). Lemma 1 in the appendix shows that (6) has an optimal solution $(\tau^{\star},\mathbf{S}^{\star})$ that satisfies $\mathbf{S}^{\star}=\alpha\mathbb{I}+\beta\mathbf{e}\mathbf{e}^{\top}$ for some $\alpha,\beta\in\mathbb{R}$. Moreover, Lemma 2 in the appendix shows that any matrix of the form $\mathbf{S}=\alpha\mathbb{I}+\beta\bm{1}\bm{1}^{\top}$ is positive semidefinite if and only if $\alpha\geq 0$ and $\alpha+n\beta\geq 0$. We thus conclude that (6) can be reformulated as

[TABLE]

where the first constraint follows from the fact that for any set $I_{j}$ in (6), we have

[TABLE]

and $\sum_{j\in J}(|I_{j}|-1)\alpha=(n-m)\alpha$ since $|I|=n$ and $|J|=m$. The statement of the theorem now follows since problem (7) is optimized by $\tau^{\star}=\frac{n(n-m)}{(n-1)}$, $\alpha^{\star}=\frac{n}{n-1}$ and $\beta^{\star}=\frac{-1}{n-1}$. ∎

The proof of Theorem 2.2 shows that the upper bound $\sqrt{\frac{n-m}{n-1}}$ on the worst-case type-$2$ Wasserstein distance $\overline{C}_{2}(n,m)$ is tight whenever there is an empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ whose scenarios $\bm{\xi}_{1},\ldots,\bm{\xi}_{n}$ correspond to a Gram matrix $\mathbf{S}=[\bm{\xi}_{1},\ldots,\bm{\xi}_{n}]^{\top}[\bm{\xi}_{1},\ldots,\bm{\xi}_{n}]=\frac{n}{n-1}\mathbb{I}-\frac{1}{n-1}\mathbf{e}\mathbf{e}^{\top}$, which implies $\|\bm{\xi}_{i}\|_{2}=\sqrt{s_{ii}}=1$ for all $i\in I$. We now show that such an empirical distribution exists when $d\geq n-1$.

Proposition 1

For $d\geq n-1$, there is $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ with $\left\lVert\bm{\xi}\right\rVert_{2}\leq 1$ for all $\bm{\xi}\in\text{\emph{supp}}(\hat{\mathbb{P}}_{n})$ such that $C_{2}(\hat{\mathbb{P}}_{n},m)=\sqrt{\frac{n-m}{n-1}}$.

Proof

Assume first that $d=n$ and consider the empirical distribution $\hat{\mathbb{P}}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{\bm{\xi}_{i}}$ defined through

[TABLE]

A direct calculation reveals that $\mathbf{S}=[\bm{\xi}_{1},\ldots,\bm{\xi}_{n}]^{\top}[\bm{\xi}_{1},\ldots,\bm{\xi}_{n}]=\frac{n}{n-1}\mathbb{I}-\frac{1}{n-1}\mathbf{e}\mathbf{e}^{\top}$.

To prove the statement for $d=n-1$, we note that the $n$ scenarios in (8) lie on the $(n-1)$-dimensional subspace $\mathcal{H}$ orthogonal to $\mathbf{e}\in\mathbb{R}^{n}$. Thus, there exists a rotation that maps $\mathcal{H}$ to $\mathbb{R}^{n-1}\times\{0\}$. As the Gram matrix is invariant under rotations, the rotated scenarios give rise to an empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{n-1},n)$ satisfying the statement of the proposition. Likewise, for $d>n$ the linear transformation $\bm{\xi}_{i}\mapsto(\mathbb{I},\bm{0})^{\top}\bm{\xi}_{i}$, $\mathbb{I}\in\mathbb{R}^{n\times n}$ and $\bm{0}\in\mathbb{R}^{n\times(d-n)}$, generates an empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ that satisfies the statement of the proposition. ∎

Proposition 1 requires that $d\geq n-1$, which appears to be restrictive. We note, however, that this condition is only sufficient (and not necessary) to guarantee the tightness of the bound from Theorem 2.2. Moreover, we will observe in Section 2.3 that the bound of Theorem 2.2 provides surprisingly accurate guidance for the Wasserstein distance between practice-relevant empirical distributions $\hat{\mathbb{P}}_{n}$ and their continuously reduced optimal distributions.

2.2 Fundamental Limits for the Type-1 Wasserstein Distance

In analogy to the previous section, we now derive a revised upper bound on $\overline{C}_{l}(n,m)$ for the type-$1$ Wasserstein distance.

Theorem 2.3

The worst-case type-$1$ Wasserstein distance satisfies $\overline{C}_{1}(n,m)\leq\sqrt{\frac{n-m}{n-1}}$.

Note that this bound is identical to the bound of Theorem 2.2 for $l=2$.

Proof

of Theorem 2.3   Leveraging again Theorem 2.1 and Remark 1, we obtain that

[TABLE]

We show that $\overline{C}_{1}(n,m)\leq\overline{C}_{2}(n,m)$ for all $n$ and $m=1,\ldots,n$, which in turn proves the statement of the theorem by virtue of Theorem 2.2. To this end, we observe that

[TABLE]

where the first inequality follows from the definition of the geometric median, which ensures that

[TABLE]

and the second inequality is due to the arithmetic-mean quadratic-mean inequality (Steele 2004, Exercise 2.14). The statement of the theorem now follows from the observation that the optimal value of (14) is identical to $\overline{C}_{2}(n,m)$. ∎

In the next proposition we derive a lower bound on $\overline{C}_{1}(n,m)$.

Proposition 2

For $d\geq n-1$, the worst-case type-1 Wasserstein distance satisfies $\overline{C}_{1}(n,m)\geq\sqrt{\frac{(n-m)(n-m+1)}{n(n-1)}}$.

Proof

Assume first that $d=n$, and consider the empirical distribution $\hat{\mathbb{P}}_{n}$ with scenarios defined as in (8). Let $\{I_{j}\}\in\mathfrak{P}(I,m)$ be an arbitrary $m$-set partition of $I$ and note that $\text{gmed}(I_{j})=\text{mean}(I_{j})$ for every $j\in J$ due to the permutation symmetry of the $\bm{\xi}_{i}$. This is indeed the case because $\bm{0}\in\partial f_{j}(\text{mean}(I_{j}))$ for each $f_{j}(\bm{\zeta})=\sum_{i\in I_{j}}\left\|y\mathbf{e}+(x-y)\mathbf{e}_{i}-\bm{\zeta}\right\|_{2}$, $j\in J$. Thus, we have

[TABLE]

where the last equality follows from the definitions of $x$ and $y$ in (8). By Theorem 2.1 and Remark 1 we therefore obtain

[TABLE]

By introducing auxiliary variables $z_{j}=|I_{j}|-1\in\mathbb{N}_{0}$, $j\in J$, we find that determining $C_{1}(\hat{\mathbb{P}}_{n},m)$ is tantamount to solving

[TABLE]

Observe that the objective function of $z_{1}=n-m$ and $z_{2}=\ldots=z_{m}=0$ evaluates to $\sqrt{(n-m)(n-m+1)}$, which implies that $C_{1}(\hat{\mathbb{P}}_{n},m)\leq\sqrt{\frac{(n-m)(n-m+1)}{n(n-1)}}$. Hence, it remains to establish the reverse inequality. To this end, we note that

[TABLE]

and thus the claim follows for $d=n$. The cases $d=n-1$ and $d>n$ can be reduced to the case $d=n$ as in Proposition 1. Details are omitted for brevity. ∎

Proposition 2 asserts that $\overline{C}_{1}(n,m)\gtrsim\frac{n-m}{n-1}=\overline{C}{}^{2}_{2}(n,m)$ whenever $d\geq n-1$. Together with Theorem 2.3, we thus obtain the following relation between the worst-case Wasserstein distances of types $l=1$ and $l=2$:

[TABLE]

We conjecture that the lower bound is tighter, but we were not able to prove this.

2.3 Discussion

Theorems 2.2 and 2.3 imply that $\overline{C}_{l}(n,m)\lesssim\sqrt{1-p}$ for large $n$ and for $l\in\{1,2\}$, where $p=\frac{m}{n}$ represents the desired reduction factor. The significance of this result is that it offers a priori guidelines for selecting the number $m$ of support points in the reduced distribution. To see this, consider any empirical distribution $\hat{\mathbb{P}}_{n}=\frac{1}{n}\sum_{i\in I}\delta_{\bm{\xi}_{i}}$, and denote by $r\geq 0$ and $\bm{\mu}\in\mathbb{R}^{d}$ the radius and the center of any (ideally the smallest) ball enclosing $\bm{\xi}_{1},\ldots,\bm{\xi}_{n}$, respectively. In this case, we have

[TABLE]

where the inequality holds because $\|(\bm{\xi}_{i}-\bm{\mu})/r\|_{2}\leq 1$ for every $i\in I$. Note that (15) enables us to find an upper bound on the smallest $m$ guaranteeing that $C_{l}(\hat{\mathbb{P}}_{n},m)$ falls below a prescribed threshold (i.e., guaranteeing that the reduced $m$-point distribution remains within some prescribed distance from $\hat{\mathbb{P}}_{n}$).

Even though the inequality in (15) can be tight, which has been established in Proposition 1, one might suspect that typically $C_{l}(\hat{\mathbb{P}}_{n},m)$ is significantly smaller than $r\cdot\sqrt{1-p}$ when the points $\bm{\xi}_{i}\in\mathbb{R}^{d}$, $i\in I$, are sampled randomly from a standard distribution, e.g., a multivariate uniform or normal distribution. However, while the upper bound (15) can be loose for low-dimensional data, Proposition 3 below suggests that it is surprisingly tight in high dimensions—at least for $l=2$.

Proposition 3

For any $\epsilon>0$ and $\delta>0$ there exist $c>0$ and $d\in\mathbb{N}$ such that

[TABLE]

where $p=\frac{m}{n}$, and the support points $\bm{\xi}_{i}$, $i\in I$, are sampled independently from the normal distribution $\mathbb{P}$ with mean $\bm{0}\in\mathbb{R}^{d}$ and covariance matrix $(\sqrt{d-1}+c)^{-2}\mathbb{I}\in\mathbb{S}^{d}$.

Proposition 3 can be paraphrased as follows. Sampling the $\bm{\xi}_{i}$, $i\in I$, independently from the normal distribution $\mathbb{P}$ yields a (random) empirical distribution $\hat{\mathbb{P}}_{n}$ that is feasible and $\delta$-suboptimal in (1) with probability $1-\epsilon$. The intuition behind this result is that, in high dimensions, samples drawn from $\mathbb{P}$ are almost orthogonal and close to the surface of the unit ball with high probability. Indeed, these two properties are shared by the worst case distribution (8) in high dimensions.

Proof

of Proposition 3   Theorem 2.1 and Remark 1 imply that

[TABLE]

From the proof of Theorem 2.2 we further know that (17) can be expressed as a continuous function $f(\mathbf{S})$ of the Gram matrix $\mathbf{S}=[\bm{\xi}_{1},\ldots,\bm{\xi}_{n}]^{\top}[\bm{\xi}_{1},\ldots,\bm{\xi}_{n}]$, that is,

[TABLE]

An elementary calculation shows that $f(\mathbb{I})=\frac{n-m}{n}=1-p$. Thus, by the continuity of $f(\cdot)$, there exists $\eta\in(0,1)$ with $\sqrt{f(\mathbf{S})}\geq\sqrt{1-p}-\delta$ whenever $\|\mathbf{S}-\mathbb{I}\|_{\rm max}\leq\eta$.

We are now ready to construct $\mathbb{P}$. First, select $c>0$ large enough to ensure that

[TABLE]

Then, select $d\in\mathbb{N}$ large enough such that

[TABLE]

where $\mathrm{\Phi}(\cdot)$ denotes the univariate standard normal distribution function. Observe that the distribution $\mathbb{P}$ is completely determined by $c$ and $d$. Next, we find that

[TABLE]

where the first and second inequalities follow from the construction of $\eta$ and $c$, respectively, while the third inequality exploits the union bound. We also have

[TABLE]

where the equality holds due to the independence of the $\bm{\xi}_{i}$, while the first and second inequalities follow from Lemma 2.8 in Hopcroft and Kannan (2012) and the construction of $c$, respectively. By the choice of $d$, we finally obtain

[TABLE]

The equality in (20) holds due to the rotation symmetry of $\mathbb{P}$, which implies that

[TABLE]

The claim then follows by substituting (19) and (20) into (18). ∎

Figure 1 compares $\overline{C}_{2}(m,n)$ with $C_{2}(\hat{\mathbb{P}}_{n},m)$ for $n=100$, $m\in\{10,\ldots,100\}$ and $d\in\{25,50,75,100\}$. The $n$ original support points are sampled randomly from the uniform distribution on the unit ball (left panel) and the normal distribution from Proposition 3 with $c=2.97$, which ensures that $\|\bm{\xi}_{i}\|_{2}\leq 1$ with 95% probability (right panel). Note that $C_{2}(\hat{\mathbb{P}}_{n},m)$ is random. Thus, all shown values are averaged across 100 independent trials. Figure 1 confirms that $C_{2}(\hat{\mathbb{P}}_{n},m)$ approaches the worst-case bound $\overline{C}_{2}(n,m)$ as the dimension $d$ increases.

3 Guarantees for Discrete Scenario Reduction

For $n$-point empirical distributions $\hat{\mathbb{P}}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{\bm{\xi}_{i}}$ supported on $\mathbb{R}^{d}$, we now study the loss of optimality incurred by solving the discrete scenario reduction problem instead of its continuous counterpart. More precisely, we want to determine the point-wise largest lower bound $\underline{\kappa}_{l}(n,m)$ and the point-wise smallest upper bound $\overline{\kappa}_{l}(n,m)$ that satisfy

[TABLE]

for the Wasserstein distances of type $l\in\{1,2\}$. Note that the existence of finite bounds $\underline{\kappa}_{l}(n,m)$ and $\overline{\kappa}_{l}(n,m)$ is not a priori obvious as they do not depend on the dimension $d$. Moreover, while it is clear that $\underline{\kappa}_{l}(n,m)\geq 1$ if it exists, it does not seem easy to derive a naïve upper bound on $\overline{\kappa}_{l}(n,m)$.

Our derivations in this section will use the following result, which is the analogue of Theorem 2.1 for the discrete scenario reduction problem.

Theorem 3.1

For any type-$l$ Wasserstein distance induced by any norm $\|\cdot\|$, the discrete scenario reduction problem can be reformulated as

[TABLE]

Proof

The proof is similar to the proof of Theorem 2.1 and is therefore omitted. ∎

The remainder of this section derives lower and upper bounds on $\underline{\kappa}_{l}(n,m)$ and $\overline{\kappa}_{l}(n,m)$ for Wasserstein distances of type $l=2$ (Section 3.1) and $l=1$ (Section 3.2), respectively. To eliminate trivial cases, we assume throughout this section that $n\geq 2$, $m\in\{1,\ldots,n-1\}$ and $d\geq 2$.

3.1 Guarantees for the Type-2 Wasserstein Distance

We first bound $\overline{\kappa}_{2}(n,m)$ in equation (21) from above (Theorem 3.2) and below (Proposition 4).

Theorem 3.2

The upper bound $\overline{\kappa}_{2}(n,m)$ in (21) satisfies $\overline{\kappa}_{2}(n,m)\leq\sqrt{2}$ for all $n,m$.

Proof

The proof proceeds in two steps. We first show that $\overline{\kappa}_{2}(n,m)\leq\sqrt{2}$ for all $n$ when $m=1$ (Step 1). Then we extend the result to all $n$ and $m$ (Step 2).

Step 1:

Fix any $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$. W.l.o.g., we can assume that $\text{mean}(I)=\bm{0}$ and $\frac{1}{n}\sum_{i\in I}\|\bm{\xi}_{i}\|_{2}^{2}=1$ by re-positioning and scaling the atoms $\bm{\xi}_{i}$ appropriately. Note that the re-positioning does not affect $C_{2}(\hat{\mathbb{P}}_{n},1)$ or $D_{2}(\hat{\mathbb{P}}_{n},1)$, and the positive homogeneity of the Wasserstein distance implies that the scaling affects both $C_{2}(\hat{\mathbb{P}}_{n},1)$ and $D_{2}(\hat{\mathbb{P}}_{n},1)$ in the same way and thus preserves their ratio $\overline{\kappa}_{2}(n,1)$. Theorem 2.1 and Remark 1 then imply that

[TABLE]

Step 1 is thus complete if we can show that $D_{2}(\hat{\mathbb{P}}_{n},1)\leq\sqrt{2}$. Indeed, we have

[TABLE]

where the first equality is due to Theorem 3.1, the fourth follows from $\sum_{i\in I}\bm{\xi}_{i}=n\cdot\text{mean}(I)=\bm{0}$, and the inequality holds since $\min_{j\in I}\|\bm{\xi}_{j}\|_{2}^{2}\leq\frac{1}{n}\sum_{i\in I}\|\bm{\xi}_{i}\|_{2}^{2}=1$.

Step 2:

Fix any $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$. Theorem 2.1 and Remark 1 imply that

[TABLE]

For an optimal partition $\{I^{\star}_{j}\}$ to this problem, $C_{2}(\hat{\mathbb{P}}_{n},m)$ can be expressed as

[TABLE]

From our discussion in Step 1 we know that $C_{2,j}$ represents the type-$2$ Wasserstein distance between the conditional empirical distribution $\hat{\mathbb{P}}_{n}^{\,j}=\frac{1}{|I^{\star}_{j}|}\sum_{i\in I^{\star}_{j}}\delta_{\bm{\xi}_{i}}$ and its closest Dirac distribution, that is, $C_{2}(\hat{\mathbb{P}}_{n}^{\,j},1)$. Analogously, we obtain that

[TABLE]

where the first inequality holds since the optimal partition $\{I^{\star}_{j}\}$ for $C_{2}(\hat{\mathbb{P}}_{n},m)$ is typically suboptimal in $D_{2}(\hat{\mathbb{P}}_{n},m)$, the second inequality follows from the fact that $D_{2,j}=D_{2}(\hat{\mathbb{P}}_{n}^{\,j},1)$ and $D_{2}(\hat{\mathbb{P}}_{n}^{\,j},1)\leq\sqrt{2}C_{2}(\hat{\mathbb{P}}_{n}^{\,j},1)$ due to Step 1, and the identity follows from the definition of $C_{2,j}$. The statement now follows. ∎

Proposition 4

There is $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ with $D_{2}(\hat{\mathbb{P}}_{n},m)=\sqrt{2}C_{2}(\hat{\mathbb{P}}_{n},m)$ for all $n,m$.

Proof

In analogy to the proof of Theorem 3.2, we first show the statement for $m=1$ (Step 1) and then extend the result to $m>1$ (Step 2).

Step 1:

The first step in the proof of Theorem 3.2 shows that $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ satisfies $D_{2}(\hat{\mathbb{P}}_{n},1)=\sqrt{2}C_{2}(\hat{\mathbb{P}}_{n},1)$ if $\sum_{i\in I}\bm{\xi}_{i}=\bm{0}$ and $\|\bm{\xi}_{1}\|=\ldots=\|\bm{\xi}_{n}\|=1$. For an even number $n=2k$, $k\in\mathbb{N}$, both conditions are satisfied if we place $\bm{\xi}_{1},\ldots,\bm{\xi}_{k}$ on the surface of the unit ball in $\mathbb{R}^{d}$ and then choose $\bm{\xi}_{k+i}=-\bm{\xi}_{i}$ for $i=1,\ldots,k$ (see left panel of Figure 2 for an illustration in $\mathbb{R}^{3}$). Likewise, for an odd number $n=2k+3$, $k\in\mathbb{N}_{0}$, we can place $\bm{\xi}_{1},\ldots,\bm{\xi}_{k}$ on the surface of the unit ball, choose $\bm{\xi}_{k+i}=-\bm{\xi}_{i}$ for $i=1,\ldots,k$ and fix $\bm{\xi}_{2k+1}=\mathbf{e}_{1}$, $\bm{\xi}_{2k+2}=-\frac{1}{2}\mathbf{e}_{1}+\frac{\sqrt{3}}{2}\mathbf{e}_{2}$ and $\bm{\xi}_{2k+3}=-\frac{1}{2}\mathbf{e}_{1}-\frac{\sqrt{3}}{2}\mathbf{e}_{2}$.

Step 2:

To prove the statement for $m>1$, we construct an empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ whose atoms satisfy $\text{supp}(\hat{\mathbb{P}}_{n})=\mathrm{\Xi}_{1}\cup\mathrm{\Xi}_{2}$ with $|\mathrm{\Xi}_{1}|=n-m+1$ and $|\mathrm{\Xi}_{2}|=m-1$. The atoms $\bm{\xi}_{1},\ldots,\bm{\xi}_{n-m+1}$ in $\mathrm{\Xi}_{1}$ are selected according to the recipe outlined in Step 1, whereas the atoms $\bm{\xi}_{n-m+2},\ldots,\bm{\xi}_{n}$ in $\mathrm{\Xi}_{2}$ satisfy $\bm{\xi}_{n-m+1+i}=(1+iM)\mathbf{e}_{1}$, $i=1,\ldots,m-1$, for any number $M$ satisfying $M>2\,\sqrt{n-m+1}$. A direct calculation then shows that the atoms in $\mathrm{\Xi}_{2}$ are sufficiently far away from those in $\mathrm{\Xi}_{1}$ as well as from each other so that any optimal partition $\{I_{j}^{\star}\}$ to the discrete scenario reduction problem in Theorem 3.1 as well as the continuous scenario reduction problem in Theorem 2.1 consists of the sets $\{i:\bm{\xi}_{i}\in\mathrm{\Xi}_{1}\}$ and $\{i\}$, $\bm{\xi}_{i}\in\mathrm{\Xi}_{2}$. The result then follows from the fact that either problem accumulates a Wasserstein distance of [math] over the atoms in $\mathrm{\Xi}_{2}$, whereas the Wasserstein distance of $D_{2}(\hat{\mathbb{P}}_{n},m)$ is a factor of $\sqrt{2}$ bigger than the Wasserstein distance of $C_{2}(\hat{\mathbb{P}}_{n},m)$ over the atoms in $\mathrm{\Xi}_{1}$ (see Step 1). ∎

Theorem 3.2 and Proposition 4 imply that $\overline{\kappa}_{2}(n,m)=\sqrt{2}$ for all $n$ and $m$, that is, the bound is indeed independent of both the number of atoms $n$ in the empirical distribution and the number of atoms $m$ in the reduced distribution. We now show that the naïve lower bound of $1$ on the approximation ratio is essentially tight.

Proposition 5

The lower bound $\underline{\kappa}_{2}(n,m)$ in (21) satisfies $\underline{\kappa}_{2}(n,m)=1$ whenever $n\geq 3$ and $m\in\{1,\ldots,n-2\}$, while $\underline{\kappa}_{2}(n,n-1)=\sqrt{2}$ always.

Proof

We first prove $\underline{\kappa}_{2}(n,m)=1$ when $m=1$ and $n\geq 3$ (Step 1) and when $m\in\{2,\ldots,n-2\}$ (Step 2). Then, we show $\underline{\kappa}_{2}(n,n-1)=\sqrt{2}$ (Step 3).

Step 1:

Choose $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ such that the first $n-1$ atoms $\bm{\xi}_{1},\ldots,\bm{\xi}_{n-1}$ are selected according to the recipe outlined in Step 1 in the proof of Proposition 4 and $\bm{\xi}_{n}=\bm{0}$. We thus have $\text{mean}(I)=\bm{0}$, and Theorem 2.1 and Remark 1 imply that the optimal continuous scenario reduction is given by the Dirac distribution $\delta_{\bm{0}}$. Since $\bm{0}\in\text{supp}(\hat{\mathbb{P}}_{n})$, we have $C_{2}(\hat{\mathbb{P}}_{n},1)=D_{2}(\hat{\mathbb{P}}_{n},1)$ and the result follows.

Step 2:

To prove the statement for $m>1$, we proceed as in Step 2 in the proof of Proposition 4. In particular, we construct an empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ whose atoms satisfy $\text{supp}(\hat{\mathbb{P}}_{n})=\mathrm{\Xi}_{1}\cup\mathrm{\Xi}_{2}$ with $|\mathrm{\Xi}_{1}|=n-m+1$ and $|\mathrm{\Xi}_{2}|=m-1$. The atoms $\bm{\xi}_{1},\ldots,\bm{\xi}_{m-m+1}$ in $\mathrm{\Xi}_{1}$ are selected according to the recipe outlined in Step 1 of this proof, whereas the remaining atoms $\bm{\xi}_{n-m+1},\ldots,\bm{\xi}_{n}$ in $\mathrm{\Xi}_{2}$ satisfy $\bm{\xi}_{n-m+1+i}=(1+iM)\mathbf{e}_{1}$, $i=1,\ldots,m-1$, for any $M>2\,\sqrt{n-m+1}$. A similar argument as in the proof of Proposition 4 then shows that $C_{2}(\hat{\mathbb{P}}_{n},m)=D_{2}(\hat{\mathbb{P}}_{n},m)$.

Step 3:

Fix any $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$. W.l.o.g., assume that $\{\bm{\xi}_{n-1},\bm{\xi}_{n}\}$ is the closest pair of atoms in terms of Euclidean distance, and let $d_{\text{min}}=\|\bm{\xi}_{n}-\bm{\xi}_{n-1}\|_{2}$. One readily verifies that the partition $I_{j}^{\star}=\{j\}$, $j=1,\ldots,n-2$, and $I_{n-1}^{\star}=\{n-1,n\}$ optimizes both the discrete scenario reduction problem in Theorem 3.1 as well as the continuous scenario reduction problem in Theorem 2.1. We thus have $C_{2}(\hat{\mathbb{P}}_{n},n-1)=\frac{1}{\sqrt{2n}}d_{\text{min}}$ and $D_{2}(\hat{\mathbb{P}}_{n},n-1)=\frac{1}{\sqrt{n}}d_{\text{min}}$, which concludes the proof. ∎

Hence, for any empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ the type-$2$ Wasserstein distance between the minimizer of the discrete scenario reduction problem and $\hat{\mathbb{P}}_{n}$ exceeds the Wasserstein distance between the minimizer of the continuous scenario reduction problem and $\hat{\mathbb{P}}_{n}$ by up to $41.4\%$, and the bound is attainable for any $n,m$.

3.2 Guarantees for the Type-1 Wasserstein Distance

In analogy to Section 3.1, we first bound $\overline{\kappa}_{1}(n,m)$ from above (Theorem 3.3) and below (Proposition 6). In contrast to the previous section, we consider an arbitrary norm $\|\cdot\|$, and we adapt the definition of the geometric median accordingly.

Theorem 3.3

The upper bound $\overline{\kappa}_{1}(n,m)$ in (21) satisfies $\overline{\kappa}_{1}(n,m)\leq 2$ whenever $m\in\{2,\ldots,n-2\}$ as well as $\overline{\kappa}_{1}(n,1)\leq 2\left(1-\frac{1}{n}\right)$ and $\overline{\kappa}_{1}(n,n-1)\leq 1$.

Proof

We first prove the statement for $m=1$ (Step 1) and then extend the result to $m\in\{2,\ldots,n-2\}$ (Step 2) and $m=n-1$ (Step 3).

Step 1:

Fix any $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$. As in the proof of Theorem 3.2, we can assume that $\text{gmed}(I)=\bm{0}$ and $\frac{1}{n}\sum_{i\in I}\|\bm{\xi}_{i}\|=1$ by re-positioning and scaling the atoms $\bm{\xi}_{i}$ appropriately. Theorem 2.1 and Remark 1 then imply that for $m=1$, we have

[TABLE]

Step 1 is thus complete if we can show that $D_{1}(\hat{\mathbb{P}}_{n},1)\leq 2\left(1-\frac{1}{n}\right)$. Indeed, we have

[TABLE]

where the two inequalities follow from the triangle inequality and the fact that $\min_{j\in I}\|\bm{\xi}_{j}\|\leq\frac{1}{n}\sum_{i\in I}\|\bm{\xi}_{i}\|=1$, respectively.

Step 2:

Fix any $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$. Theorem 2.1 and Remark 1 then imply that

[TABLE]

Let $\{I_{j}^{\star}\}$ be an optimal partition for this problem. The same arguments as in the proof of Theorem 3.2 show that

[TABLE]

and the last expression is equal to $2\,C_{1}(\hat{\mathbb{P}}_{n},m)$ by definition of $C_{1,j}$.

Step 3:

For $n=2$ and $m=n-1=1$, Step 1 shows that $\overline{\kappa}_{1}(2,1)\leq 2\left(1-\frac{1}{2}\right)=1$. For $n>2$ and $m=n-1$, the statement can be derived in the same way as the third step in the proof of Proposition 5. We omit the details for the sake of brevity. ∎

Proposition 6

There is $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ such that $D_{1}(\hat{\mathbb{P}}_{n},m)=2\left(1-\frac{m}{n}\right)C_{1}(\hat{\mathbb{P}}_{n},m)$ under the $1$-norm for all $n$ divisible by $2m$, all $m$ and all $d\geq\frac{n}{2m}$.

Proof

We first prove the statement for $m=1$ (Step 1) and then extend the result to $m>1$ (Step 2). Throughout the proof, we set $k=\frac{n}{2m}$ and consider w.l.o.g. the case where $d=k$.

Step 1:

Fix $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ with the atoms $\bm{\xi}_{i}=+\mathbf{e}_{i}$ as well as $\bm{\xi}_{k+i}=-\mathbf{e}_{i}$, $i=1,\ldots,k$. The symmetric placement of the atoms implies that $\text{gmed}(I)=\bm{0}$ and hence $C_{1}(\hat{\mathbb{P}}_{n},1)=1$. Furthermore, we note that $\|\bm{\xi}_{i}-\bm{\xi}_{j}\|_{1}=2$ for all $i\neq j$, that is, any two atoms are equidistant from another (see right panel of Figure 2 for an illustration in $\mathbb{R}^{3}$). By Theorem 3.1, any $1$-point discrete scenario reduction results in a Wasserstein distance of $2\,\frac{n-1}{n}$ to $\hat{\mathbb{P}}_{n}$.

Step 2:

To prove the statement for $m>1$, we construct an empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ whose atoms satisfy $\text{supp}(\hat{\mathbb{P}}_{n})=\bigcup_{j=1}^{m}(\mathrm{\Xi}_{j}^{+}\cup\mathrm{\Xi}_{j}^{-})$ with $|\mathrm{\Xi}_{j}^{+}|=|\mathrm{\Xi}_{j}^{-}|=k$, $j=1,\ldots,m$. The atoms $\bm{\xi}_{2(j-1)k+1},\ldots,\bm{\xi}_{2(j-1)k+k}$ in $\mathrm{\Xi}_{j}^{+}$ satisfy $\bm{\xi}_{2(j-1)k+i}=+\mathbf{e}_{i}+jM\mathbf{e}_{1}$, $i=1,\ldots,k$, whereas the atoms $\bm{\xi}_{2(j-1)k+k+1},\ldots,\bm{\xi}_{2jk}$ in $\mathrm{\Xi}_{j}^{-}$ satisfy $\bm{\xi}_{2(j-1)k+k+i}=-\mathbf{e}_{i}+jM\mathbf{e}_{1}$, $i=1,\ldots,k$, for any number $M$ satisfying $M>2n+2$. The same arguments as in the proof of Theorem 3.2 show that any optimal partition $\{I_{j}^{\star}\}$ to the discrete scenario reduction problem in Theorem 3.1 as well as the continuous scenario reduction problem in Theorem 2.1 consists of the sets indexing the atoms in $\mathrm{\Xi}_{j}^{+}\cup\mathrm{\Xi}_{j}^{-}$, $j=1,\ldots,m$. Step 1 shows that the continuous scenario reduction problem accumulates a Wasserstein distance of $1$ over each set, whereas the discrete scenario reduction problem accumulates a Wasserstein distance of $2\,\frac{2k-1}{2k}$ over each set. The result then follows from the fact that there are $m$ such sets and hence the ratio of the respective overall Wasserstein distances amounts to m\big{(}2\frac{2k-1}{2k}\big{)}/m=2\left(1-\frac{m}{n}\right). ∎

Theorem 3.3 and Proposition 6 imply that $\overline{\kappa}_{1}(n,m)\in[2(1-m/n),2]$ for all $n$ and $m\in\{2,\ldots,n-2\}$. For the small ratios $m:n$ commonly used in practice, we thus conclude that the bound is essentially independent of both the number of atoms $n$ in the empirical distribution and the number of atoms $m$ in the reduced distribution. We close with an analysis of the lower bound $\underline{\kappa}_{1}(n,m)$.

Proposition 7

The lower bound $\underline{\kappa}_{1}(n,m)$ in (21) satisfies $\underline{\kappa}_{1}(n,m)=1$ for all $n,m$.

Proof

The proof widely parallels that of Proposition 5, with the difference that the atoms $\bm{\xi}_{1},\ldots,\bm{\xi}_{n}$ of the empirical distribution $\hat{\mathbb{P}}_{n}$ are placed such that a geometric median (as opposed to the mean) of each subset in the optimal partition coincides with one of the atoms in that subset. This allows both continuous and discrete scenario reduction to choose the same support points for the reduced distribution, hence incurring the same Wasserstein distance. Details are omitted for brevity. ∎

In conclusion, for any $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ the type-$1$ Wasserstein distance between the minimizer of the discrete scenario reduction problem and $\hat{\mathbb{P}}_{n}$ exceeds the Wasserstein distance between the minimizer of the continuous scenario reduction problem and $\hat{\mathbb{P}}_{n}$ by up to $100\%$, and this bound is asymptotically attained for decreasing ratios $m:n$.

4 Solution Methods

We now review existing and propose new solution schemes for the discrete and continuous scenario reduction problems. More precisely, we will study two heuristics for discrete and continuous scenario reduction, respectively, that do not come with approximation guarantees (Section 4.1), we will propose a constant-factor approximation scheme for both the discrete and the continuous scenario reduction problem (Section 4.2), and we will discuss two exact reformulations of these problems as mixed-integer optimization problems (Section 4.3).

In the remainder of this section, we denote by $D_{l}(\hat{\mathbb{P}}_{n},\mathrm{\Xi})$ the type-$l$ Wasserstein distance between $\hat{\mathbb{P}}_{n}$ and its closest distribution supported on the finite set $\mathrm{\Xi}$. Moreover, for an algorithm providing an upper bound $\overline{D}_{l}(\hat{\mathbb{P}}_{n},m)$ on the discrete scenario reduction problem in $\mathbb{R}^{d}$, we define the algorithm’s approximation ratio as the maximum fraction $\overline{D}_{l}(\hat{\mathbb{P}}_{n},m)/D_{l}(\hat{\mathbb{P}}_{n},m)$, where the maximum is taken over all $n$ and $m$, as well as all empirical distributions $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$.

4.1 Heuristics for the Discrete Scenario Reduction Problem

We review in Section 4.1.1 a popular heuristic for the discrete scenario reduction problem due to Dupačová et al (2003). We will show that despite the simplicity and efficiency of the algorithm, there is no finite upper bound on the algorithm’s approximation ratio. In Section 4.1.2 we adapt a widely used clustering heuristic to the continuous scenario reduction problem, and we show that this algorithm’s approximation ratio cannot be bounded from above either.

4.1.1 Dupačová et al.’s Algorithm

We outline Dupačová et al.’s algorithm for the problem $D_{l}(\hat{\mathbb{P}}_{n},m)$ below.

Dupačová et al.’s algorithm for $D_{l}(\hat{\mathbb{P}}_{n},m)$:

Initialize the set of atoms in the reduced set as $R\leftarrow\emptyset$.

Select the next atom to be added to the reduced set as

$\bm{\zeta}\in\underset{\bm{\zeta}\in\text{supp}(\hat{\mathbb{P}}_{n})}{\mathop{\text{arg\,min}}}\ D_{l}(\hat{\mathbb{P}}_{n},R\cup\{\bm{\zeta}\})$

and update $R\leftarrow R\cup\{\bm{\zeta}\}$.

Repeat Step 2 until $|R|=m$.

Given an empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$, the algorithm iteratively populates the reduced set $R$ containing the atoms of the reduced distribution $\mathbb{Q}$. Each atom $\bm{\zeta}\in\text{supp}(\hat{\mathbb{P}}_{n})$ is selected greedily so as to minimize the Wasserstein distance between $\hat{\mathbb{P}}_{n}$ and the closest distribution supported on the augmented reduced set $R\cup\{\bm{\zeta}\}$. After termination, the distribution $\mathbb{Q}$ can be recovered from the reduced set $R$ as follows. Let $\{I_{\bm{\zeta}}\}\in\mathfrak{P}(I,m)$ be any partition of $\text{supp}(\hat{\mathbb{P}}_{n})$ into sets $I_{\bm{\zeta}}$, $\bm{\zeta}\in R$, such that $I_{\bm{\zeta}}$ contains all elements of $\text{supp}(\hat{\mathbb{P}}_{n})$ that are closest to $\bm{\zeta}$ (ties may be broken arbitrarily). Then $\mathbb{Q}=\sum_{\bm{\zeta}\in R}q_{\bm{\zeta}}\delta_{\bm{\zeta}}$, where $q_{\bm{\zeta}}=|I_{\bm{\zeta}}|/n$.

Theorem 4.1

For every $d\geq 2$ and $l,p\geq 1$, the approximation ratio Dupačová et al.’s algorithm is unbounded.

Proof

The proof constructs a specific distribution $\hat{\mathbb{P}}_{n}$ (Step 1), bounds $D_{l}(\hat{\mathbb{P}}_{n},m)$ from above (Step 2) and bounds the Wasserstein distance between $\hat{\mathbb{P}}_{n}$ and the output $\mathbb{Q}$ of Dupačová et al.’s algorithm from below (Step 3).

Step 1:

Fix $d\geq 2$, $l,p\geq 1$ and $m=4$, and consider the empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$ with $n=4z+1$ for some positive integer $z$ as well as $\text{supp}(\hat{\mathbb{P}}_{n})=\mathrm{\Xi}_{1}\cup\cdots\cup\mathrm{\Xi}_{4}\cup\{\bm{\xi}_{4z+1}\}$ with $\mathrm{\Xi}_{j}=\{\bm{\xi}_{(j-1)z+1},\ldots,\bm{\xi}_{jz}\}$, $j=1,\ldots,4$, and

[TABLE]

and $\bm{\xi}_{4z+1}=\bm{0}$, where $\mathcal{B}_{\epsilon}(\bm{x})=\{\bm{\xi}\in\mathbb{R}^{d}\,:\,\|\bm{\xi}-\bm{x}\|_{p}\leq\epsilon\}$ denotes the $\epsilon$-ball around $\bm{x}$. Here, $\epsilon>0$ is small enough so that each atom in $\mathrm{\Xi}_{i}$ is closer to $\bm{0}$ than to any atom in any of the other sets $\mathrm{\Xi}_{j}$. The triangle inequality then implies that

[TABLE]

Step 2:

By construction, we have that

[TABLE]

where the first inequality holds because $\bm{\xi}_{z},\bm{\xi}_{2z},\bm{\xi}_{3z},\bm{\xi}_{4z}\in\text{supp}(\hat{\mathbb{P}}_{n})$, the second inequality holds since moving the atoms in $\mathrm{\Xi}_{j}$ to $\bm{\xi}_{jz}$, $j=1,\ldots,4$, and $\bm{\xi}_{4z+1}$ to $\bm{\xi}_{z}$ represents a feasible transportation plan, and the third inequality is due to (23) and the triangle inequality.

Step 3:

We first show that for a sufficiently small $\epsilon>0$, Dupačová et al.’s algorithm adds $\bm{\xi}_{4z+1}=\bm{0}$ to the reduced set $R$ in the first iteration. We then show that under this selection, the output $\mathbb{Q}$ of Dupačová et al.’s algorithm can be arbitrarily worse than the bound on $D_{l}(\hat{\mathbb{P}}_{n},4)$ determined in the previous step.

To show the first point, the symmetry inherent in $\text{supp}(\hat{\mathbb{P}}_{n})$ implies that it suffices to show that $d_{l}(\hat{\mathbb{P}}_{n},\delta_{\bm{0}})<d_{l}(\hat{\mathbb{P}}_{n},\delta_{\bm{\xi}_{1}})$. To this end, we note that

[TABLE]

due to equation (23), while at the same time

[TABLE]

As $\epsilon$ tends to [math], we have that

[TABLE]

where the strict inequality is due to $l\geq 1$. As a consequence, we may conclude that there indeed exists an $\epsilon>0$ such that Dupačová et al.’s algorithm adds $\bm{\xi}_{4z+1}=\bm{0}$ to the reduced set $R$ in the first iteration.

As for the second point, we note that after adding $\bm{\xi}_{4z+1}=\bm{0}$ to the reduced set $R$, there must be at least one subset $\mathrm{\Xi}_{j}$, $j\in\{1,\ldots,4\}$, such that no $\bm{\xi}_{i}\in\mathrm{\Xi}_{j}$ is contained in the final reduced set $R$. Assume w.l.o.g. that this is the case for $j=1$. We then have

[TABLE]

and combining this with the result of Step 2, we can conclude that the approximation ratio $d_{l}(\hat{\mathbb{P}}_{n},\mathbb{Q})/D_{l}(\hat{\mathbb{P}}_{n},4)$ approaches $\infty$ as $z\to\infty$ and $z\epsilon^{l}\to 0$. ∎

We remark that the algorithm of Dupačová et al. can be improved by adding multiple atoms to the reduced set $R$ in Step 2. Nevertheless, a similar argument as in the proof of Theorem 4.1 shows that the resulting improved algorithm does not allow for a finite upper bound on the approximation ratio either.

4.1.2 $k$-Means Clustering Algorithm

The $k$-means clustering algorithm has first been proposed in 1957 for a pulse-code modulation problem (Lloyd 1982), and it has since then become a widely used heuristic for various classes of clustering problems. It aims to partition a set of observations $\bm{x}_{1},\ldots,\bm{x}_{n}$ into $m$ clusters $S_{1},\ldots,S_{m}$ such that the intra-cluster sums of squared distances are minimized. By generalizing the algorithm to arbitrary powers and norms, we can adapt the algorithm to our continuous scenario reduction problem as follows.

$k$-means clustering algorithm for $C_{l}(\hat{\mathbb{P}}_{n},m)$:

Initialize the reduced set $R=\{\bm{\zeta}_{1},\ldots,\bm{\zeta}_{m}\}\subseteq\text{supp}(\hat{\mathbb{P}}_{n})$ arbitrarily.

Let $\{I_{j}\}\in\mathfrak{P}(I,m)$ be any partition whose sets $I_{j}$, $j\in J$, contain all atoms of $\text{supp}(\hat{\mathbb{P}}_{n})$ that are closest to $\bm{\zeta}_{j}$ (ties may be broken arbitrarily).

For each $j\in J$, update $\bm{\zeta}_{j}\leftarrow\mathop{\text{arg\,min}}\{\sum_{i\in I_{j}}\|\bm{\xi}_{i}-\bm{\zeta}\|^{l}\,:\,\bm{\zeta}\in\mathbb{R}^{d}\}$.

Repeat Steps 2 and 3 until the reduced set $R$ no longer changes.

For the empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$, the algorithm iteratively updates the reduced set $R$ containing the atoms of the reduced distribution $\mathbb{Q}$ through a sequence of assignment (Step 2) and update (Step 3) steps. Step 2 assigns each atom $\bm{\xi}_{i}\in\text{supp}(\hat{\mathbb{P}}_{n})$ of the empirical distribution to the closest atom in the reduced set, and Step 3 updates each atom in the reduced set so as to minimize the sum of $l$-th powers of the distances to its assigned atoms from $\text{supp}(\hat{\mathbb{P}}_{n})$. After termination, the continuously reduced distribution $\mathbb{Q}$ can be recovered from the reduced set $R$ in the same way as in the previous subsection.

Remark 1 implies that for $l=2$ and $\|\cdot\|=\|\cdot\|_{2}$, Step 3 reduces to $\bm{\zeta}_{j}\leftarrow\frac{1}{|I_{j}|}\sum_{i\in I_{j}}\bm{\xi}_{i}$, in which case we recover the classical $k$-means clustering algorithm. Although the algorithm terminates at a local minimum, Dasgupta (2008) has shown that for $l=2$ and $\|\cdot\|=\|\cdot\|_{2}$, the solution determined by the algorithm can be arbitrarily suboptimal. We now generalize this finding to generic type-$l$ Wasserstein distances induced by arbitrary $p$-norms.

Theorem 4.2

If initialized randomly in Step 1, the approximation ratio of the $k$-means clustering algorithm is unbounded for every $d,l,p\geq 1$ with significant probability.

Proof

In analogy to the proof of Theorem 4.1, we construct a specific distribution $\hat{\mathbb{P}}_{n}$ (Step 1), bound $D_{l}(\hat{\mathbb{P}}_{n},m)$ from above (Step 2) and bound the Wasserstein distance between $\hat{\mathbb{P}}_{n}$ and the output $\mathbb{Q}$ of the $k$-means algorithm from below (Step 3).

Step 1:

Fix $d,l,p\geq 1$ and m = 3, and consider the empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},m)$ with $n=3z+1$ for some positive integer $z$ as well as $\text{supp}(\hat{\mathbb{P}}_{n})=\mathrm{\Xi}_{1}\cup\mathrm{\Xi}_{2}\cup\{\bm{\xi}_{3z+1}\}$ with $\mathrm{\Xi}_{1}=\{\bm{\xi}_{1},\ldots,\bm{\xi}_{2z}\}$, $\mathrm{\Xi}_{2}=\{\bm{\xi}_{2z+1},\ldots,\bm{\xi}_{3z}\}$ and

[TABLE]

as well as $\bm{\xi}_{3z+1}=\mathbf{e}_{1}$, where again $\mathcal{B}_{\epsilon}(\bm{x})=\{\bm{\xi}\in\mathbb{R}^{d}\,:\,\|\bm{\xi}-\bm{x}\|_{p}\leq\epsilon\}$. By construction, the distance between any pair of atoms in $\mathrm{\Xi}_{j}$ is bounded above by $2\epsilon<\frac{1}{2}$, $j=1,2$, whereas the distance between two atoms from $\mathrm{\Xi}_{1}$ and $\mathrm{\Xi}_{2}$ is bounded from below by $1-2\epsilon>\frac{1}{2}$.

Step 2:

As similar argument as in the proof of Theorem 4.1 shows that

[TABLE]

where the last inequality follows from (24).

Step 3:

We first show that with significant probability, the algorithm chooses a reduced set $R$ containing two atoms from $\mathrm{\Xi}_{1}$ and one atom from $\mathrm{\Xi}_{2}$ in the first step. We then show that under this initialization, the output $\mathbb{Q}$ of the algorithm can be arbitrarily worse than the bound on $C_{l}(\hat{\mathbb{P}}_{n},3)$ determined above.

In view of the first point, we note that the probability of the reduced set $R$ containing two atoms from $\mathrm{\Xi}_{1}$ and one atom from $\mathrm{\Xi}_{2}$ after the first step is $\binom{2z}{2}\binom{z}{1}/\binom{3z+1}{3}$ and approaches 44.44% as $z\to\infty$. In the following, we thus assume w.l.o.g. that $R=\{\bm{\zeta}_{1},\bm{\zeta}_{2},\bm{\zeta}_{3}\}$ with $\bm{\zeta}_{1},\bm{\zeta}_{2}\in\mathrm{\Xi}_{1}$ and $\bm{\zeta}_{3}\in\mathrm{\Xi}_{2}$ after the first step.

As for the second point, we note that Step 2 of the algorithm assigns the atoms $\bm{\xi}_{i}\in\mathrm{\Xi}_{1}$ to either $\bm{\zeta}_{1}$ or $\bm{\zeta}_{2}$, whereas the atoms $\bm{\xi}_{i}\in\mathrm{\Xi}_{2}\cup\{\bm{\xi}_{3z+1}\}$ are assigned to $\bm{\zeta}_{3}$. Hence, the update of the reduced set $R$ in the next iteration satisfies $\bm{\zeta}_{1},\bm{\zeta}_{2}\in\mathcal{B}_{\epsilon}(-\mathbf{e}_{1})$, whereas $\bm{\zeta}_{3}$ is chosen with respect to the set $\mathrm{\Xi}_{2}\cup\{\bm{\xi}_{3z+1}\}$. The algorithm then terminates in the third iteration as the reduced set $R$ no longer changes. We thus find that

[TABLE]

Recall that $\bm{\xi}_{3z}\in\mathcal{B}_{\epsilon}(\bm{0})$ and $\bm{\xi}_{3z+1}=\mathbf{e}_{1}$, which implies that

[TABLE]

and that at least one of the two terms $\|\bm{\xi}_{3z}-\bm{\zeta}_{3}\|_{p}$ or $\|\bm{\xi}_{3z+1}-\bm{\zeta}_{3}\|_{p}$ is greater than $\frac{1-\epsilon}{2}$. We thus conclude that

[TABLE]

which by virtue Step 2 implies that $d_{l}(\hat{\mathbb{P}}_{n},\mathbb{Q})/C_{l}(\hat{\mathbb{P}}_{n},3)\rightarrow\infty$ as $\epsilon\to 0$. ∎

4.2 Constant-Factor Approximation for the Scenario Reduction Problem

We now propose a simple approximation scheme for the discrete scenario reduction problem under the type-$1$ Wasserstein distance whose approximation ratio is bounded from above by $5$. We also show that this algorithm gives rise to an approximation scheme for the continuous scenario reduction problem with an approximation ratio of $10$. To our best knowledge, we describe the first constant-factor approximations for the discrete and continuous scenario reduction problems.

Our algorithm follows from the insight that the discrete scenario reduction problem under the type-$1$ Wasserstein distance is equivalent to the $k$-median clustering problem. The $k$-median clustering problem is a variant of the $k$-means clustering problem described in Section 4.1.2, where the $l$-th power of the norm terms is dropped (i.e., $l=1$). In the following, we adapt a well-known local search algorithm (Arya et al 2004) to our discrete scenario reduction problem:

Local search algorithm for $D_{l}(\hat{\mathbb{P}}_{n},m)$:

Initialize the reduced set $R\subseteq\text{supp}(\hat{\mathbb{P}}_{n})$, $|R|=m$, arbitrarily.

Select the next exchange to be applied to the reduced set as

(\bm{\zeta},\bm{\zeta}^{\prime})\in\mathop{\text{arg\,min}}\left\{D_{l}(\hat{\mathbb{P}}_{n},R\cup\{\bm{\zeta}\}\setminus\{\bm{\zeta}^{\prime}\})\,:\,(\bm{\zeta},\bm{\zeta}^{\prime})\in\big{(}\text{supp}(\hat{\mathbb{P}}_{n})\setminus R\big{)}\times R\right\},

and update $R\leftarrow R\cup\{\bm{\zeta}\}\setminus\{\bm{\zeta}^{\prime}\}$ if $D_{l}(\hat{\mathbb{P}}_{n},R\cup\{\bm{\zeta}\}\setminus\{\bm{\zeta}^{\prime}\})<D_{l}(\hat{\mathbb{P}}_{n},R)$.

Repeat Step 2 until no further improvement is possible.

For an empirical distribution $\hat{\mathbb{P}}_{n}\in\mathcal{P}_{\mathrm{E}}(\mathbb{R}^{d},n)$, the algorithm constructs a sequence of reduced sets $R$ containing the atoms of the reduced distribution $\mathbb{Q}$. In each iteration, Step 2 selects the exchange $R\cup\{\bm{\zeta}\}\setminus\{\bm{\zeta}^{\prime}\}$, $\bm{\zeta}\in\text{supp}(\hat{\mathbb{P}}_{n})$ and $\bm{\zeta}^{\prime}\in R$, that maximally reduces the Wasserstein distance $D_{l}(\hat{\mathbb{P}}_{n},R)$. For performance reasons, this ‘best fit’ strategy can also be replaced with a ‘first fit’ strategy which conducts the first exchange $R\cup\{\bm{\zeta}\}\setminus\{\bm{\zeta}^{\prime}\}$ found that leads to a reduction of $D_{l}(\hat{\mathbb{P}}_{n},R)$. After termination, the reduced distribution $\mathbb{Q}$ can be recovered from the reduced set $R$ in the same way as in Section 4.1.1.

It follows from Arya et al (2004) that the above algorithm (with either ‘best fit’ or ‘first fit’) has an approximation ratio of $5$ for the discrete scenario reduction problem for all $d$. We now show that the algorithm also provides solutions to the continuous scenario reduction problem with an approximation ratio of at most $10$.

Corollary 1

The problems $D_{l}(\hat{\mathbb{P}}_{n},m)$ and $C_{l}(\hat{\mathbb{P}}_{n},m)$ are related as follows.

Any approximation algorithm for $D_{2}(\hat{\mathbb{P}}_{n},m)$ under the $2$-norm with approximation ratio $\alpha$ gives rise to an approximation algorithm for $C_{2}(\hat{\mathbb{P}}_{n},m)$ under the $2$-norm with approximation ratio $\sqrt{2}\alpha$. 2. 2.

Any approximation algorithm for $D_{1}(\hat{\mathbb{P}}_{n},m)$ under any norm with approximation ratio $\alpha$ gives rise to an approximation algorithm for $C_{1}(\hat{\mathbb{P}}_{n},m)$ under the same norm with approximation ratio $2\alpha$.

Proof

The two statements follow directly from Theorems 3.2 and 3.3, respectively. ∎

As presented, the local search algorithm is not guaranteed to terminate in polynomial time. This can be remedied by a variant of the algorithm that only accepts exchanges $R\cup\{\bm{\zeta}\}\setminus\{\bm{\zeta}^{\prime}\}$ that reduce the Wasserstein distance $D_{l}(\hat{\mathbb{P}}_{n},R)$ by at least $\epsilon/((n-m)m)$ for some constant $\epsilon>0$. It follows from Arya et al (2004) that for any $\epsilon$, this variant terminates in polynomial time and provides a $(5+\epsilon)$-approximation for the discrete scenario reduction problem. The algorithm can also be extended to accommodate multiple swaps in every iteration, which lowers the approximation ratio to $3+\epsilon$ at the expense of additional computations.

We remark that there is a wealth of algorithms for the $k$-median problem that can be adapted to the discrete scenario reduction problem. For example, Charikar and Li (2012) present a rounding scheme for the $k$-median problem that gives rise to a polynomial-time algorithm for $D_{1}(\hat{\mathbb{P}}_{n},R)$ with an approximation ratio of $3.25$. Likewise, Kanungo et al (2004) propose a local search algorithm for the $k$-median problem that gives rise to a polynomial-time algorithm for $D_{2}(\hat{\mathbb{P}}_{n},R)$ under the $2$-norm with an approximation ratio of $9+\epsilon$. In both cases, Corollary 1 allows us to extend these guarantees to the corresponding versions of the continuous scenario reduction problem.

4.3 Mixed-Integer Reformulations of the Discrete and Continuous Scenario Reduction Problems

We first review a well-known mixed-integer linear programming (MILP) reformulation of the discrete scenario reduction problem $D_{l}(\hat{\mathbb{P}}_{n},m)$:

Theorem 4.3

The discrete scenario reduction problem can be formulated as the MILP

[TABLE]

with $\mathbf{D}\in\mathbb{S}^{n}$ and $d_{ij}=\|\bm{\xi}_{i}-\bm{\xi}_{j}\|^{l}$ encoding the distances among the atoms in $\text{supp}(\hat{\mathbb{P}}_{n})$.

Proof

See e.g. Heitsch and Römisch (2003). ∎

In problem (25), the decision variable $\pi_{ij}$ determines how much of the probability mass of atom $\bm{\xi}_{i}$ in $\hat{\mathbb{P}}_{n}$ is shifted to the atom $\bm{\zeta}_{j}$ in the reduced distribution $\mathbb{Q}$, whereas the decision variable $\lambda_{j}$ determines whether the atom $\bm{\xi}_{j}\in\text{supp}(\hat{\mathbb{P}}_{n})$ is contained in the support of $\mathbb{Q}$. A solution $(\mathbf{\Pi}^{\star},\bm{\lambda}^{\star})$ to problem (25) allows us to recover the reduced distribution via $\mathbb{Q}=\frac{1}{n}\sum_{j=1}^{n}\mathbf{e}^{\top}\mathbf{\Pi}^{\star}\mathbf{e}_{j}\cdot\delta_{\bm{\xi}_{j}}$. Problem (25) has $n$ binary and $n^{2}$ continuous variables as well as $n^{2}+n+1$ constraints.

We now consider the continuous scenario reduction problem. Due to its bilinear objective function, which involves products of transportation weights $\pi_{ij}$ and the distances $\|\bm{\xi}_{i}-\bm{\zeta}_{j}\|^{l}$ containing the continuous decision variables $\bm{\zeta}_{j}$, this problem may not appear to be amenable to a reformulation as a mixed-integer convex optimization problem. We now show that such a reformulation indeed exists.

Theorem 4.4

The continuous scenario reduction problem can be formulated as the mixed-integer convex optimization problem

[TABLE]

where $M=\max_{i,j\in I}\|\bm{\xi}_{i}-\bm{\xi}_{j}\|^{l}$ denotes the diameter of the support of $\hat{\mathbb{P}}_{n}$.

In problem (26), the decision variable $\pi_{ij}$ determines whether or not the probability mass of atom $\bm{\xi}_{i}$ in the empirical distribution $\hat{\mathbb{P}}_{n}$ is shifted to the atom $\bm{\zeta}_{j}$ in the reduced distribution $\mathbb{Q}$, whereas the decision variable $c_{i}$ records the cost of moving the atom $\bm{\xi}_{i}$ under the transportation plan $\mathbf{\Pi}$.

Proof

of Theorem 4.4   We prove the statement by showing that optimal solutions to problem (26) correspond to feasible solutions in problem (2) with the same objective function value in their respective problems and vice versa.

Fix a minimizer $(\mathbf{\Pi}^{\star},\bm{c}^{\star},\bm{\zeta}_{1}^{\star},\ldots,\bm{\zeta}_{m}^{\star})$ to problem (26), which corresponds to a feasible solution $(\{I_{j}\},\bm{\zeta}_{1}^{\star},\ldots,\bm{\zeta}_{m}^{\star})$ in problem (2) if we set $I_{j}=\{i\in I\,:\,\pi_{ij}^{\star}=1\}$ for all $j\in J$. Note that $c_{i}^{\star}=\|\bm{\xi}_{i}-\bm{\zeta}_{j}^{\star}\|^{l}$ for $j\in J$ and $i\in I_{j}$. Thus, both solutions adopt the same objective value.

Conversely, fix a minimizer $(\{I_{j}^{\star}\},\bm{\zeta}_{1}^{\star},\ldots,\bm{\zeta}_{m}^{\star})$ to problem (2). This solution corresponds to a feasible solution $(\mathbf{\Pi},\bm{c},\bm{\zeta}_{1}^{\star},\ldots,\bm{\zeta}_{m}^{\star})$ to problem (26) if we set $\pi_{ij}=1$ if $i\in I_{j}^{\star}$ and $\pi_{ij}=0$ otherwise for all $j\in J$, as well as $c_{i}=\|\bm{\xi}_{i}-\bm{\zeta}_{j}^{\star}\|^{l}$ for all $i\in I_{j}$ and $j\in J$. By construction, both solutions adopt the same objective value. ∎

A solution $(\mathbf{\Pi}^{\star},\bm{c}^{\star},\bm{\zeta}^{\star}_{1},\ldots,\bm{\zeta}^{\star}_{m})$ to problem (26) allows us to recover the reduced distribution via $\mathbb{Q}=\frac{1}{n}\sum_{j=1}^{m}\mathbf{e}^{\top}\mathbf{\Pi}^{\star}\mathbf{e}_{j}\cdot\delta_{\bm{\zeta}^{\star}_{j}}$. Problem (26) has $nm$ binary and $n+md$ continuous variables as well as $nm+n$ constraints. We now show that (26) typically reduces to an MILP or a mixed-integer second-order cone program (MISOCP).

Proposition 8

For the type-$1$ Wasserstein distance induced by $\|\cdot\|_{1}$ or $\|\cdot\|_{\infty}$, problem (26) reduces to an MILP. For any type-$l$ Wasserstein distance induced by $\|\cdot\|_{p}$, where $l\geq 1$ and $p\geq 1$ are rational numbers, problem (26) reduces to an MISOCP.

Proof

In view of the first statement, we note that $\|\bm{\xi}_{i}-\bm{\zeta}_{j}\|_{1}\leq c_{i}+M(1-\pi_{ij})$ is satisfied if and only if there is $\bm{\phi}_{ij}\in\mathbb{R}^{d}$ such that

[TABLE]

Likewise, $\|\bm{\xi}_{i}-\bm{\zeta}_{j}\|_{\infty}\leq c_{i}+M(1-\pi_{ij})$ holds if and only if there is $\phi_{ij}\in\mathbb{R}$ with

[TABLE]

As for the second statement, we note that $\|\bm{\xi}_{i}-\bm{\zeta}_{j}\|_{p}^{l}\leq c_{i}+M(1-\pi_{ij})$ is satisfied if and only if there is $\phi_{ij}\in\mathbb{R}$ such that

[TABLE]

For rational $l,p\geq 1$, both inequalities can be expressed through finitely many second-order cone constraints (Alizadeh and Goldfarb 2003, Section 2.3). ∎

5 Numerical Experiment: Color Quantization

Color quantization aims to reduce the color palette of a digital image without compromising its visual appearance. In the standard RGB24 model colors are encoded by vectors of the form $(r,g,b)\in\{0,1,\ldots,255\}^{3}$. This means that the RGB24 model can represent a vast number of 16,777,216 distinct colors. Consequently, color quantization serves primarily as a lossy image compression method.

In the following we interpret the color quantization problem as a discrete scenario reduction problem using the type-1 Wasserstein distance induced by the 1-norm on $\mathbb{R}^{3}$. Thus, we can solve color quantization problems via Dupačová’s greedy heuristic, the local search algorithm or the exact MILP reformulation (25). In our experiment we aim to compress all 24 pictures from the Kodak Lossless True Color Image Suite (http://r0k.us/graphics/kodak/) to $m=2^{1},\ldots,2^{9}$ colors. As the MILP reformulation scales poorly with $n$, we first reduce each image to $n\,\lesssim\,$1,024 colors using the Linux command “convert -colors”, which is distributed through ImageMagick (https://www.imagemagick.org). We henceforth refer to the resulting 1,024-color images as the originals.

In all experiments we use an efficient variant of Dupačová’s algorithm due to Heitsch and Römisch (2003) (DPCV), and we initialize the local search algorithm either with the color palette obtained from Dupačová’s algorithm (LOC-1) or naïvely with the $m$ most frequent colors of the original image (LOC-2). The MILP (25) is solved with GUROBI 7.0.1 (MILP). All algorithms are implemented in C++, and all experiments are executed on a 3.40GHz i7 CPU machine with 16GB RAM. We report the average and the worst-case runtimes in Table 1. Note that DPCV, LOC-1 and LOC-2 all terminate in less than 14 seconds across all instances, while MILP requires substantially more time (the maximum runtime was set to ten hours). Moreover, warmstarting the local search algorithm with the color palette obtained from DPCV can significantly reduce the runtimes.

As an example, Figure 3 shows the image “kodim15.png” as well as the results of different color quantization algorithms for $m=16$. While the outputs of LOC-1, LOC-2 and MILP are almost indistinguishable, the output of DPCV has ostensible deficiencies (e.g., it misrepresents the yellow color around the subject’s eye). For comparison, we also show the output of the color quantization routine in Microsoft Paint (MS Paint). Figure 4 visualizes the optimality gaps of DPCV, LOC-1 and LOC-2 relative to MILP (i.e., their respective approximation ratio $-$ $1$). Our experiment suggests that the local search algorithm is competitive with MILP in terms of output quality but at significantly reduced runtimes. Moreover, the local search algorithm LOC-1 warmstarted with the color palette obtained from DPCV is guaranteed to outperform DPCV in terms of optimality gaps.

Acknowledgements

This research was funded by the SNSF grant BSCGI0_157733 and the EPSRC grants EP/M028240/1 and EP/M027856/1.

Appendix: Auxiliary Results

The proofs of Theorem 2.2 relies on the following two lemmas.

Lemma 1

The semidefinite program (6) admits an optimal solution $(\tau,\mathbf{S})$ with $\mathbf{S}=\alpha\mathbb{I}+\beta\bm{11}^{\top}$ for some $\alpha,\beta\in\mathbb{R}$.

Proof

Let $(\tau,\mathbf{S}^{\star})$ be any optimal solution to (6), which exists because (6) has a continuous objective function and a compact feasible set, and denote by $\mathfrak{S}$ the set of all permutations of $I$. For any $\sigma\in\mathfrak{S}$, the permuted solution $(\tau,\mathbf{S}^{\sigma})$, with $s^{\sigma}_{ij}=s^{\star}_{\sigma(i)\sigma(j)}$ is also optimal in (6). Note first that $(\tau,\mathbf{S}^{\sigma})$ is feasible in (6) because

[TABLE]

where the index sets $I^{\sigma}_{j}=\{\sigma(i):~{}i\in I_{j}\}$ for $j\in J$ form an $m$-set partition from within $\mathfrak{P}(I,m)$, and because $\mathbf{S}^{\sigma}\succeq\bm{0}$ and $s^{\sigma}_{ii}=s^{\star}_{\sigma(i)\sigma(i)}\leq 1$ for all $i\in I$ by construction. Moreover, it is clear that $(\tau,\mathbf{S}^{\sigma})$ and $(\tau,\mathbf{S}^{\star})$ share the same objective value in (6). Thus, $(\tau,\mathbf{S}^{\sigma})$ is optimal in (6) for every $\sigma\in\mathfrak{S}$.

The convexity of problem (6) implies that $(\tau,\mathbf{S})$ with $\mathbf{S}=\frac{1}{n!}\sum_{\sigma\in\mathfrak{S}}\mathbf{S}^{\sigma}$ is also optimal in (6). The claim follows by noting that $\mathbf{S}$ is invariant under permutations of the coordinates and thus representable as $\alpha\mathbb{I}+\beta\bm{11}^{\top}$ for some $\alpha,\beta\in\mathbb{R}$. ∎

Lemma 2

For $\alpha,\beta\in\mathbb{R}$ the eigenvalues of $\mathbf{S}=\alpha\mathbb{I}+\beta\bm{11}^{\top}\in\mathbb{S}^{n}$ are given by $\alpha+n\beta$ (with multiplicity 1) and $\alpha$ (with multiplicity $n-1$).

Proof

Note that $\mathbf{S}$ is a circulant matrix, meaning that each of its rows coincides with the preceding row rotated by one element to the right. Thus, the eigenvalues of $\mathbf{S}$ are given by $\alpha+\beta(1+\rho_{j}^{1}+\ldots\rho_{j}^{n-1})$, $j=0,\ldots,n-1$, where $\rho_{j}=e^{2\pi ij/n}$ and $i$ denotes the imaginary unit; see e.g. Gray (2006). For $j=0$ we then obtain the eigenvalue $\alpha+n\beta$, and for $j=1,\ldots,n-1$ we obtain the other $n-1$ eigenvalues, all of which equal $\alpha$ because $\sum_{k=0}^{n-1}e^{2\pi ijk/n}=(1-e^{2\pi ij})/(1-e^{2\pi ij/n})=0$. ∎