An LP-based, Strongly-Polynomial 2-Approximation Algorithm for Sparse Wasserstein Barycenters

TL;DR

This paper introduces a strongly-polynomial 2-approximation algorithm for computing sparse Wasserstein barycenters, significantly improving efficiency over exact linear programming methods while maintaining low approximation error.

Contribution

It presents a novel strongly-polynomial 2-approximation algorithm for sparse Wasserstein barycenters that combines support union computation with support refinement for non-mass splitting transports.

Findings

Algorithm runs in strongly-polynomial time

Achieves a 2-approximation with sparse support

Demonstrates significant speedup over exact LP methods

Abstract

Discrete Wasserstein barycenters correspond to optimal solutions of transportation problems for a set of probability measures with finite support. Discrete barycenters are measures with finite support themselves and exhibit two favorable properties: there always exists one with a provably sparse support, and any optimal transport to the input measures is non-mass splitting. It is open whether a discrete barycenter can be computed in polynomial time. It is possible to find an exact barycenter through linear programming, but these programs may scale exponentially. In this paper, we prove that there is a strongly-polynomial 2-approximation algorithm based on linear programming. First, we show that an exact computation over the union of supports of the input measures gives a tight 2-approximation. This computation can be done through a linear program with setup and solution in strongly-polynomial time. The resulting measure is sparse, but an optimal transport may split mass. We then devise a second, strongly-polynomial algorithm to improve this measure to one with a non-mass splitting transport of lower cost. The key step is an update of the possible support set to resolve mass split. Finally, we devise an iterative scheme that alternates between these two algorithms. The algorithm terminates with a 2-approximation that has both a sparse support and an associated non-mass splitting optimal transport. We conclude with some sample computations and an analysis of the scaling of our algorithms, exhibiting vast improvements in running time over exact LP-based computations and low practical errors.

Tables

Table 1. Table 1: Average numbers (error, time) for an initial approximation P ¯ org subscript ¯ 𝑃 org \bar{P}_{\text{org}} , first P ¯ ′ superscript ¯ 𝑃 ′ \bar{P}^{\prime} , full runs of Algorithm 3 , and an exact computation in a grid setting. The numbers in each row were derived from 100 100 100 random samples of 16 × 16 16 16 16\times 16 digits from the MNIST data set.

	first ${\overline{P}}_{\text{org}}$		first ${\overline{P}}^{\prime }$	full run of Alg. 3			exact
no. of measures	error	time (s)	error	error	time (s)	iterations	time (s)
4	14.8%	4.2	3.8%	3.1%	9.9	2.2	120
5	15.2%	5.7	4.4%	4.1%	16.4	2.8	204
6	15.5%	8.5	4.5%	3.9%	22.3	2.7	540
7	15.1%	12.1	4.6%	4.2%	29.8	3.1	1602
8	16.2%	16.3	5.2%	4.8%	36.7	3.1	4330
9	n/a	23.0	n/a	n/a	45.2	3.4	–
12	n/a	39.1	n/a	n/a	74.8	3.3	–
16	n/a	58.4	n/a	n/a	99.3	3.7	–
20	n/a	90.3	n/a	n/a	169.2	4.5	–
40	n/a	298.7	n/a	n/a	557.3	5.0	–
70	n/a	681.5	n/a	–	–	–	–
100	n/a	1198.2	n/a	–	–	–	–

Table 2. Table 2: A comparison (error, time) of an initial approximation P ¯ org subscript ¯ 𝑃 org \bar{P}_{\text{org}} and first P ¯ ′ superscript ¯ 𝑃 ′ \bar{P}^{\prime} to the result Q ¯ ¯ 𝑄 \bar{Q} of a widely-used regularization-based approximation algorithm ( Cuturi and Doucet ( 2014 ) ). The error gap shows the ratio ϕ ​ ( Q ¯ ) / ϕ ​ ( P ¯ ) italic-ϕ ¯ 𝑄 italic-ϕ ¯ 𝑃 \phi(\bar{Q})/\phi(\bar{P}) of objective function values for Q ¯ ¯ 𝑄 \bar{Q} and P ′ ¯ ¯ superscript 𝑃 ′ \bar{P^{\prime}} .

	first ${\overline{P}}_{\text{org}}$		first ${\overline{P}}^{\prime }$	regularization-based $\overline{Q}$
no. of measures	error	time (s)	error	error	time (s)	error gap $\varphi \left(\overline{Q}\right)/\varphi \left(\overline{P}\right)$
4	14.8%	4.2	3.8%	22.4%	6.8	1.179
6	15.5%	8.5	4.5%	21.9%	10.2	1.167
8	16.2%	16.3	5.2%	23.6%	13.7	1.175
12	n/a	39.1	n/a	n/a	20.4	1.182
16	n/a	58.4	n/a	n/a	27.7	1.188
20	n/a	90.3	n/a	n/a	33.0	1.179
40	n/a	298.7	n/a	n/a	67.9	1.191
70	n/a	681.5	n/a	n/a	115.5	1.185
100	n/a	1198.2	n/a	n/a	164.9	1.182

Table 3. Table 3: Average numbers (error, time) for an initial approximation P ¯ org subscript ¯ 𝑃 org \bar{P}_{\text{org}} , first P ¯ ′ superscript ¯ 𝑃 ′ \bar{P}^{\prime} , and full runs of Algorithm 3 , and an exact computation for data in general position. The numbers in each row were derived from 100 100 100 runs on a set of measures with 9 9 9 support points of randomly chosen masses.

	first ${\overline{P}}_{\text{org}}$		first ${\overline{P}}^{\prime }$	full run of Alg. 3			exact
no. of measures	error	time (s)	error	error	time (s)	iterations	time (s)
8	10.1%	0.7	2.0%	1.6%	0.9	1.4	505
12	9.8%	1.1	2.4%	1.9%	3.0	1.8	12400
50	n/a	2.3	n/a	n/a	8.0	2.2	–
100	n/a	3.5	n/a	n/a	26.6	2.3	–
200	n/a	5.3	n/a	n/a	112.8	2.9	–
500	n/a	7.9	n/a	n/a	730.4	3.5	–
1000	n/a	18.5	n/a	–	–	–	–
5000	n/a	92.0	n/a	–	–	–	–
10000	n/a	229.2	n/a	–	–	–	–
20000	n/a	963.4	n/a	–	–	–	–

Equations

$$\varphi (\bar{P}):=\sum _{i=1}^N{\lambda }_i{W}_2(\bar{P},P_i{)}^2=\underset{P\in {\mathcal{P}}^2({\mathbb{R}}^d)}{\inf }\sum _{i=1}^N{\lambda }_i{W}_2(P,P_i{)}^2,$$

$$\varphi (\bar{P})=\sum _{i=1}^N{\lambda }_i{W}_2(\bar{P},P_i{)}^2=\underset{P\in {\mathcal{P}}^2({\mathbb{R}}^d)}{\inf }\sum _{i=1}^N{\lambda }_i{W}_2(P,P_i{)}^2.$$

$$S:=\left\{\sum _{i=1}^N{\lambda }_i{x}_i:x_i\in \text{supp}(P_i)\right\}.$$

$$\varphi (\bar{P})=\underset{P\in {\mathcal{P}}_{\mathcal{S}}^2({\mathbb{R}}^d)}{\inf }\sum _{i=1}^N{\lambda }_i{W}_2(P,P_i{)}^2.$$

$$\underset{y}{\min }$$

$$\text{s.t.}\sum _{k=1}^{|P_i|}{y}_{ijk}=$$

$$\sum _{j=1}^{|P_0|}{y}_{ijk}=$$

$$y_{ijk}\ge$$

$$\underset{z,y}{\min }$$

$$\text{s.t.}\sum _{k=1}^{|P_i|}{y}_{ijk}=$$

$$\sum _{j=1}^{|S_0|}{y}_{ijk}=$$

$$y_{ijk}\ge$$

$$z_j\in$$

$$|\bar{P}|\le \sum _{i=1}^N|P_i|-N+1.$$

$$|{\bar{P}}_0|\le \sum _{i=1}^N|P_i|-N+1.$$

$$\sum _{i=1}^N{\lambda }_i\Vert (s+c)-x_i{\Vert }^2=s^{T}s-c^{T}c+\sum _{i=1}^N{\lambda }_i{x}_i^T{x}_i,$$

$$S:=\left\{\sum _{i=1}^N{\lambda }_i{x}_i:x_i\in \text{supp}(P_i)\right\},$$

$$S_{\text{org}}:=\bigcup _{i=1}^N\text{supp}(P_i),$$

$$\varphi ({\bar{P}}_{\text{org}})\le 2\cdot \varphi (\bar{P})$$

$$\sum _{i=1}^N{\lambda }_i\Vert s-x_i{\Vert }^2=s^{T}s-2c^{T}s+\sum _{i=1}^N{\lambda }_i{x}_i^T{x}_i=(s^{T}s-2c^{T}s+c^{T}c)+$$

$$+(c^{T}c-2c^{T}c+\sum _{i=1}^N{\lambda }_i{x}_i^T{x}_i)=\sum _{i=1}^N{\lambda }_i(\Vert s-c{\Vert }^2+\Vert c-x_i{\Vert }^2)$$

$$\sum _{i=1}^N{\lambda }_i\Vert s-x_i{\Vert }^2=\sum _{i=1}^N{\lambda }_i(\Vert s-c{\Vert }^2+\Vert c-x_i{\Vert }^2)\le 2\cdot \sum _{i=1}^N{\lambda }_i\Vert c-x_i{\Vert }^2.$$

$$\varphi (\bar{P})={\lambda }_1\cdot \Vert c-x_{11}{\Vert }^2+{\lambda }_2\cdot \Vert c-x_{21}{\Vert }^2={\lambda }_1\cdot \Vert {\lambda }_2(x_{21}-x_{11}){\Vert }^2+{\lambda }_2\cdot \Vert {\lambda }_1(x_{11}-x_{21}){\Vert }^2=$$

$$={\lambda }_1{\lambda }_2({\lambda }_2+{\lambda }_1)\Vert x_{21}-x_{11}{\Vert }^2={\lambda }_1{\lambda }_2\Vert x_{21}-x_{11}{\Vert }^2.$$

$$\underset{z,y}{\min }$$

$$\text{s.t.}\sum _{k=1}^{|P_i|}{y}_{ijk}=$$

$$\sum _{j=1}^{|S_{\text{org}}|}{y}_{ijk}=$$

$$y_{ijk}\ge$$

$$z_j\in$$

$$\sum _{i=1}^N{\lambda }_i\Vert s-x_{i1}{\Vert }^2=\sum _{i=1}^N{\lambda }_i(\Vert c-x_{i1}{\Vert }^2+\Vert s-c{\Vert }^2)<2\cdot \sum _{i=1}^N{\lambda }_i\Vert c-x_{i1}{\Vert }^2,$$

$$\varphi (\bar{P})\le \varphi ({\bar{P}}^{\prime })\le \varphi ({\bar{P}}_{\text{org}}).$$

$$\frac{\varphi ({\bar{P}}_{\text{org}})}{\varphi (\bar{P})}=\frac{\frac{1}{2}{ϵ}^2}{\frac{3}{16}+\frac{1}{4}{ϵ}^2}\to 2.$$

$$|S_{\text{org}}|+|S_{\text{org}}|\cdot |P_{\text{sum}}|\text{ variables }\text{ and }N\cdot |S_{\text{org}}|+|P_{\text{sum}}|\text{ constraints.}$$

$$\sum _{i=1}^N{\lambda }_i\Vert s-x_{i{q}_i}^l{\Vert }^2=\sum _{i=1}^N{\lambda }_i(\Vert s-c{\Vert }^2+\Vert c-x_{i{q}_i}^l{\Vert }^2),$$

$$\sum _{i=1}^N{\lambda }_i\Vert s_j-x_{i{q}_i}^l{\Vert }^2=\sum _{i=1}^N{\lambda }_i(\Vert s_j-c{\Vert }^2+\Vert c-x_{i{q}_i}^l{\Vert }^2)=\sum _{i=1}^N{\lambda }_i\Vert s_l-x_{i{q}_i}^l{\Vert }^2.$$

$$|{\bar{P}}^{\prime }|=\sum _{l=1}^{|{\bar{P}}_{\text{org}}|}|{\bar{P}}^l|\le \sum _{l=1}^{|{\bar{P}}_{\text{org}}|}(\sum _{i=1}^N|P_i^l|-N+1)\le \sum _{l=1}^{|{\bar{P}}_{\text{org}}|}(\sum _{i=1}^N|P_i|-N+1)\le (\sum _{i=1}^N|P_i|-N+1{)}^2.$$

$$\varphi ({\bar{P}}^l)=d_l\cdot \sum _{i=1}^N{\lambda }_i\Vert c-x_{i1}^l{\Vert }^2<d_l\cdot \sum _{i=1}^N{\lambda }_i\Vert s-x_{i1}^l{\Vert }^2.$$

Full text

It is open whether a discrete barycenter can be computed in polynomial time. It is possible to find an exact barycenter through linear programming, but these programs may scale exponentially. In this paper, we prove that there is a strongly-polynomial $2$-approximation algorithm based on linear programming. First, we show that an exact computation over the union of supports of the input measures gives a tight $2$-approximation. This computation can be done through a linear program with setup and solution in strongly-polynomial time. The resulting measure is sparse, but an optimal transport may split mass. We then devise a second, strongly-polynomial algorithm to improve this measure to one with a non-mass splitting transport of lower cost. The key step is an update of the possible support set to resolve mass split.

Finally, we devise an iterative scheme that alternates between these two algorithms. The algorithm terminates with a $2$-approximation that has both a sparse support and an associated non-mass splitting optimal transport. We conclude with some sample computations and an analysis of the scaling of our algorithms, exhibiting vast improvements in running time over exact LP-based computations and low practical errors.

Keywords: discrete barycenter, optimal transport, $2$-approximation, linear programming

MSC: 90B80, 90C05, 90C46, 90C90

1 Introduction

Transportation problems for several marginals arise in applications ranging from finance and economics (Beiglböck et al., 2013; Galichon et al., 2014; Rabin et al., 2012; Pass, 2014) to physics (Buttazzo et al., 2012; Cotar et al., 2013), economics (Carlier and Ekeland, 2010; Chiaporri et al., 2010), statistics (Bigot and Klein, 2017; Boissard et al., 2015; Mileyko et al., 2011), and data analytics (del Barrio et al., 2019; Gadat et al., 2018). The so-called Wasserstein barycenters correspond to optimal solutions to these problems, and have seen much recent attention. Barycenters are intimately connected to Fréchet means in Euclidean space (Munch et al., 2015; Trouvé and Younes, 2005; Turner et al., 2014; Zemel and Panaretos, 2019), which is one of the origins of this field of research and the reason why statistical and probability notation is commonly used.

Given probability measures $P_{1},\ldots,P_{N}$ on $\mathbb{R}^{d}$ and a weight vector $\lambda=(\lambda_{1},\dots,\lambda_{N})\in\mathbb{R}^{N}_{>0}$ with $\sum_{i=1}^{N}\lambda_{i}=1$, a ($\lambda$-weighted) Wasserstein barycenter is a probability measure $\bar{P}$ on $\mathbb{R}^{d}$ which satisfies

[TABLE]

where $W_{2}$ is the quadratic Wasserstein distance and $\mathcal{P}^{2}(\mathbb{R}^{d})$ is the set of all probability measures on $\mathbb{R}^{d}$ with finite second moments. We recommend the monographs Villani (2003, 2009), and the more recent Panaretos and Zemel (2019); Peyré and Cuturi (2019), for a review of the Wasserstein distance and an overview of the literature on optimal transport problems.

1.1 Exact, approximate and heuristic algorithms

Exact barycenter computations for a set of continuous measures are intractable outside of some special cases, in particular because an evaluation of the Wasserstein distance between two given continuous measures itself already is challenging. Because of this, the literature uses several types of simplifications to facilitate practical computations. The arguably most important one is an input of discrete data:

In many applications, data is given as a set of discrete (probability) measures $P_{1},\ldots,P_{N}$ having finite support in $\mathbb{R}^{d}$. A discrete Wasserstein barycenter is a probability measure $\bar{P}$ which satisfies Eq. (1) for such measures. Discrete probability measures arise naturally in applications in operations research. The finite support often corresponds to a set of geographical locations (customers, facilities, service providers) and the measures represent data that varies over time. In computer graphics or image science, discrete probability measures are supported on a grid or arise through a discretization of the underlying space.

In Anderes et al. (2016), some theoretical results were developed for these discrete barycenters. They mirror the continuous case, established in Agueh and Carlier (2011), with a few notable exceptions. First, unlike in the continuous case, there may exist several discrete barycenters for the same set of measures. All of them have finite support and there always exists a discrete barycenter with provably sparse support. Analogously to the continuous case, a discrete barycenter always has a non-mass splitting optimal transport to each discrete marginal, i.e., each barycenter support point transports its whole mass to just a single support point for each measure. These results were proven for uniform $\lambda_{i}=\frac{1}{n}$ in Anderes et al. (2016), but are readily transferred to general $\lambda_{i}$ (Miller (2016)).

Both sparsity and non-mass split are crucial to applications. Sparsity is desirable in many applications of operations research such as facility location. Non-mass split is often imposed by physical limitations of applications, such as the design of deformable templates (Boissard et al., 2015; Jain et al., 1998; Trouvé and Younes, 2005). For example, in metal shaping, sheets of metal have to be pressed into a collection of different shapes. Each of these shapes is modeled as a measure. A ‘mean deformation’ (barycenter) is a best shape for the initial sheet of metal with respect to the energy required to mold (transport) it into all required shapes. Only because of the existence of a non-mass splitting transport, the mean deformation can indeed be transformed into each shape through only bending and stretching. See Boissard et al. (2015) for more details on such applications. We are interested in the computation of a discrete barycenter that exhibits these beneficial properties.

It is open whether the computation of an exact discrete barycenter can be done in polynomial time. It is well-known that linear programming can be used to approximate or solve optimal transport problems (Srivastava et al., 2018; Ye et al., 2017). Most importantly, exact discrete barycenters can be computed through linear programming (Anderes et al., 2016; Borgwardt and Patterson, 2020; Carlier et al., 2015). However, these programs may scale exponentially in the number of measures $N$ (see Section 2.2) and thus have not been widely considered for practical use. Much larger and faster practical computations are possible through various heuristics that fundamentally differ from LP-based approaches. Most of these algorithms are based on simplifications to the Wasserstein distance to obtain an easier objective function.

The arguably most popular tools are entropic regularization techniques, which are used to make the objective function smooth and strictly convex (Cuturi, 2013). In recent years there has been significant progress on these techniques and they led to a flurry of competitive algorithms for good approximations of barycenters in practice. For example, regularization is used in a well-behaved implementation of a gradient descent algorithm that uses information from both smoothed primal and dual optimal transport formulations (Cuturi and Doucet, 2014). Further, regularization not only greatly simplifies the underlying optimal transport problem itself, but the regularized barycenter problem then also allows the efficient computation of iterative Bregman projections (Benamou et al., 2015). Iterative Bregman methods have proven to be a competitive approach for large-scale computations (Yang et al., 2019; Ye et al., 2017); their number of variables scales roughly linearly in the number of marginals. The entropy regularized Wasserstein distance converges towards the actual Wasserstein distance in $O(\frac{1}{w})$, where $w$ is the entropic regularization factor, and the non-regularized transport cost computed with a regularized transport plan converges towards the Wasserstein distance in $O(e^{-w})$ (Benamou et al., 2015; Carlier et al., 2017; Luise et al., 2018). The factor $w$ usually is chosen empirically.

The great scalability of regularization-based methods, see for example Solomon et al. (2015), comes at the cost of a few drawbacks. First, they typically require a fixed support over which a barycenter approximation is to be computed. For grid-structured data, approaches in the literature usually just specify the underlying grid as the support set. One of the main results in this paper is that doing so, by itself, leads to an approximation error of up to $2$ – an exact optimum over the original support is only a (possibly tight) $2$-approximation. An exact barycenter for grid-structured data lies in an $N$-times finer grid. An additional challenge lies in scenarios where measures with sparse support are spread out over a large region, or where they lie in high dimension, such that it is not feasible to discretize the whole underlying space. These restrictions have led to interest in (different) approaches where the support for a barycenter approximation (Claici et al., 2018; Frogner et al., 2019; Luise et al., 2019) or a discretization of the input measures (Staib et al., 2017) is not part of the input. Second, regularization leads to fully dense solutions, which are considered an undesirable ‘blur’ in many applications. This contrasts with the search for sparse exact or approximate barycenters as done in this paper. Post-processing could be used to ‘sparsify’ a dense solution, but it is open whether this can be done in a way to retain a provable approximation guarantee or to obtain a non-mass splitting transport. And third, they sometimes exhibit poor numerical behavior as regularization decreases (Kroshnin et al., 2019); typically a fixed number of iterations is hard-coded.

There are other successful algorithms that do not rely on smoothing or regularization: for example, a non-smooth optimization algorithm based on quasi-Newton steps and the fast computation of super-gradients performs well in practice (Carlier et al., 2015). Many further examples are based on other types of simplifications of the Wasserstein distance. The so-called Radon barycenters and Sliced barycenters (Bonneel et al., 2015; Rabin et al., 2012) are restricted to special instances (Radon barycenters are only practical for data on a grid, Sliced barycenters deal with support points of uniform mass), but provide good results in low dimension. The idea is to use a Radon transform to obtain $1$-dimensional projections of the support points to lines sampled randomly, from which an expectation of the Wasserstein distance can be devised. Further, a use of the simpler $W_{1}$-distance instead of the $W_{2}$-distance leads to the so-called Beckmann problem, which allows various efficient approaches (Auricchio et al., 2019; Essid and Solomon, 2017; Solomon et al., 2014).

1.2 Contributions

In this paper, we study LP-based approaches to the discrete barycenter problem. In Section 2, we introduce some notation and recall previous related work on linear programming for the problem. In Section 3, we present and discuss our main contributions.

First, we show that an optimal measure for Eq. (1), when restricted to the union of supports of the original measures, gives a tight $2$-approximation for the barycenter problem (Theorem 3.1). This result has an immediate implication for algorithms in the literature that compute an approximate barycenter for grid-structured data: if the computation is done over the original grid itself, the algorithm does not converge to an exact barycenter, but to a best approximation of it over the grid, which can give up to a $2$-error.

Next, we exhibit that a restriction to the support of the original measures allows us to trade this small, provable approximation error for a dramatic improvement in the size of barycenter LPs: we obtain an LP-based $2$-approximation algorithm that can be set up and solved in strongly-polynomial time (Algorithm 1, Theorem 3.2), i.e., polynomial in the number of variables and constraints in the input (and their actual size does not matter). The algorithm finds a sparse approximation; to the best of the author’s knowledge, this is the first algorithm with an approximation guarantee for a sparse solution to the problem. The result shows that the barycenter problem can be efficiently approximated for any data.

The output of Algorithm 1 may not allow for a non-mass splitting transport. (Recall that the existence of such a transport is a property of all exact barycenters and is important for many applications.) Next, we present a second algorithm that improves an approximate barycenter as computed through Algorithm 1 to another measure with a non-mass splitting transport of lower cost, and prove that this computation also runs in strongly-polynomial time (Algorithm 2, Theorems 3.3 and 3.4, proofs in Appendix 0.A). This algorithm achieves the improvement by moving mass out of the union of original supports to a new, updated support set. Both Algorithms 1 and 2 work for any input and do not require the a priori specification of a fixed support set.

Finally, we use the two algorithms as the building blocks of an iterative scheme alternating between them (Algorithm 3). We prove that it terminates with a $2$-approximation with both sparse support and an associated non-mass splitting optimal transport at the same time (Theorem 3.5, proof in Appendix 0.B). The theoretical running time of this third algorithm remains open at this time; in practical computations we observe a low number of iterations (often just two to four) before termination. This behavior is reminiscent of the well-known $k$-means algorithm (Lloyd, 1982; MacQueen, 1967). Further, while we exhibit an example that shows the $2$-approximation bound for Algorithm 3 is tight in theory, we have not observed more than a $20\%$ error, respectively a multiplicative $1.2$ error, in the computations in this paper.

We conclude with sample computations and an analysis of the scaling of the algorithms in Section 4. We provide comparisons to exact, LP-based computations, and observe dramatic improvements in running time. We also provide a brief comparison to one of the most popular regularization-based algorithms in the literature (Cuturi and Doucet (2014)) for grid-structured data. However, one has to be careful in such a comparison; our algorithms have some properties that come at a significant computational cost: Algorithm 1 already is an exact, sparse solution over the original support, and its result is further refined through Algorithms 2 and 3. Further, the algorithms are numerically stable, work for any data, and work without specification of a fixed support set for the solution. As expected, simplifications to the objective function (like entropy regularized input), and the search for dense approximations over the original support, leads to worse approximation errors, but much faster running times. In fact, we discuss why grid data is an especially poor setting for our methods in view of computational speed: the size of the original support scales quadratically with the density of the grid.

We close the discussion with a ‘best-case’ example for computational speed: a large number of measures of small, overlapping support (without grid-structure). Due to a linear scaling of the LPs in the number of measures, we are then able to find solutions for thousands of measures. Such data is common in operations research applications that involve geographical locations, but algorithms in the literature are not designed to work (well) in this setting.

2 Preliminaries

We begin by recalling some background on LP-based approaches to the discrete barycenter problem, following Anderes et al. (2016) and Miller (2016). We are given a set of discrete probability measures $P_{1},\dots,P_{N}$, i.e., they have finite support in $\mathbb{R}^{d}$ and their total mass sums up to $1$. For a simple wording, we call them measures in this paper, or discrete measures to stress that they have a finite support. A set of support points with associated total mass less than $1$ will be called a partial measure. We denote the support of $P_{i}$ as $\text{supp}(P_{i})$ and the corresponding number of support points as $|P_{i}|=|\text{supp}(P_{i})|$. $|P_{i}|$ is called the size of $P_{i}$. Further, we are given a weight vector $\lambda=(\lambda_{1},\dots,\lambda_{N})\in\mathbb{R}^{N}_{>0}$ with $\sum_{i=1}^{N}\lambda_{i}=1$.

The general definition of a Wasserstein barycenter refers to a measure $\bar{P}$ on $\mathbb{R}^{d}$ which satisfies Eq. (1), i.e.,

[TABLE]

For discrete measures $P_{1},\dots,P_{N}$, one can show (Anderes et al., 2016) that all optimizers of Eq. (1) must be supported in the finite set $S\subset\mathbb{R}^{d}$ defined as

[TABLE]

$S$ is the set of weighted centroids for all possible combinations of support points, one from each measure $P_{i}$. Note that $S$ does not have to overlap with the support sets $\text{supp}(P_{i})$.

2.1 Linear programs for discrete barycenters

Setting $\mathcal{P}_{S}^{2}(\mathbb{R}^{d}):=\{P\in\mathcal{P}^{2}(\mathbb{R}^{d})|\text{ }\text{supp}(P)\subseteq S\}$, the infinite-dimensional problem in Eq. (1) can be solved by replacing the requirement $P\in\mathcal{P}^{2}(\mathbb{R}^{d})$ with $P\in\mathcal{P}_{S}^{2}(\mathbb{R}^{d})$ to obtain

[TABLE]

This yields a finite-dimensional minimization problem, which can be solved through linear programming (Anderes et al., 2016; Borgwardt and Patterson, 2020; Carlier et al., 2015). We recall this construction in two steps: we begin with the computation of the value $\phi(P_{0})=\sum_{i=1}^{N}\lambda_{i}W_{2}(P_{0},P_{i})^{2}$, i.e., the cost of an optimal transport from $P_{0}$ to all the $P_{i}$ for a given $P_{0}$. Then we make $P_{0}$ part of the optimization, too.

Let $P_{1},\dots,P_{N}$ be a set of discrete measures and let \text{supp}(P_{i})=\{x_{ik}\big{|}k=1,...,|P_{i}|\}. Further, let $P_{0}$ be another (fixed) discrete measure and let \text{supp}(P_{0})=\{x_{j}\big{|}j=1,...,|P_{0}|\}. Finally, let $d_{ik}$ be the mass of support point $x_{ik}$ in $P_{i}$ and $d_{j}$ be the mass of support point $x_{j}$ in $P_{0}$. Then we can find the value of $\phi(P_{0})$ by solving the following LP:

[TABLE]

Note that we not only find the optimal objective function value $\phi(P_{0})$, but also a corresponding (optimal) transport $y=(y_{ijk})_{i\leq N,j\leq|P_{0}|,k\leq|P_{i}|}$ between $P_{0}$ and the $P_{1},\dots,P_{N}$.

Next, the mass becomes part of the optimization. Instead of just searching for an optimal transport from a fixed measure $P_{0}$, we use a set $S_{0}$ of possible support points with associated variables that represent mass on them. By introducing variables $z=(z_{j})_{j\leq|S_{0}|}$ for the points in a given set S_{0}=\{x_{j}\big{|}j=1,...,|S_{0}|\} to denote the possible mass at $x_{j}\in S_{0}$, we obtain an LP that both finds an optimal measure $P_{0}$ supported on $S_{0}$, as well as a corresponding optimal transport:

[TABLE]

Note that the variables $z_{j}$ satisfy $z_{j}\geq 0$ and $\sum_{j=1}^{|S_{0}|}z_{j}=1$ because of satisfaction of the other constraints and because $\sum_{i=1}^{|P_{i}|}d_{ik}=1$ for all $i\leq N$. Thus, it suffices to specify $z_{j}\in\mathbb{R}$.

The above LP computes a measure represented by $z$ and a corresponding optimal transport $y$. For $S_{0}=S$, the returned $(z,y)$ represents a discrete barycenter by $z$ and a corresponding optimal transport by $y$. For $S_{0}\neq S$, we call the measure represented by $z$ an $S_{0}$-barycenter, an approximation of the barycenter in $S_{0}$, or simply an approximate barycenter when the context is clear.

2.2 Scaling of the LPs

Let us consider the size of LP (2.1). It consists of $|S_{0}|+|S_{0}|\cdot\sum_{i=1}^{N}|P_{i}|$ variables and $N\cdot|S_{0}|+\sum_{i=1}^{N}|P_{i}|$ equality constraints. For the computation of an exact barycenter, we set $S_{0}=S$. In this case, we get a worst-case bound of $|S_{0}|=\prod_{i=1}^{N}|P_{i}|$. Let now $|P_{\max}|=\max_{i=1,\dots,N}|P_{i}|$. If all measures have the same number of support points, we get $\sum_{i=1}^{N}|P_{i}|=N\cdot|P_{\max}|$ and $\prod_{i=1}^{N}|P_{i}|=|P_{\max}|^{N}$. So we have an LP of up to $|P_{\max}|^{N}+|P_{\max}|^{N}\cdot N\cdot|P_{\max}|$ variables and $N\cdot|P_{\max}|^{N}+N\cdot|P_{\max}|$ equality constraints.

A refined analysis reveals that some of the variables and constraints can be redundant. For example, if the measures overlap in some of their support points, then $|S_{0}|$ and consequently the size of the LP becomes smaller. In fact, LP (2.1) is always of polynomial size for data on a grid (Borgwardt and Patterson, 2020). However, in general one cannot rule out a scaling of the size of the LP for $S_{0}=S$ that is exponential in $N$ even if $|P_{\max}|$ is fixed, and a polynomial scaling in $|P_{\max}|$ even if $N$ is fixed. The main reason why it was possible to compute an exact barycenter for the example in Anderes et al. (2016) with only $8$ measures of $9$ support points was the fact that all measures had the same small support, which had a dramatic effect in reducing $|S_{0}|$. This highlights the potential benefit from performing an approximate computation where one replaces $S$ by a smaller set $S_{0}$.

2.3 Sparsity and non-mass split

The feasible regions of LPs (2.1) and (2.1) are bounded, and thus standard arguments of linear programming show that there always exists an optimal vertex. In a vertex, an inclusion-maximal set of variables is set to [math]. By a careful analysis of which of the variables $z_{j},y_{ijk}$ are equal to [math], it is possible to show a first favorable property: in contrast to the large number $|S|$ of possible support points, which can be up to $\prod_{i=1}^{N}|P_{i}|$, there always exists a barycenter that assigns nonzero mass to at most $\sum_{i=1}^{N}|P_{i}|-N+1$ points (Anderes et al., 2016).

Proposition 1

Let $P_{1},\ldots,P_{N}$ be discrete measures. Then for any weights $\lambda\in\mathbb{R}^{n}_{>0}$, there exists a barycenter $\bar{P}$ of these measures such that the size $|\bar{P}|$ satisfies

[TABLE]

We call a measure $\bar{P}$ that satisfies $|\bar{P}|\leq\sum_{i=1}^{N}|P_{i}|-N+1$ sparse. Proposition 1 states that there always exists a sparse barycenter. A proof is based on the existence of an optimal vertex of the polyhedron for LP (2.1) (Anderes et al., 2016; Miller, 2016). The argument also works if a support set $S_{0}\neq S$ is used. LP (2.1) then optimizes the objective function in Eq. (1) over the set $\mathcal{P}_{S_{0}}^{2}(\mathbb{R}^{d})$ of all measures $P$ with support in $S_{0}$. For these different support sets, we have the following generalization of Proposition 1.

Corollary 1

Let $P_{1},\ldots,P_{N}$ be discrete measures in $\mathbb{R}^{d}$, let $S_{0}=\{x_{j}:j=1,\dots,|S_{0}|\}\subset\mathbb{R}^{d}$, and let $\mathcal{P}_{S_{0}}^{2}(\mathbb{R}^{d})$ be the set of all measures $P$ with support in $S_{0}$. Then for any weights $\lambda\in\mathbb{R}^{n}_{>0}$, there exists an approximate barycenter $\bar{P}_{0}$ in $S_{0}$ such that the size $|\bar{P}_{0}|$ satisfies

[TABLE]

Further, for any exact barycenter $\bar{P}$ there exists a non-mass splitting optimal transport from $\bar{P}$ to the $P_{1},\dots,P_{N}$ (Anderes et al., 2016; Miller, 2016). This means that for all $x_{j}\in\text{supp}(\bar{P})$ with mass $d_{j}$ and for each $i$, there is exactly one $k$ with $y_{ijk}=d_{j}$, while $y_{ijk^{\prime}}=0$ for all $k^{\prime}\neq k$. Each support point of a barycenter only transports mass to exactly one support point in each measure. In this case, we say that a support point does not split mass or that a support point is non-mass splitting.

In fact, *any * optimal transport from a discrete barycenter $\bar{P}$ to the corresponding set of measures is non-mass splitting. While this has not been stated explicitly in (Anderes et al., 2016), it is not hard to prove: recall that the (weighted) centroid $c$ of a set of points $x_{1},\dots,x_{n}$ is the unique minimizer of a functional that measures the (weighted) summed-up squared Euclidean distances of a single point to all points in the set. This can be seen through a simple transformation

[TABLE]

which is minimal for $s^{T}s=0$, so $s=0$. If there was a barycenter support point splitting mass, it could be split into two (or more) centroids of support points in the measures of the same total mass, and the cost of transport would be strictly lower. We formally state this observation.

Proposition 2

Let $P_{1},\ldots,P_{N}$ be discrete measures, and let $\bar{P}$ be a barycenter for these measures. Then any optimal transport from $\bar{P}$ to $P_{1},\ldots,P_{N}$ is non-mass splitting.

3 Main Results

In this paper, we study approximations of the discrete barycenter problem where the set $S$, required to find an exact barycenter, is replaced by a much smaller set $S_{0}$. This is motivated by the unfavorable scaling of LP (2.1) with respect to $|S|$, respectively $|S_{0}|$; see Section 2.2.

3.1 A strongly-polynomial $2$-approximation

Recall that the set of possible support points of a discrete barycenter is

[TABLE]

which may consist of up to $\prod\limits_{i=1}^{N}|P_{i}|$ points. This is a much larger number than the size of the union of supports of the measures

[TABLE]

which satisfies $|S_{\text{org}}|\leq\sum\limits_{i=1}^{N}|P_{i}|$ with equality if and only if the supports are disjoint.

Note that the maximal size of $S_{\text{org}}$ only barely exceeds the bound given in Proposition 1. First, we show that the approximation error from searching for an approximate barycenter in $S_{\text{org}}$, i.e., setting $S_{0}=S_{\text{org}}$ in LP (2.1), can be bounded by a factor of two. This bound is tight.

Theorem 3.1

*Let $\bar{P}$ be a barycenter and let $\bar{P}_{\text{org}}$ be an approximate barycenter in $S_{\text{org}}$. Then *

[TABLE]

and this bound can become tight, i.e., there is a set of measures $P_{1},\dots,P_{N}$ and a set of weights $\lambda_{1},\dots,\lambda_{N}$ for which $\phi(\bar{P}_{\text{org}})=2\cdot\phi(\bar{P})$.

Proof

We denote the mass of a support point $c$ of a barycenter $\bar{P}$ by $d_{c}$. By Proposition 2, there is an optimal transport such that $c$ transports its mass to exactly one support point $x_{i}$ in each $P_{i}$ for all $i\leq N$. Due to optimality of $\bar{P}$, $c$ is the weighted centroid $c=\sum_{i=1}^{N}\lambda_{i}x_{i}$ of these points. Recall the discussion after Corollary 1.

Each support point $c$ contributes $d_{c}\cdot\sum_{i=1}^{N}\lambda_{i}\|c-x_{i}\|^{2}$ to the corresponding value $\phi(\bar{P})$. Let $s\in S_{\text{org}}=\bigcup_{i=1}^{N}\text{supp}(P_{i})$ be such that $\|s-c\|^{2}$ is minimal and note that

[TABLE]

for any $s$. By choice of $s$ and the fact that $x_{i}\in\text{supp}(P_{i})$, we know $\|s-c\|^{2}\leq\|c-x_{i}\|^{2}$ for all $i\leq N$, so we obtain

[TABLE]

Thus the transport from $s$, instead of from $c$ itself, introduces an approximation error of $2$, i.e., each such $s$ contributes at most $2\cdot d_{c}\sum_{i=1}^{N}\lambda_{i}\|c-x_{i}\|^{2}$ to the value $\phi(\bar{P}_{\text{org}})$. As this holds for all weighted centroids $c\in\text{supp}(\bar{P})$ and corresponding closest $s\in S_{\text{org}}$, this shows the existence of a measure $\bar{P}_{\text{org}}\in\mathcal{P}_{\text{org}}^{2}(\mathbb{R}^{d})$ with approximation error $2$ with respect to $\phi$.

It remains to prove that the bound can be tight. We do so through a simple example. Let $P_{1},P_{2}$ be two measures with a single support point $x_{11}\in\text{supp}(P_{1})$, $x_{21}\in\text{supp}(P_{2})$, each of mass $1$. Then $\bar{P}$ consists of the single support point $c=\lambda_{1}x_{11}+\lambda_{2}x_{21}$ of mass $1$ and thus

[TABLE]

In contrast, the restriction of an approximate barycenter $\bar{P}_{\text{org}}$ to possible support $S_{\text{org}}=\{x_{11},x_{21}\}$ would give $\phi(\bar{P}_{\text{org}})=\min\{\lambda_{1},\lambda_{2}\}\cdot\|x_{21}-x_{11}\|^{2}$. Note $\lambda_{1}\cdot\lambda_{2}\geq\frac{1}{2}\min\{\lambda_{1},\lambda_{2}\}$, with equality if and only if $\lambda_{1}=\lambda_{2}=\frac{1}{2}$. In this case, $\phi(\bar{P}_{\text{org}})=2\cdot\phi(\bar{P})$. $\square$

The difference between the support for an exact barycenter and for an approximation in $S_{\text{org}}$ is highlighted in Figure 1: the first two rows show four handwritten digits scanned into a $16\times 16$ grid. (See LeCun et al. (1998) for some information on this data set.) These are the measures $P_{1},\dots,P_{4}$. The varying shades of grey indicate different masses at the support points of the grid (the darker, the larger the mass). The masses for each measure add up to $1$. The bottom row depicts an exact barycenter and a $2$-approximation in the original $16\times 16$ grid (for all $\lambda_{i}=\frac{1}{4}$). The support grid for the exact barycenter is four times finer, a $(4\cdot 16-3)\times(4\cdot 16-3)=61\times 61$ grid.

We formally denote the choice of $S_{\text{org}}$ in LP (2.1), as performed for Theorem 3.1, as Algorithm 1. Note that the algorithm is stated to compute an optimal vertex of the feasible region. (For a convenient wording, we will say that Algorithm 1 is used with a given different support $S_{0}$ as input when we require an optimal vertex, and not just any optimal solution, of LP (2.1).) The search for a vertex guarantees that the sparsity condition stated in Corollary 1 is satisfied, so the returned measure is not only an approximate barycenter in $S_{\text{org}}$, but also sparse. However, it is possible that any corresponding optimal transport splits mass (which would not happen for an exact barycenter). Here is an example of this type.

Example 1

Consider the two measures $P_{1},P_{2}$ depicted at the top of Figure 2. The different radii of the filled circles represent the different masses on the support points. Let $\lambda_{1}=\lambda_{2}=\frac{1}{2}$

Measure $\bar{P}_{\text{org}}\in\mathcal{P}_{\text{org}}^{2}(\mathbb{R}^{d})$ (second row, left) is an optimal barycenter approximation in $S_{\text{org}}$. It consists of only two support points, while $P_{1}$ and $P_{2}$ have three support points. Thus, there exists a support point of $\bar{P}_{\text{org}}$ that splits mass in any transport, including the unique optimal one (second row, right): the top support point of $\bar{P}_{\text{org}}$ transports its mass $\frac{1}{2}$ in two parts $\frac{1}{4}$ to two support points of $P_{1}$; the same happens in the bottom part with respect to $P_{2}$. Such a split of mass does not happen for an exact barycenter (third row). $\square$

We would like to note that the $2$-bound in Theorem 3.1 can only be tight in very special cases. Let $s\in S_{\text{org}}$ be such that $\|s-c\|^{2}$ is minimal for a given weighted centroid $c\notin S_{\text{org}}$ transporting to $x_{i1},\dots,x_{N1}$ with $x_{i1}\in P_{i}$ . Then the approximation error $2$ is not tight if $\|c-x_{i1}\|^{2}\neq\|c-x_{j1}\|^{2}$ for any $i\neq j$. This holds because then

[TABLE]

as there has to be an $i\leq N$ with $\|c-x_{i1}\|^{2}>\|s-c\|^{2}$. Further, it is easy to give examples where even an exact barycenter is actually contained in $S_{\text{org}}$. For example, take arbitrary measures $P_{1},\dots,P_{N}$, compute their barycenter $\bar{P}$, and consider a new collection $P_{1},\dots,P_{N},P_{N+1}$ with $P_{N+1}=\bar{P}$. Then one has $\phi(\bar{P}_{\text{org}})=\phi(\bar{P})$.

Next, we prove that Algorithm 1 runs in strongly-polynomial time. Recall that LPs are generally solvable in weakly-polynomial time, i.e., the number of arithmetic operations is polynomial in the length of a bit representation of the input. This means polynomiality in the number of variables and constraints, as well as in the logarithm of absolute values of numbers in the input. In contrast, a strongly-polynomial running time restricts polynomiality to only the number of variables and constraints.

Theorem 3.2

For all rational input, a $2$-approximate barycenter can be computed in strongly-polynomial time.

Proof

A proof of strong polynomiality for Algorithm 1 is based on exhibiting that LP (2.1) is of strongly-polynomial size, and that its parameters can be computed in strongly-polynomial time. General LPs are known to be solvable in weakly-polynomial time. However, it suffices to restrict the dependency of the running time only to the parameters that appear in the constraint matrix; the numbers in the objective function or the right-hand side of the constraints do not matter (Tardos, 1986).

First, note that the constraint matrix of LP (2.1) for $S_{\text{org}}$ only consists of entries in $\{-1,0,1\}$. For the claim of strongly-polynomial solvability, it only remains to prove that the number of variables and constraints of the LP is strongly-polynomial in the size of the input, and that the parameters that appear in the objective function and right-hand sides can be computed from the original input in strongly-polynomial time.

Let $\mathcal{I}$ be an instance of the problem and let $|\mathcal{I}|$ be the number of bits to represent the input. Any representation of the input $\mathcal{I}$ has to satisfy $|\mathcal{I}|\geq\sum_{i=1}^{N}|P_{i}|$. As $|S_{0}|=|S_{\text{org}}|\leq\sum_{i=1}^{N}|P_{i}|\leq|\mathcal{I}|$, LP (2.1) indeed has a strongly-polynomial number of constraints and variables.

The actual numbers that appear in the LP are of types $\lambda_{i}$, $d_{ik}$, or $\|x_{j}-x_{ik}\|^{2}$. The $\lambda_{i}$ and $d_{ik}$ appear directly in the input, and so do the vectors $x_{j}$ and $x_{ik}$. As we use rational input, $\|x_{j}-x_{ik}\|^{2}=(x_{j}-x_{ik})^{T}(x_{j}-x_{ik})$ is a rational number derived by the sum over products of pairs of coefficients in $x_{j}$ and $x_{ik}$. This implies that $\|x_{j}-x_{ik}\|^{2}$ can be computed in strongly-polynomial time (polynomial in $\text{log }x_{j}+\text{log }x_{ik}$) and represented in a number of bits that is strongly-polynomial in the number of bits of the original representation of $x_{j},x_{ik}$. This proves the claim. $\square$

3.2 Recovery of Non-Mass Split

Next, we design an algorithm that begins with a (sparse) $2$-approximate barycenter computed by Algorithm 1. The algorithm improves it to another measure supported on a subset of $S$ (instead of $S_{\text{org}}$), for which there exists a non-mass splitting transport of lower cost, i.e., the approximation error can only become better. Algorithm 2 sums up the approach in pseudocode. We here describe the algorithm in some detail; additional technical details are given in the proof of Theorem 3.3 in Appendix 0.A.

The algorithm greedily breaks up each support point (that splits mass) of the approximate barycenter into several non-mass splitting support points (Steps $1-3$). In the end, all of the non-mass splitting support points are combined to a new measure (Step $4$). The preprocessing performed in Step $2$ guarantees that the non-mass split property for each support point in Step $3$ transfers to a non-mass splitting transport for the new measure constructed in Step $4$. Figure 3 shows a run of the algorithm, which is discussed in more detail as Example 2 at the end of the section.

Step 1. First, the approximate barycenter $\bar{P}_{\text{org}}$ and measures $P_{1},\dots,P_{N}$ are broken up into disjoint parts; each part corresponds to a support point $s_{l}=x_{t_{l}}$ in the approximate barycenter. By construction, each $P_{i}^{l}$ consists of those support points in $P_{i}$ to which $s_{l}$ transports mass. The mass of a support point in $P_{i}^{l}$ equals the mass it receives as transport from $s_{l}$. Then we assign new indices to the support points in $P_{i}^{l}$ and their masses for a simpler notation, so we do not have to refer to the original $z$ or $y$ in the other steps.

Lexicographic Ordering. Step $2$ and Step $3$ are based on the construction of so-called lexicographically maximal vectors. A vector $a=(a_{1},\dots,a_{n})$ is lexicographically larger than a vector $b=(b_{1},\dots,b_{n})$ if there is an index $j\leq n$ such that $a_{j}>b_{j}$, and $a_{i}\geq b_{i}$ for all $i<j$. For example, the vector $a=(2,2,0,1)$ is lexicographically larger than $b=(2,1,5,10)$. Lexicographic maximality with respect to a set states that there is no lexicographically larger vector in the set. Note that the term gives rise to a total ordering.

The intuition for the construction of lexicographically maximal vectors is to resolve ties. This is necessary in two different settings: in Step 2, as much mass as possible is greedily shifted to support points of lower indices; in Step 3, a lexicographically decreasing sequence of weighted centroids is created from each support point. Together, these two steps make sure that all the weighted centroids that are merged to form $\bar{P}^{\prime}$ in Step 4 are distinct and only transport to a single support point in each measure, implying the existence of a non-mass splitting transport.

Step 2. Step $2$ iteratively transforms $(d_{1},\dots,d_{r})$ to be lexicographically larger and larger while retaining an approximate barycenter supported in $\text{supp}(\bar{P}_{\text{org}})$ (that is, the cost of an optimal transport does not increase). It does so via a greedy scheme, where mass is moved to support points in $\text{supp}(\bar{P}_{\text{org}})$ with the lowest indices, until this is not possible anymore. We call a $(d_{1},\dots,d_{r})$ that is not altered by Step $2$ (anymore) greedily lexicographically maximal. Note that such a vector need not be lexicographically maximal among all approximate barycenters with the same support, but this is enough for our purposes.

The two loops for $l$ and $j$ establish an order for checking whether mass can be moved from $s_{l}$ to $s_{j}$ while keeping optimality over $\text{supp}(\bar{P}_{\text{org}})$. The indices $q_{i}=\text{arg}\max_{q\leq|P_{i}^{l}|}(s_{j}-s_{l})^{T}x^{l}_{iq}$ selected in $2a)$ identify support points in the $P_{i}^{l}$ that lie the furthest in direction of $s_{j}-s_{l}$.

Note $c$ satisfies $\|c-s_{l}\|^{2}\leq\|c-s_{j}\|^{2}$ (because of optimality of $\bar{P}_{\text{org}}$) and it is a maximizer of $\|c-s_{l}\|^{2}-\|c-s_{j}\|^{2}\leq 0$. If $\|c-s_{l}\|^{2}=\|c-s_{j}\|^{2}$, which is checked in $2b)$, then mass can be shifted from $s_{l}$ to $s_{j}$ to make $(d_{1},\dots,d_{r})$ lexicographically larger, while keeping optimality. The remainder of $2b)$ is a technical description of this shift of mass.

Step 3. Next, we perform a (greedy) routine to spread out the mass of each $s_{l}$ to several support points. We do so by picking a set of lexicographically maximal support points $x^{l}_{iq_{i}}$ in each $P_{i}^{l}$ (i.e., we pick an $x^{l}_{iq_{i}}$ with a largest first coordinate, and among those one with a largest second coordinate, and so on). Then we move mass $d_{\text{min}}$ to the weighted centroid $c=\sum_{i=1}^{N}\lambda_{i}x^{l}_{iq_{i}}$, where $d_{\text{min}}$ is the minimal mass among the $d^{l}_{q_{i}}$. We repeat this scheme until the whole mass of $s_{l}$ has been spread out. The result is a partial measure $\bar{P}^{l}$ that has a non-mass splitting transport by construction. Then we continue with the next support point.

Step 4. Finally, we combine the partial measures $\bar{P}^{l}$ from Step $3$ to a new measure. It is at least as good an approximation of an exact barycenter as $\bar{P}_{\text{org}}$: Step $1$ and $2$ do not change the cost of transport. In Step $3$, for any chosen set of support points $x^{l}_{iq_{i}}$ we put the corresponding mass on their weighted centroid, which is best-possible (at least as good as transport from $s_{l}$).

We sum up the favorable properties of the algorithm in Theorem 3.3. In addition to the existence of a non-mass splitting transport, and keeping a $2$-approximation error, we are able to bound the size of the support by the square of the bound in Proposition 1. We do not prove that the returned measure is an approximate barycenter (which implies optimality over the given support by definition). The associated non-mass splitting transport is trivial to construct, but we do not prove that this transport is optimal. Due to this, we have to be careful in the wording of the following statements (Theorems 3.3 and 3.4). A detailed proof is given in Appendix 0.A.

Theorem 3.3

Algorithm 2 returns a measure $\bar{P}^{\prime}$ supported on a subset of $S$ with $\phi(\bar{P}^{\prime})\leq 2\cdot\phi(\bar{P})$ and there is a non-mass splitting transport realizing this bound. Further $|\bar{P}^{\prime}|\leq(\sum_{i=1}^{N}|P_{i}|-N+1)^{2}$.

Let us discuss a small example for Steps $2-4$ of the algorithm.

Example 2

We revisit the measures $P_{1}$ and $P_{2}$ used for Example 1 and again let $\lambda_{1}=\lambda_{2}=\frac{1}{2}$. They receive their mass transported from two fixed support points $s_{1},s_{2}$ of mass $d_{1}=\frac{1}{4}$, $d_{2}=\frac{3}{4}$ (second row, left). Note that $s_{1},s_{2}\notin S_{\text{org}}$, which may happen in later iterations of Algorithm 3 (Section 3.3), where Algorithm 2 is used as a subroutine. (For this example, this does not matter.)

The two central points, which receive their mass from $s_{2}$, have a centroid $c$ that is equally far from $s_{1}$ and $s_{2}$ (second row, right). These two points would be selected in Step $2a)$ of Algorithm 2 and their mass shifted from $s_{2}$ to $s_{1}$ in Step $2b)$. Then $d_{1}=\frac{3}{4}$, $d_{2}=\frac{1}{4}$ (third row, left).

In Step $3$, the mass of $s_{1}$ and $s_{2}$ is spread out to a set of centroids that transport to just a single support point in each measure. The result for $s_{1}$ is depicted in the third row (right). By lexicographically maximal choice of the points in the measures, the central support point of mass $\frac{1}{2}$ is constructed first, followed by the left one of mass $\frac{1}{4}$. $s_{2}$ is not changed, because it already is the centroid of a set of single support points in each measure (fourth row, left).

These partial measures are combined to form measure $\bar{P}^{\prime}$ in Step $4$ (fourth row, right) and the algorithm stops. We actually found an exact barycenter, which is not the case in general. $\square$

We close our discussion of Algorithm 2 by showing that it runs in strongly-polynomial time. The quite technical proof is given in Appendix 0.A.

Theorem 3.4

For all rational input, a measure can be computed in strongly-polynomial time that is a $2$-approximation of a barycenter and for which there is a non-mass splitting transport realizing this bound.

3.3 An Iterative Improvement

Finally, we combine Algorithms 1 and 2 to an iterative scheme, denoted as Algorithm 3. The algorithm begins by computing an approximate barycenter in $S_{\text{org}}$ using Algorithm 1. Then Algorithm 2 is used to spread out its support points to find a new measure $\bar{P}^{\prime}$ of better approximation error and that allows for a non-mass splitting transport. We set $S_{\text{org}}=\text{supp}(\bar{P}^{\prime})$ and repeat Algorithm 1 to find an optimal approximate barycenter over this new support (in other words, an optimal vertex of LP (2.1) over the new support is found). Then its support points are spread out again using Algorithm 2. This scheme is repeated until there is no improvement anymore.

After a finite number of iterations, the algorithm terminates with a sparse $2$-approximate barycenter supported on a subset of $S$, and with a non-mass splitting optimal transport. This is a provable approximation that possesses both favorable properties of an exact barycenter, sparsity and non-mass split, at the same time.

Theorem 3.5

Algorithm 3 returns an approximate barycenter $\bar{P}^{\prime}$ supported on a subset of $S$ for which $\phi(\bar{P}^{\prime})\leq 2\cdot\phi(\bar{P})$, where $\bar{P}$ is a barycenter, and there is a non-mass splitting optimal transport realizing this bound. Further $|\bar{P}^{\prime}|\leq\sum_{i=1}^{N}|P_{i}|-N+1$.

We prove Theorem 3.5 in Appendix 0.B. Let us take a closer look at the approximation error of Algorithm 3. We distinguish three different measures: $\bar{P}$ is an exact barycenter, $\bar{P}_{\text{org}}$ is an approximate barycenter in $\mathcal{P}_{\text{org}}^{2}(\mathbb{R}^{d})$ , and $\bar{P}^{\prime}$ is the solution of Algorithm 3. By optimality of $\bar{P}$ and $\bar{P}_{\text{org}}$ with respect to $\phi$ in their respective support, we have

[TABLE]

We are particularly interested in the gap between $\phi(\bar{P})$ and $\phi(\bar{P}^{\prime})$. Theorem 3.1 states $\phi(\bar{P}_{\text{org}})\leq 2\cdot\phi(\bar{P})$. Thus the whole sequence of inequalities is bounded by a total approximation factor of $2$. This implies that if $\alpha\phi(\bar{P}^{\prime})=\phi(\bar{P}_{\text{org}})$ for some $\alpha\geq 1$, then $\phi(\bar{P}^{\prime})\leq\frac{2}{\alpha}\phi(\bar{P})$. Informally, Algorithm 3 already begins with a provable $2$-approximation and any improvement towards $\bar{P}^{\prime}$ allows for the statement of a better approximation guarantee.

In practice, one obtains a strictly better approximation factor than $2$ for essentially all real-world problems using Algorithm 3. But there exist worst-case examples, such as the following, that show the bound is tight.

Example 3

Consider the example depicted in Figure 4. Four measures $P_{1},...,P_{4}$ are shown in the top row, $P_{2}$ and $P_{3}$ are depicted in the center. Note $P_{2}=P_{3}$. Each of the measures consists of two support points of mass $\frac{1}{2}$. Let $\epsilon>1$. $P_{1}$ is supported on coordinates $(-\epsilon,0)$ and $(\epsilon,1)$, $P_{2}$ and $P_{3}$ are supported on $(0,0)$ and $(0,1)$, and $P_{4}$ is supported on $(-\epsilon,1)$ and $(\epsilon,0)$, where $\epsilon>1$. For increasing $\epsilon$, the horizontal distance of the support points of $P_{1}$ and $P_{4}$ to those of $P_{2},P_{3}$ increases proportionally (second row).

Let $\lambda_{i}=\frac{1}{4}$ for $i=1,\dots,4$. Independently of $\epsilon$, an approximate barycenter $\bar{P}_{\text{org}}$ in $S_{\text{org}}$ is identical to $P_{2}=P_{3}$, (third row, left). A corresponding optimal transport sends the mass to the support points in the same ‘layer’ (third row, middle). Note that the support points are already the (weighted) centroids of the points they transport to, and that the transport is non-mass splitting. Because of this, Algorithm 3 stops without any change to $\bar{P}_{\text{org}}$ at the end of the first iteration.

The cost of transport for $\bar{P}_{\text{org}}$ is $\phi(\bar{P}_{\text{org}})=\frac{1}{4}\cdot 2\epsilon^{2}=\frac{1}{2}\epsilon^{2}$. (Recall $\lambda_{i}=\frac{1}{4}$ for all $i$.) An exact barycenter $\bar{P}$ (third row, right) and a corresponding optimal transport (fourth row) are strictly better. The coordinates for the two support points are $(-\frac{1}{2}\epsilon,\frac{3}{4})$ and $(\frac{1}{2}\epsilon,\frac{1}{4})$. The cost of transport is $\phi(\bar{P})=\frac{1}{4}\cdot(\frac{3}{4}+\epsilon^{2})=\frac{3}{16}+\frac{1}{4}\epsilon^{2}$. For $\epsilon\rightarrow\infty$,

[TABLE]

Thus the error bound goes to $2$. $\square$

In Section 4 we conclude the paper with a discussion of the theoretical scaling of our algorithms and some observations on practical computations. In our implementation, we use some tweaks for a speed-up of Algorithm 3. First, we perform Step $3$ of Algorithm 2 as the exact computation of a barycenter $\bar{P}^{l}$ when the number $N_{l}$ of support points to which a given $s_{l}$ transports is low. This leads to a better approximation bound for $\bar{P}^{\prime}$ at the end of each iteration of Algorithm 3 and a lower number of iterations overall. Further, this leads to a smaller support for the LPs in the second iteration and beyond.

Recall $N_{l}$ is bounded below by $N$, as the support point transports mass to at least one support point in each measure. We tried different values for $N_{l}$ for the cutoff to an exact barycenter computation $\bar{P}^{l}$. There is a tradeoff between each run of Algorithm 2 taking longer and a reduction in the total number of iterations. We observed good results for $N_{l}$ between $2\cdot N$ and $4\cdot N$ when $N\ll|S_{\text{org}}|$ (common for grid-structured data) and $N_{l}$ between $N+\log{N}$ and $N+4\log{N}$ when $N\gg|S_{\text{org}}|$ (common for data in general position). For these values, each run of Algorithm 2, respectively Step $2$ of Algorithm 3, takes slightly longer, but the running time remains negligible compared to Step $1$ in each iteration. For the computations in Section 4, we chose cutoffs of $2\cdot N$ and $N+\log{N}$. We observed a noticable positive impact, dropping the total running time of Algorithm 3 by about $20\%$ on average.

Second, we explicitly construct the associated transport devised in Algorithm 2 and use it for a warm-start of the subsequent Step $1$ of Algorithm 3. The LPs of Step $1$ are solved through a primal simplex method. By construction, the transport from the previous iteration is not only a feasible vertex of the new primal LP, but already close to the new optimum. (It is a set of weighted centroids transporting to their respective support points, after all.) Thus, the subsequent primal LP can be warm-started and finding the exact optimum over the new support is much faster than solving the LP from scratch. This is a crucial part of the implementation, as the LPs in later iterations can have millions of variables and otherwise would be slow to solve. With this tweak, the setup of the LP in iteration $2$ remains as a bottleneck. The LPs in later iterations can be set up through an update of the previous one. Because of the warm-start, the actual solution time of the LPs is negligible in comparison to the setup time.

4 Sample Computations and Scaling

We implemented Algorithms 1, 2, and 3 in the Julia language using Clp as linear programming solver. Julia is a modern programming language for high-performance numerical computing that provides a competitive tradeoff between efficient, but cumbersome low-level languages (C, C++) and easy-to-use, but typically slow high-level languages (Python, Matlab) (Bezanson et al., 2017; Lubin and Dunning, 2015). A primal simplex method is called in Clp for the availability to warm-start iterations. The algorithms were run on a standard laptop (Win 10, 32GB memory, i7-6820).

Our sample computations are on two representative types of data: the MNIST database of handwritten digits (widely-used for benchmarking) (LeCun et al., 1998) for grid-structured data, as well as the firehouse example from Anderes et al. (2016) (and randomly generated larger instances) for data in general position, i.e., where $S$ would be exponentially-sized. Together, these two settings cover most applications in practice: grid-structured data is common in machine learning and there is a wealth of algorithms for this setting. Data in general position typically arises when geographical locations are involved, like in many applications of operations research. The algorithms in this paper work on any data, and for any choice of $\lambda$. In contrast, working with data in general position, differing masses on the support points, or a non-uniform weight vector makes most of the algorithms in the literature impractical. Our goal is to identify for which data the practical performance of our algorithms is the most favorable. As we will see, it is data without an underlying grid-structure, with a small support and a large number of measures, that forms a best-case scenario. (We also treat the less favorable, but common grid setting in detail, for the sake of completeness.) We start with some sample runs for both types of data and then turn to the theoretical and practical scaling of computations.

4.1 Sample runs of Algorithm 3 and Observations

We begin by performing a sample run of Algorithm 3 in a grid setting, using the four digits representing number six in a $16\times 16$ grid depicted in Figure 5. They have a barycenter depicted at the bottom of the figure (for all $\lambda_{i}=\frac{1}{4}$).

Figure 6 shows the stages of a run of Algorithm 3 for this input. Each row shows one of the iterations. The approximate barycenter in the original support is already a $1.142$–approximation of the exact barycenter (top left), i.e., $\phi(\bar{P}_{\text{org}})\leq 1.142\cdot\phi(\bar{P})$, which we denote as an additive $14.2\%$-error in the figure. The first split-up using Algorithm 2 (Steps $2$ to $4$) gives an improvement to a $4.3\%$-error (top right). This is further improved to a $2.0\%$-error (in Step $1$ of Iteration $2$) by computing an optimum over the support of the previous approximation (bottom left). Now the algorithm terminates, because all of the support points of this approximate barycenter are already the weighted centroids of the support points to which they transport mass, and there is no mass split. Algorithm 3 completes in about $10$ seconds on average for a set of four measures ($9.6$ seconds for the above example). In contrast, the computation of an exact barycenter takes roughly $120$ seconds.

Next, we perform a sample run for data in general position. Here we begin with the input for the computations in Anderes et al. (2016). There are $8$ measures with the same $9$ support points of varying masses. Figure 7 shows two of the measures. Circles of larger radius indicate higher mass.

Figure 8 shows the first $\bar{P}_{\text{org}}$ and the result of the first run of Algorithm 2. The approximate barycenter in the original support is a $1.102$–approximation of the exact barycenter, i.e., there is a $10.2\%$-error. This is improved to a $1.9\%$-error in the split-up using Algorithm 2. The second iteration of Algorithm 3 does not improve the solution anymore and the algorithm terminates.

This run completes in less than a second. In contrast, the computation of an exact barycenter for such a small problem size already is surprisingly hard: despite all measures having the same support, the set $S$ is of exponential size. The LP for an exact computation has $939510$ variables and $103032$ constraints and takes roughly $500$ seconds to solve (Anderes et al., 2016).

Observations. Both of these sample runs are representative in a couple of ways. The approximation error for $\bar{P}_{\text{org}}$ and the first improvement to $\bar{P}^{\prime}$ using Algorithm 2 are already much better than the guaranteed bound of $2$. In the computations in Section 4.2, we have not encountered a run with an approximation error worse than $20\%$ for the initial $\bar{P}_{\text{org}}$ or $8.7\%$ for the initial improvement to $\bar{P}^{\prime}$. The improvement between $\bar{P}_{\text{org}}$ and $\bar{P}^{\prime}$ in the first iteration is significant. However, the additional iterations of Algorithm 3 only perform minor improvements on the approximation factor. In fact, the example in Figure 6 shows one of the largest improvements after the first iteration observed in all our computations.

Only two parts contribute significantly to the total running time of Algorithm 3: the first run of Algorithm 1 and the setup of the LP for Step $1$ in the second iteration. Together, these accounted for more than $80\%$ of the total running time. Only the first LP is run in full; the later LPs can be warm-started. Further, the setup of LPs in iterations $3$ or later is a simple update from the previous iteration. The split-up of mass in Algorithm 2 is negligible in running time.

An initial run of Algorithm 1 is unavoidable in all situations, but is efficient through the use of support set $S_{\text{org}}$. The setup of the second LP can be computationally expensive, because $S_{\text{org}}$ is replaced by the larger $\text{supp}(\bar{P}^{\prime})$. However, we observed that in practice $|\bar{P}^{\prime}|$ does not only satisfy the guaranteed bound $|\bar{P}^{\prime}|\leq(\sum_{i=1}^{N}|P_{i}|-N+1)^{2}$ (Theorem 3.3), but remains close to $\sum_{i=1}^{N}|P_{i}|-N+1$, the bound for $|P_{\text{org}}|$ (Theorem 3.5).

The main benefit of a full run of Algorithm 3 is the recovery of the combination of a sparse support and a non-mass splitting transport. If these are not crucial properties for an application, we recommend performing just a single iteration of Algorithm 3, i.e., a single run of Algorithms 1 and 2, for faster computations.

4.2 Scaling of Algorithm 1 and runs for Algorithm 3

Finally, we study the scaling of Algorithm 1 and the practical running time of our algorithms. Algorithm 1 is the main pillar of the two viable approaches in this paper for practical computations: either a single run of Algorithm 1 and 2 or a full run of Algorithm 3, if the problem size allows. It is the only LP without a warm-start. Algorithm 1 is based on an LP formulation using $S_{\text{org}}$ as the set of possible support points. Before we turn to more computations, let us take a closer look at the number of variables and constraints in this LP to set up proper expectations. Using $|P_{\text{sum}}|=\sum_{i=1}^{N}|P_{i}|$, this LP has

[TABLE]

Two types of scaling are of interest: scaling the number $N$ of measures and scaling $|S_{\text{org}}|$. Note $N\cdot|S_{\text{org}}|\geq|P_{\text{sum}}|$, so the number of constraints scales linearly in $N$ and $|S_{\text{org}}|$. Further, note that the number of constraints is always lower than the number of variables ($|P_{\text{sum}}|\geq N$), and often dramatically so ($|P_{\text{sum}}|\gg N$).

The dominating factor in the number of variables is $|S_{\text{org}}|\cdot|P_{\text{sum}}|$. If the support sets of the $P_{i}$ are disjoint, then $|S_{\text{org}}|=|P_{\text{sum}}|$ and $|S_{\text{org}}|\cdot|P_{\text{sum}}|=|S_{\text{org}}|^{2}=|P_{\text{sum}}|^{2}$; the number of variables scales quadratically in $|S_{\text{org}}|$. For the sake of a simple analysis, we assume all measures have the same number of support points $|P_{\text{max}}|$. Then $|P_{\text{sum}}|=N\cdot|P_{\text{max}}|$ and the scaling of the number of variables is quadratic in $N$.

4.2.1 Scaling for grid-structured data

First, we consider scaling for (two-dimensional) grid-structured data. Let $K$ denote the number of grid points in each direction. The grid has $K^{2}$ points, which is an upper bound on $|P_{\text{max}}|$ and on $|S_{\text{org}}|$. Thus the number of variables is bounded by $K^{2}+K^{2}\cdot(N\cdot K^{2})$ and the number of constraints is bounded by $N\cdot K^{2}+N\cdot K^{2}=2(N\cdot K^{2})$. Note that the bound on the number of variables is roughly $K^{4}\cdot N$.

While the actual sizes of $|P_{\text{max}}|$ and $|S_{\text{org}}|$ are usually significantly smaller than $K^{2}$, they typically remain a linear fraction of $K^{2}$ (for MNIST digits between $\frac{1}{5}$ and $\frac{1}{3}$), and so these bounds immediately imply two types of consequences. First, the scaling with the number $N$ of measures is linear. Second, doubling the density $K$ of the underlying grids will multiply the number of variables by $16$.

Let us turn to some computations. In Table 1, we report on average errors and completion times for a large set of runs of Algorithm 3. For each of these runs, we used random samples of $16\times 16$ digits as the measures. Each row is based on the data from a total of $100$ runs, $10$ for each digit $0,1,\dots,9$. The table lists the initial error and time for the computation of an approximate barycenter $\bar{P}_{\text{org}}$ in $S_{\text{org}}$ (Algorithm 1), the error for the first $\bar{P}^{\prime}$ (Algorithm 2), the error, time, and number of iterations for a full run of Algorithm 3, and the time for an exact barycenter computation.

We have been able to run Algorithm 3 for up to $40$ measures in less than ten minutes, and Algorithm 1 for up to $100$ measures in less than twenty minutes. The big difference is the setup of the LP for the second iteration. In contrast, we have not been able to perform the computation of an exact barycenter for more than $8$ measures (within a fixed time limit of four hours), even using some refinements to an exact barycenter computation (Borgwardt and Patterson, 2020). This is the reason for the ”n/a” entries in the table, where an approximation error is not available because of the lack of an exact solution. For increasing number of measures, the difference between the running times of Algorithm 3 and an exact computation becomes dramatic, even though grid-structured data, in fact, is a scenario where the exact LP does not scale exponentially (recall an exact barycenter is contained in an $N$-times finer grid).

The first row shows numbers on random samples of four measures, as in the example depicted in Figures 5 and 6. We observed a termination of Algorithm 3 after an average of $2.2$ iterations. This low number of iterations is not surprising because of the low initial error and our practical implementation of Step $3$ of Algorithm 2; see Section 3.3. The same effects extend to larger computations, where the approximation error of the initial $\bar{P}_{\text{org}}$ is already low, most of the further improvement already happens towards the first $\bar{P}^{\prime}$, and less than $5$ iterations were necessary on average. The times include setup of the problems. We did not observe a clear pattern with respect to the errors for the first $\bar{P}_{\text{org}}$ and $\bar{P}^{\prime}$ or the final approximation, but the average number of iterations of Algorithm 3 increases slightly with the number of measures.

Computations in denser grids quickly become impossible for the algorithms in this paper, due to the quadratic scaling of the underlying LP with respect to $|S_{\text{org}}|$. Recall that doubling the density $K$ of a grid multiplies the number of variables by $16$. Figure 9 shows the results of Algorithm 1 for four measures in a $32\times 32$ grid and a $64\times 64$ grid. The computations took about $5$ minutes, respectively $92$ minutes. The $64\times 64$ example exceeds $10$ million variables and is only solvable because of the extremely low number of constraints.

For grid-structured data, there are many algorithms in the literature that are much faster options to tackle larger problem instances (Bonneel et al., 2015; Cuturi, 2013; Cuturi and Doucet, 2014; Peyré and Cuturi, 2019). In Table 2, we compare approximation errors and computational speed for the first $\bar{P}_{\text{org}}$ and $\bar{P}^{\prime}$ to the widely-used, regularization-based algorithm from Cuturi and Doucet (2014). The algorithm uses Sinkhorn distances to simplify the objective function and leads to dense approximations $\bar{Q}$ over $S_{\text{org}}$. Recall that $\bar{P}_{\text{org}}$ is already an exact, sparse solution over $S_{\text{org}}$, and is further refined to $\bar{P}^{\prime}$. Thus, the approximation error for $\bar{Q}$ always has to be worse than for $\bar{P}_{\text{org}}$, and the gap becomes signficant when compared to $\bar{P}^{\prime}$.

For regularization-based algorithms, there are several parameters that allow for a tradeoff of computational speed and quality of result. We used the recommended settings for MNIST data from Cuturi and Doucet (2014). Table 2 shows exact approximation errors for the small instances for which an exact barycenter computation was possible (up to $8$ measures), as well as the ratio $\phi(\bar{Q})/\phi(\bar{P})$ between cost of transport for $\bar{Q}$ and $\bar{P}$ for all instances. Our computations revealed a ratio of about $1.18$, respectively a gap of $18\%$, throughout. In view of computational speed, however, the regularization-based algorithm scales dramatically better with the number of measures than our approach, as expected. We observe a near-linear increase of running time, which sharply contrasts with a linear scaling of the size of the LP for Algorithm 1.

Of course, one has to be careful in this comparison. In addition to the better approximation error, the additional computational cost of the algorithms in this paper leads to several favorable properties that are hard to measure quantitatively: the guarantee of a $2$-approximation, sparsity, non-mass split, numerical stability, and support in $S$ and not only in $S_{\text{org}}$. For grid-structured data, the tradeoff to obtain these properties may often not be worth the additional time in practice. However, our algorithms run without specification of a fixed support set for the solution, and thus have the ability to work for any data. Next, we turn to a best-case type of data for our algorithms, where we are able to scale our computations to thousands of measures. Notably, it is data in general position, for which many algorithms in the literature do not work at all.

4.2.2 Scaling for data in general position

Applications in operations research often are based on a small set of geographical locations that do not exhibit an obvious structure. In this section, we consider data sets as depicted in Figure 7 - a set of $N$ measures that all have the same, small support of size $|P_{\text{max}}|$. This configuration leads to a best-case scenario in that $|S_{\text{org}}|=|P_{\text{max}}|$, i.e., the union of original supports is the same as any of the supports. This has an extremely positive effect on the size of the LP for Algorithm 1, respectively the first iteration of Algorithm 3: the number of variables is $|P_{\text{max}}|+N\cdot|P_{\text{max}}|^{2}$ and the number of constraints is $2(N\cdot|P_{\text{max}}|)$. In comparison, if the support points of the $P_{i}$ did not overlap, one would have $N\cdot|P_{\text{max}}|+(N\cdot|P_{\text{max}}|)^{2}$ variables and $N\cdot|P_{\text{max}}|+N^{2}\cdot|P_{\text{max}}|^{2}$ constraints. The advantage is a factor of about $N$ in the number of both variables and constraints. Note, however, that the same does not transfer to the LPs in later iterations of Algorithm 3, where $\text{supp}(\bar{P}^{\prime})$ is used as the new support, and is only guaranteed to satisfy the bound in Theorem 3.3.

In Table 3, we report on average errors and completion times for runs of Algorithm 1 and Algorithm 3. The numbers in each row were derived from $100$ runs for the given number of measures, all consisting of the same $|P_{\text{max}}|=9$ support points in general position, randomly generated for each run. (We discuss the effect of scaling to larger $|P_{\text{max}}|$ below; here we exhibit a best-case scenario.) The measures were constructed through a random assignment of varying masses to the support. We also chose the weights $\lambda_{i}$ randomly. Figure 10 (top) shows some approximate barycenters computed in this setting.

The combination of data in general position, which makes a discretization of the underlying space unavailable, and a non-uniform weight vector makes for an impractical setting for algorithms in the literature; see Section 1.1. This is why our comparisons are restricted to an exact, LP-based solution. However, exact computations in this setting are extremely hard due to the (always) exponential scaling of $S$ and the corresponding LPs. (Here, $|S|=|P_{\text{max}}|^{N}$.) The largest number of measures for which we successfully found an exact solution is $12$ (in a bit less than 3.5 hours). Because of this, approximation errors are not available for more than $12$ measures. At the same time, the LP for Algorithm 1 for such an instance still is of trivial size: it has $12\cdot 9^{2}$ variables and $2(12\cdot 9)$ constraints. The speed-up over an exact computation for a set of $8$ measures (the example in Section 4.2) is a factor of more than $600$. This factor escalates quickly - for $12$ measures our algorithm is already faster by a factor of more than $4000$.

We have been able to run Algorithm 1 for up to $20000$ measures in about 15 minutes. Instances up to $1000$ measures solve in less than $20$ seconds. One of the main reasons for these low running times is the extremely low number of constraints. The same scalability of Algorithm 3 cannot be expected, as the size of the support for LPs in later iterations is equal to $|\bar{P}^{\prime}|$, which is between linear and quadratic in $N\cdot|P_{\text{max}}|$ (Corollary 5, Theorem 3.3). We were able to run it for up to $500$ measures in a bit more than $12$ minutes. Scaling further, one quickly reaches a point where the LP for the second iteration cannot be constructed anymore. (For $1000$ measures, it would have close to $100$ million variables.)

The times reported in Table 3 include the setup of the problems. Again, we observe that $\bar{P}_{\text{org}}$ is significantly improved in the first step towards $\bar{P}^{\prime}$, and that further iterations do not change it noticably anymore. The average number of iterations increases slowly with the number of measures. The increase is slower than in the grid-structured setting. Informally, less repetition in the weighted centroids for different combinations of support points means that later iterations of Algorithm 3 are less likely to further improve the solution.

Unlike for grid-structured data, a scaling of $|P_{\text{max}}|$ for Algorithm 1 is easier in this setting. Figure 10 (bottom) show two examples with a larger support. Recall that the number of variables for Algorithm 1 is quadratic in $|P_{\text{max}}|$. Doubling the number of support points increases the number of variables in the problem by factor $4$, the same effect as increasing the number of measures by factor $4$. The number of constraints is doubled (and remains extremely low in comparison to the number of variables). An instance with $2500$ measures of $18$ support points exhibits a similar running time to an instance with $10000$ measures of $9$ support points (line 9 in Table 3).

Summing up, the combination of few, overlapping support points in general position and a large number of measures is a best-case scenario for the combination of a single run of Algorithms 1 and 2. Applications in operations research, such as facility location problems, often fall into this category. In this setting, exact computations are impossible for any reasonable problem size. Algorithms in the literature are not designed to work (or work well) for such data. In this situation, we recommend use of the presented methods for computational speed. In all other situations, the favorable properties of the output - like sparsity, non-mass split, and a guaranteed error bound - have to be crucial to the application to be worth the significant additional computational cost over popular heuristics.

Acknowledgments

The author would like to thank Ethan Anderes for the support with implementations in the Julia language, and Jacob Miller for the helpful discussions. The author gratefully acknowledges support through the Collaboration Grant for Mathematicians Polyhedral Theory in Data Analytics of the Simons Foundation.

Appendix

Appendix 0.A - Proofs of Theorems 3.3 and 3.4

We begin by proving Theorem 3.3.

Theorem 3

Algorithm 2 returns a measure $\bar{P}^{\prime}$ supported on a subset of $S$ with $\phi(\bar{P}^{\prime})\leq 2\cdot\phi(\bar{P})$ and there is a non-mass splitting transport realizing this bound. Further $|\bar{P}^{\prime}|\leq(\sum_{i=1}^{N}|P_{i}|-N+1)^{2}$.

Proof

First, note that the $P_{i}^{l}$ constructed in Step $1$ satisfy $\text{supp}(P_{i}^{l})\subset\text{supp}(P_{i})$. Thus $\text{supp}(\bar{P}^{l})\subset S$, and consequently $\text{supp}(\bar{P}^{\prime})\subset S$. Further, $\bar{P}^{\prime}=\sum_{l=1}^{r}\bar{P^{l}}$ is a measure. This holds because $\sum_{l=1}^{r}d_{l}=\sum_{l=1}^{r}z_{t_{l}}=1$, because Step $2$ does not affect this sum, and because the total mass in $\bar{P^{l}}$ equals $d_{l}$ by construction. Thus, $\bar{P}^{\prime}$ is a measure supported in $S$.

Second, we prove correctness of Step $2$. We will show that a greedily lexicographically maximal $(d_{1},\dots,d_{r})$ is created while retaining an approximate barycenter in $\text{supp}(P_{\text{org}})$. In particular, we have to show that the objective function value $\phi(\bar{P}_{\text{org}})$ does not change during the shift of mass. For a simple wording, let $\bar{P}_{\text{lex}}$ be the measure corresponding to $(d_{1},\dots,d_{r})$ after Step $2$. We will prove $\phi(\bar{P}_{\text{org}})=\phi(\bar{P}_{\text{lex}})$.

Let $x_{iq_{i}}^{l}\in P_{i}^{l}$ for $i\leq N$ and $c=\sum_{i=1}^{N}\lambda_{i}x_{iq_{i}}^{l}$, as in Step $2a)$. Then $\|c-s_{l}\|\leq\|c-s_{j}\|$ for all $j\neq l$. To see this, recall

[TABLE]

as demonstrated in the proof of Theorem 3.1. If $\|c-s_{l}\|>\|c-s_{j}\|$ for some $j\neq l$, $\bar{P}_{\text{org}}$ would not have been optimal.

By $q_{i}=\text{arg}\max_{q\leq|P_{i}^{l}|}(s_{j}-s_{l})^{T}x^{l}_{iq}$ in Step $2a)$, we pick the $x^{l}_{iq_{i}}$ such that their weighted centroid $c=\sum_{i=1}^{N}\lambda_{i}x^{l}_{iq_{i}}$ maximizes the difference $\|c-s_{l}\|^{2}-\|c-s_{j}\|^{2}\leq 0$. Only if $\|c-s_{l}\|^{2}=\|c-s_{j}\|^{2}$, mass is shifted from $s_{l}$ to $s_{j}$. But then the approximation error does not change, because

[TABLE]

Thus, the objective function value does not change during Step $2$; we have $\phi(\bar{P}_{\text{org}})=\phi(\bar{P}_{\text{lex}})$.

By definition of the running indices $l$ and $j$, mass can only be moved from support points of higher index $l$ to support points of lower index $i$. For each pair of $l$ and $j$, we repeat this shift of mass until there is no weighted centroid with $\|c-s_{l}\|=\|c-s_{j}\|$ anymore. Due to decreasing $l$ in the outer loop and increasing $j$ in the inner loop, $(d_{1},\dots,d_{r})$ is transformed to be greedily lexicographically maximal and the corresponding measure remains an approximate barycenter.

Next, we prove correctness of Steps $3$ and $4$. We show that $\phi(\bar{P}_{\text{org}})\geq\phi(\bar{P}^{\prime})$. Further, we show that for each constructed partial measure $\bar{P}^{l}$ there is a non-mass splitting transport to the $P_{i}^{l}$, and that they combine to a $\bar{P}^{\prime}$ that allows for a non-mass splitting transport that is at least as good as an optimal transport for $\bar{P}_{\text{org}}$. Finally, we show $|\bar{P}^{\prime}|\leq(\sum_{i=1}^{N}|P_{i}|-N+1)^{2}$.

Recall that in Step $3$, the mass of each $s_{l}$ is spread out to a set of weighted centroids to obtain $\bar{P}^{l}$. Independently of how the $x^{l}_{iq_{i}}$ are picked from the $P^{l}_{i}$ for all for all $i\leq N$, their weighted centroid $c=\sum_{i=1}^{N}\lambda_{i}x^{l}_{iq_{i}}$ satisfies $\sum_{i=1}^{N}\lambda_{i}\|c-x_{iq_{i}}^{l}\|^{2}\leq\sum_{i=1}^{N}\lambda_{i}\|s_{l}-x_{iq_{i}}^{l}\|^{2}$. By construction of $\bar{P}^{\prime}$ from the $\bar{P}^{l}$ (Step $4$), this already implies $\phi(\bar{P}^{\prime})\leq\phi(\bar{P}_{\text{org}})$. The algorithm started with a $2$-approximation, and thus it is guaranteed to return a $\bar{P}^{\prime}$ with $\phi(\bar{P}^{\prime})\leq 2\cdot\phi(\bar{P})$.

The existence of a non-mass splitting transport from $\bar{P}^{\prime}$ to $P_{1},\dots,P_{N}$, and the fact that this transport realizes the above bound, is a consequence of two reasons. First, each $\bar{P}^{l}$ itself allows for a non-mass splitting transport to the $P_{i}^{l}$ by lexicographically maximal choice of the $x_{iq_{i}}^{l}$ in Step $3a)$: due to this choice, the first constructed weighted centroid $c$ is lexicographically maximal among all (possible) weighted centroids that can be constructed from any $x_{iq}^{l}$ in the $P^{l}_{i}$. Further, by reducing the mass at each used support point by $d_{\text{min}}$ in Step $3b)$, at least one of the $d^{l}_{iq_{i}}$ becomes [math]. The corresponding support point is removed from $P^{l}_{i}$ (followed by some reindexing) and thus cannot be used for the construction of a weighted centroid in further iterations. Thus, the second centroid constructed in the inner loop is lexicographically strictly smaller than the first one. The same holds for all subsequent ones.

Second, any two partial measures $\bar{P}^{l_{1}}$, $\bar{P}^{l_{2}}$ from Step $3$ satisfy $\text{supp}(\bar{P}^{l_{1}})\cap\text{supp}(\bar{P}^{l_{2}})=\emptyset$ for $l_{1}\neq l_{2}$, because of the earlier preprocessing in Step $2$: weighted centroids that would be equally distant from both $s_{l_{1}}$ and $s_{l_{2}}$ cannot exist, because this would have caused a shift of mass to the lower index in Step $2$ to create a lexicographically larger $(d_{1},\dots,d_{r})$. Summing up, $\bar{P}^{\prime}$ consists of a set of distinct support points, for which it is trivial to give a non-mass splitting transport to the $P_{i}$ that is at least as good as an optimal transport for $\bar{P}_{\text{org}}$: this transport just sends the whole mass of each support point in $\bar{P}^{\prime}$ to the support points in the $P_{i}$ that were used for its construction.

The removal of at least one support point from a $P^{l}_{i}$ in Step $3b)$ implies that there are at most $\sum_{i=1}^{N}|P^{l}_{i}|-N+1$ runs of $3a)$ and $3b)$ to construct a $P^{l}$: the ’go back to $a)$’ statement is applied while $d_{l}>0$; this is the case while there still is a support point in a $P^{l}_{i}$ with mass on it. In the final run of Steps $3a)$ and $3b)$ for each $P^{l}$, all the $P^{l}_{i}$ have precisely one support point with the same mass left. This gives the claimed bound, and in particular $|P^{l}|\leq\sum_{i=1}^{N}|P^{l}_{i}|-N+1$.

Due to $|P^{l}_{i}|\leq|P_{i}|$ and $|\bar{P}_{\text{org}}|\leq\sum_{i=1}^{N}|P_{i}|-N+1$, we obtain

[TABLE]

Thus $\bar{P}^{\prime}$ satisfies all claimed properties. $\square$

Next, we prove that Algorithm 2 runs in strongly-polynomial time.

Theorem 4

For all rational input, a measure can be computed in strongly-polynomial time that is a $2$-approximation of a barycenter and for which there is a non-mass splitting transport realizing this bound.

Proof

We consider the running time of each part of the algorithm. For readability, we say ‘polynomial’ in this proof in place of ‘strongly-polynomial’. We use ‘linear’ and ‘quadratic’ to refer to the bit size $\mathcal{I}$ of the input. Note that $N$, the $|P_{i}|$, and the dimension $d$ are all bounded above by $|\mathcal{I}|$.

In Step $1$, the input for the subsequent steps is created. By sparsity of $\bar{P}_{\text{org}}$, $r\leq\sum_{i=1}^{N}|P_{i}|-N+1$. For each of the $r$ support points $s_{l}$, $N$ images $P^{l}_{i}$ with $|P^{l}_{i}|\leq|P_{i}|$ are created. In the application of the stated rule, each $y_{it_{l}k}$ has to be processed (at most) once. For each $y_{it_{l}k}$, a single comparison and a fixed number of elementary operations suffices to update the support point and mass in $P_{i}^{l}$. In total, data structures of polynomial size are created in polynomial time.

Step $2$ is the preprocessing of $(d_{1},\dots,d_{r})$ to be greedily lexicographically maximal. For each pair of support points $s_{l},s_{j}$ with $j<l$, we perform the inner part of the loop. Finding $q_{i}$ in $2a)$ can be done by considering all $x^{l}_{iq}\in P_{i}^{l}$ exactly once and comparing the inner products $(s_{j}-s_{l})^{T}x^{l}_{iq}$. This is possible in linear time. $c$ is created through the scaling and the sum of $N$ rational $d$-dimensional vectors.

Step $2b)$ begins with the computation of $c-s_{j}$ and $c-s_{l}$, then computes $\|c-s_{j}\|^{2}=(c-s_{j})^{T}(c-s_{j})$ and $\|c-s_{l}\|^{2}=(c-s_{l})^{T}(c-s_{l})$, and then compares the two values. This is possible in quadratic time. Picking the minimal mass among the $x_{iq_{i}}^{l}$ is possible in linear time, and so is updating the masses, performing the set operations on $P_{i}^{l}$ and $P_{i}^{j}$, and reindexing. In this update, $|P^{l}_{i}|$ is reduced by at least one, so the ’go back to $a)$’ statement is followed at most $|P^{l}_{i}|$ times. Summing up, Step $2$ runs in polynomial time.

Step $3$ performs the spread-out of the $r$ support points. Picking a lexicographically maximal support point $x_{iq_{i}}^{l}$ in $3a)$ can be done by considering all support points in $P_{i}^{l}$ once. One saves the current best support point and compares each other support point with respect to their lexicographic order. For identifying the lexicographic order of a pair of $d$-dimensional support points, (at most) all $d$ of their coefficients have to be compared to each other. This is possible in linear time. Again, $c$ is created through the scaling and the sum of $N$ rational $d$-dimensional vectors.

In $3b)$, we pick the minimal mass among the $x_{iq_{i}}^{l}$ used for the construction of $c$, which can be done in linear time. The same holds for the update of masses, the set operations on $P_{i}^{l}$, and the reindexing. By this update, the size of one of the $|P^{l}_{i}|$ is reduced by at least one, so the ’go back to $a)$’ statement is followed not more than $\sum_{i=1}^{N}|P^{l}_{i}|$ times; more precisely, there are at most $|P^{l}_{i}|-N+1$ runs of $3a)$ and $3b)$ for each $l$. Summing up, the construction of each $\bar{P}^{l}$ runs in polynomial time, and so does the construction of all the $\bar{P}^{l}$.

In Step $4$, the partial measures $\bar{P}^{l}$ are combined to obtain $\bar{P}^{\prime}$. This is the construction of a measure with the appropriate mass put on at most $|\bar{P}^{\prime}|\leq(\sum_{i=1}^{N}|P_{i}|-N+1)^{2}$ support points. Each of these support points is just a copy of a support point in one of the $\bar{P}^{l}$. Thus, all steps run in polynomial time, which proves the claim. $\square$

Appendix 0.B - Proof of Theorem 3.5

Theorem 5

Algorithm 3 returns an approximate barycenter $\bar{P}^{\prime}$ supported on a subset of $S$ for which $\phi(\bar{P}^{\prime})\leq 2\cdot\phi(\bar{P})$, where $\bar{P}$ is a barycenter, and there is a non-mass splitting optimal transport realizing this bound. Further $|\bar{P}^{\prime}|\leq\sum_{i=1}^{N}|P_{i}|-N+1$.

Proof

First, recall that the output $\bar{P}^{\prime}$ of Algorithm 2 (Step $2$) always satisfies $\text{supp}(\bar{P}^{\prime})\subset S$. Further, Algorithm 2 always returns a measure that has a corresponding non-mass splitting transport. As $\bar{P}_{\text{org}}$ from Algorithm 1 (Step $1$) is not changed in the final run of Algorithm 2, the returned non-mass splitting transport is optimal. Further, recall that all approximate barycenters $\bar{P}_{\text{org}}$ computed in Step $1$ have a support that satisfies $|\bar{P}_{\text{org}}|\leq\sum_{i=1}^{N}|P_{i}|-N+1$. This transfers to the sparsity of $\bar{P}^{\prime}$ returned by Algorithm 3.

It remains to prove termination of Algorithm 3 and the error bound. We will do so by showing that $\phi(\bar{P}^{\prime})<\phi(\bar{P}_{\text{org}})$ if $\bar{P}^{\prime}\neq\bar{P}_{\text{org}}$ for $\bar{P}_{\text{org}},\bar{P}^{\prime}$ from the same iteration. This leads to a strictly decreasing sequence of values $\phi(\bar{P}^{\prime})$ as long as the algorithm keeps running. The first approximate barycenter in this sequence already is a $2$-approximation and it can only become better. This immediately gives $\phi(\bar{P}^{\prime})\leq 2\cdot\phi(\bar{P})$. At the end of each Step $2$, we update $S_{\text{org}}=\text{supp}(\bar{P}^{\prime})\subset S$ before going back to Step $1$, where an exact optimum over this new support, a subset of $S$, is computed. Because of this, and the fact that there are only finitely many subsets of $S$, the sequence of values $\phi(\bar{P}^{\prime})$ is finite.

Now, it only remains to prove that $\phi(\bar{P}^{\prime})<\phi(\bar{P}_{\text{org}})$ if $\bar{P}^{\prime}\neq\bar{P}_{\text{org}}$. We begin by considering Step $3$ of Algorithm 2. Assume $P^{l}_{i}$ consists of a single support point $x^{l}_{i1}$ for all $i\leq N$. Then the unique barycenter $\bar{P}^{l}$ of the $P^{l}_{i}$ is the weighted centroid $c=\sum_{i=1}^{N}\lambda_{i}x^{l}_{i1}$ and the cost of transport from $\bar{P}^{l}$ to all the $P^{l}_{i}$ is $\phi(\bar{P}^{l})=d_{l}\cdot\sum_{i=1}^{N}\lambda_{i}\|c-x^{l}_{i1}\|^{2}$. For all $s\neq c$, in particular for $s=s_{l}$, we get

[TABLE]

If some of the $P^{l}_{i}$ consist of more than one support point, Step $3$ selects a set of exactly one support point $x_{iq}^{l}$ from each measure $P^{l}_{i}$, forms a weighted centroid $c$ with corresponding mass $d_{c}=d_{\text{min}}$, and adds it to $\text{supp}(\bar{P}^{l})$. Then this scheme is repeated for the remaining support points and remaining mass. This means that $\bar{P}^{l}$ is constructed as a set of weighted centroids $c$ of support points $x^{l}_{iq}$ to which these centroids $c$ transport. Each of them satisfies $d_{c}\cdot\sum_{i=1}^{N}\lambda_{i}\|c-x^{l}_{iq}\|^{2}\leq d_{c}\cdot\sum_{i=1}^{N}\lambda_{i}\|s_{l}-x^{l}_{iq}\|^{2}$. By summing over all $c$ that are constructed, one obtains

[TABLE]

Informally, it is at least as costly to transport to the measures $P^{l}_{i}$ from the support point $s_{l}$ as from the set of weighted centroids (with appropriate masses) constituting $\bar{P}^{l}$. Equality in the above can only hold if the single support point $s_{l}$ itself already is the weighted centroid of single-support point measures $P^{l}_{1},\dots,P^{l}_{N}$. But this means that Step $3$ of Algorithm $2$ just copies $s_{l}$ with mass $d_{l}$ to $\bar{P}^{l}$. The algorithm stops when $\bar{P}^{\prime}=\bar{P}_{\text{org}}$. By $\phi(\bar{P}^{\prime})=\sum_{l=1}^{r}\phi(\bar{P}^{l})$, this means all $s_{l}$ have to satisfy $\phi(\bar{P}^{l})=\sum_{i=1}^{N}\lambda_{i}\sum_{q=1}^{|P_{i}^{l}|}d^{l}_{iq}\cdot\|s_{l}-x_{iq}^{l}\|^{2}.$ So all $s_{l}$ are already the weighted centroids of their single-support measures $P^{l}_{i}$.

Further, note that when a shift of mass from $s_{l}$ to $s_{j}$ with $j<l$ happens in Step $2$ of Algorithm 2, then Step $3$ is guaranteed to find a strictly better transport than before: there exists a set of support points that, before the shift, receive transport from $s_{l}$, but have a weighted centroid $c\neq s_{l}$. Such a set of support points would be moved from $P^{l}_{i}$ to $P^{j}_{i}$ (and at least one of the support points was not associated to $s_{j}$ before). Then $s_{j}$ is guaranteed to split mass and, in the following Step $3$, the cost of transport is strictly improved; see above.

Thus $\phi(\bar{P}^{\prime})<\phi(\bar{P}_{\text{org}})$ if $\bar{P}^{\prime}\neq\bar{P}_{\text{org}}$ and Algorithm 3 terminates with $\bar{P}^{\prime}=\bar{P}_{\text{org}}$ in the final iteration. $\square$