Semi-discrete optimal transport: a solution procedure for the unsquared Euclidean distance case
Abstract
EconStor is a publication server for scholarly economic literature, provided as a non-commercial public service by the ZBW.
Full text
Hartmann, Valentin; Schuhmacher, Dominic Article — Published Version Semi-discrete optimal transport: a solution procedure for the unsquared Euclidean distance case Mathematical Methods of Operations Research Provided in Cooperation with: Springer Nature Suggested Citation: Hartmann, Valentin; Schuhmacher, Dominic (2020) : Semi-discrete optimal transport: a solution procedure for the unsquared Euclidean distance case, Mathematical Methods of Operations Research, ISSN 1432-5217, Springer, Berlin, Heidelberg, Vol. 92, Iss. 1, pp. 133-163, https://doi.org/10.1007/s00186-020-00703-z This Version is available at: https://hdl.handle.net/10419/288281 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/
Mathematical Methods of Operations Research (2020) 92:133–163 https://doi.org/10.1007/s00186-020-00703-z ORIGINAL ARTICLE Semi-discrete optimal transport: a solution procedure for the unsquared Euclidean distance case Valentin Hartmann1,2 ·Dominic Schuhmacher1 Received: 22 May 2019 / Revised: 20 December 2019 / Published online: 12 February 2020 © The Author(s) 2020 Abstract We consider the problem of finding an optimal transport plan between an absolutely continuous measure and a finitely supported measure of the same total mass when the transport cost is the unsquared Euclidean distance. We may think of this problem as closest distance allocation of some resource continuously distributed over Euclidean space to a finite number of processing sites with capacity constraints. This article gives a detailed discussion of the problem, including a comparison with the much better studied case of squared Euclidean cost. We present an algorithm for computing the optimal transport plan, which is similar to the approach for the squared Euclidean cost by Aurenhammer et al. (Algorithmica 20(1):61–76, 1998) and Mérigot (Comput Graph Forum 30(5):1583–1592, 2011). We show the necessary results to make the approach work for the Euclidean cost, evaluate its performance on a set of test cases, and give a number of applications. The later include goodness-of-fit partitions, a novel visual tool for assessing whether a finite sample is consistent with a posited probability density. Keywords Monge–Kantorovich problem ·Spatial resource allocation ·Wasserstein metric ·Weighted Voronoi tessellation Mathematics Subject Classification Primary 65D18; Secondary 51N20 ·62-09 VH was partially supported by Deutsche Forschungsgemeinschaft RTG 2088. We thank Marcel Klatt for helpful discussions and three anonymous referees for comments that led to an improvement of the paper. BDominic Schuhmacher [email protected] Valentin Hartmann v[email protected] 1Institute for Mathematical Stochastics, University of Goettingen, Goldschmidtstr. 7, 37077 Goettingen, Germany 2Present Address: IC IINFCOM DLAB, EPFL, Station 14, 1015 Lausanne, Switzerland 123
134 V. Hartmann, D. Schuhmacher 1 Introduction Optimal transport and Wasserstein metrics are nowadays among the major tools for analyzing complex data. Theoretical advances in the last decades characterize existence, uniqueness, representation and smoothness properties of optimal transport plans in a variety of different settings. Recent algorithmic advances (Peyré and Cuturi 2018) make it possible to compute exact transport plans and Wasserstein distances between discrete measures on regular grids of tens of thousands of support points, see e.g. Schmitzer (2016, Sect. 6), and to approximate such distances (to some extent) on larger and/or irregular structures, see Altschuler et al. (2017) and references therein. The development of new methodology for data analysis based on optimal transport is a booming research topic in statistics and machine learning, see e.g. Sommerfeld and Munk (2018), Schmitz et al. (2018), Arjovsky et al. (2017), Genevay et al. (2018), and Flamary et al. (2018). Applications are abundant throughout all of the applied sciences, including biomedical sciences (e.g. microscopy or tomography images; Basua et al. 2014,Gramfortetal.2015), geography (e.g. remote sensing; Courty et al. 2016, Guo et al. 2017), and computer science (e.g. image processing and computer graphics; Nicolas 2016, Solomon et al. 2015). In brief: whenever data of a sufficiently complex structure that can be thought of as a mass distribution is available, optimal transport offers an effective, intuitively reasonable and robust tool for analysis. More formally, for measures μand νon Rdwith μ(Rd)=ν(Rd)<∞the Wasserstein distance of order p≥1 is defined as Wp(μ, ν) =min πRd×Rdx−ypπ(dx,dy)1/p ,(1) where the minimum is taken over all transport plans (couplings) πbetween μand ν, i.e. measures πon Rd×Rdwith marginals π(A×Rd)=μ(A)and π(Rd×A)=ν(A) for every Borel set A⊂Rd. The minimum exists by Villani (2009, Theorem 4.1) and it is readily verified, see e.g. Villani (2009, after Example 6.3), that the map Wpis a [0,∞]-valued metric on the space of measures with fixed finite mass. The constraint linear minimization problem (1) is known as Monge–Kantorovich problem (Kantorovich 1942; Villani 2009). From an intuitive point of view, a minimizing π describes how the mass of μis to be associated with the mass of νin order to make the overall transport cost minimal. Atransport map from μto νis a measurable map T:Rd→Rdsatisfying T#μ=ν, where T#denotes the push-forward, i.e. (T#μ)(A)=μ(T−1(A)) for every Borel set A⊂Rd. We say that T induces the coupling π=πTif πT(A×B)=μ(A∩T−1(B)) 123
Semi-discrete optimal transport: the unsquared Euclidean distance case 135 for all Borel sets A,B⊂Rd, and call the coupling πdeterministic in that case. It is easily seen that the support of πTis contained in the graph of T. Intuitively speaking, we associate with each location in the domain of the measure μexactly one location in the domain of the measure νto which positive mass is moved, i.e. the mass of μis not split. The generally more difficult (non-linear) problem of finding (the p-th root of) inf TRdx−T(x)pμ(dx)=inf TRd×Rdx−ypπT(dx,dy), (2) where the infima are taken over all transport maps Tfrom μto ν(and are in general not attained) is known as Monge’s problem (Monge 1781; Villani 2009). In practical applications, based on discrete measurement and/or storage procedures, we often face discrete measures μ=m i=1μiδxiand ν=n j=1νjδyj, where {x1,...,xm},{y1,...,yn}are finite collections of support points, e.g. grids of pixel centers in a grayscale image. The Monge–Kantorovich problem (1) is then simply the discrete transport problem from classical linear programming (Luenberger and Ye 2008): Wp(μ, ν) =min (πij) m i=1 n j=1 dijπij1/p ,(3) where dij =xi−yjpand any measure π=m i=1n j=1πijδ(xi,yj)is represented by the m×nmatrix (πij)i,jwith nonnegative entries πij satisfying n j=1 πij =μifor 1 ≤i≤mand m i=1 πij =νjfor 1 ≤j≤n. Due to the sheer size of mand nin typical applications this is still computationally a very challenging problem; we have e.g. m=n=106for 1000 ×1000 grayscale images, which is far beyond the performance of a standard transportation simplex or primal-dual algorithm. Recently many dedicated algorithms have been developed, such as (Schmitzer 2016), which can give enormous speed-ups mainly if p=2 and can compute exact solutions for discrete transportation problems with 105support points in seconds to a few minutes, but still cannot deal with 106or more points. Approximative solutions can be computed for this order of magnitude and p=2 by variants of the celebrated Sinkhorn algorithm (Cuturi 2013; Schmitzer 2019; Altschuler et al. 2017), but it has been observed that these approximations have their limitations (Schmitzer 2019; Klatt et al. 2019). The main advantage of using p=2 is that we can decompose the cost function as x−y2=x2+y2−2xyand hence formulate the Monge–Kantorovich problem equivalently as maxπRd×Rdxyπ(dx,dy). For the discrete problem (3)this decomposition is used in Schmitzer (2016) to construct particularly simple so-called shielding neighborhoods. But also if one or both of μand νare assumed absolutely continuous with respect to Lebesgue measure, this decomposition for p=2 has clear computational advantages. For example if the measures μand νareassumedtohave 123
136 V. Hartmann, D. Schuhmacher densities fand g, respectively, the celebrated Brenier’s theorem, which yields an optimal transport map that is the gradient of a convex function u(McCann 1995), allows to solve Monge’s problem by finding a numerical solution uto the Monge-Ampère equation det(D2u(x)) =f(x)g(∇u(x)); see Santambrogio (2015, Sect. 6.3) and the references given there. In the rest of this article we focus on the semi-discrete setting, meaning here that the measure μis absolutely continuous with respect to Lebesgue measure and the measure νhas finite support. This terminology was recently used in Wolansky (2015), Kitagawa et al. (2019), Genevay et al. (2016) and Bourne et al. (2018) among others. In the semi-discrete setting we can represent a solution to Monge’s problem as a partition of Rd, where each cell is the pre-image of a support point of νunder the optimal transport map. We refer to such a partition as optimal transport partition. In the case p=2 this setting is well studied. It was shown in Aurenhammer et al. (1998) that an optimal transport partition always exists, is essentially unique, and takes the form of a Laguerre tessellation, a.k.a. power diagram. The authors proved further that the right tessellation can be found numerically by solving a (typically high dimensional) unconstrained convex optimization problem. Since Laguerre tessellations are composed of convex polytopes, the evaluation of the objective function can be done very precisely and efficiently. Mérigot (2011) elaborates details of this algorithm and combines it with a powerful multiscale idea. In Kitagawa et al. (2019) a damped Newton algorithm is presented for the same objective function and the authors are able to show convergence with optimal rates. In this article we present the corresponding theory for the case p=1. It is shown in Sect. 2.3 of Crippa et al. (2009) and independently in Geiß et al. (2013), which both treat more general cost functions, that an optimal transport partition always exists, is essentially unique and takes the form of a weighted Voronoi tessellation, or more precisely an Apollonius diagram. We extend this result somewhat within the case p=1 in Theorems 1and 2below. We prove then in Theorem 3that the right tessellation can be found by optimizing an objective function corresponding to that from the case p=2. Since the cell boundaries in an Apollonius diagram in 2d are segments of hyperbolas, computations are more involved and we use a new strategy for computing integrals over cells and for performing line search in the optimization method. Details of the algorithm are given in Sect. 4and the complete implementation can be downloaded from Github1and is included in the latest version of the transport-package (Schuhmacher et al. 2019) for the statistical computing environment R(R Core Team 2017). Up to Sect. 4the present paper is a condensed version of the thesis (Hartmann 2016), to which we refer from time to time for more details. In the remainder we evaluate the performance of our algorithm on a set of test cases (Sect. 5), give a number of applications (Sect. 6), and provide a discussion and open questions for further research (Sect. 7). At the time of finishing the present paper, it has come to our attention that Theorem 2.1 of Kitagawa et al. (2019), which is for very general cost funtions including the Euclidean distance (although the remainder of the paper is not), has a rather large overlap with our Theorem 3. Within the case of Euclidean cost it assumes somewhat 1https://github.com/valentin-hartmann-research/semi-discrete-transport. 123
Semi-discrete optimal transport: the unsquared Euclidean distance case 137 stronger conditons than our Theorem 3, namely a compact domain Xand a bounded density for μ. In addition the statement is somewhat weaker as it does not contain our statement (c). We also believe that due to the simpler setting of p=1 our proof is accessible to a wider audience and it is more clearly visible that the additional restrictions on Xand μare in fact not needed. We end this introduction by providing some first motivation for studying the semidiscrete setting for p=1. This will be further substantiated in the application Sect. 6. 1.1 Why semi-discrete? The semi-discrete setting appears naturally in problems of allocating a continuously distributed resource to a finite number of sites. Suppose for example that a fast-food chain introduces a home delivery service. Based on a density map of expected orders (the “resource”), the management would like to establish delivery zones for each branch (the “sites”). We assume that each branch has a fixed capacity (at least in the short run), that the overall capacity matches the total number of orders (peak time scenario), and that the branches are not too densely distributed, so that the Euclidean distance is actually a reasonable approximation to the actual travel distance; see Boscoe et al. (2012). We take up this example in Sect. 6.2. A somewhat different model that adds waiting time costs to the distance-based costs instead of using capacity constraints was studied theoretically in Crippa et al. (2009). An important general class that builds on resource allocation are location-alloca- tion problems: where to position a number of sites (branches, service stations, etc.) in such a way that the sum of the resource allocation cost plus maybe further costs for installation, maintainance and waiting times is minimized, possibly under capacity and/or further constraints. See e.g. Mallozzi et al. (2019) for a flexible model, which was algorithmically solved via discretizing the continuous domain. Positioning of sites can also be competitive, involving different agents (firms), such as in Núñez and Scarsini (2016). A special case of location-allocation is the quantization problem, which consists in finding positions and capacities of sites that minimize the resulting resource allocation cost. See Bourne et al. (2018, Sect. 4) for a recent discussion using incomplete transport and p=2. As a further application we propose in Sect. 6.3 optimal transport partitions as a simple visual tool for investigating local deviations from a continuous probability distribution based on a finite sample. Since the computation of the semi-discrete optimal transport is linear in the resolution at which we consider the continuous measure (for computational purposes), it can also be attractive to use the semi-discrete setting as an approximation of either the fully continuous setting (if νis sufficiently simple) or the fully discrete setting (if μ has a large number of support points). This will be further discussed in Sect. 2. 123
138 V. Hartmann, D. Schuhmacher 1.2 Why p=1? The following discussion highlights some of the strengths of optimal transport based on an unsquared Euclidean distance (p=1), especially in the semi-discrete setting, and contrasts p=1 with p=2. From a computational point of view the case p=2 can often be treated more efficiently, mainly due to the earlier mentioned decomposability, leading e.g. to the algorithms in Schmitzer (2016) in the discrete and Aurenhammer et al. (1998), Mérigot (2011) in the semi-discrete setting. The case p=1 has the advantage that the Monge–Kantorovich problem has a particularly simple dual (Villani 2009, Particular Case 5.16), which is equivalent to Beckmann’s problem (Beckmann 1952; Santambrogio 2015, Theorem 4.6). If we discretize the measures (if necessary) to a common mesh of npoints, the latter is an optimization problem in nvariables rather than the n2variables needed for the general discrete transport formulation (3). Algorithms that make use of this reduction have been described in Solomon et al. (2014) (for general discrete surfaces) and in Schmitzer and Wirth (2019, Sect. 4) (for general incomplete transport), but their performance in a standard situation, e.g. complete optimal transport on a regular grid in Rd, remains unclear. In particular we are not aware of any performance comparisons between p=1 and p=2. In the present paper we do not make use of this reduction, but keep the source measure μtruly continuous except for an integral approximation that we perform for numerical purposes. We describe an algorithm for the semi-discrete problem with p=1 that is reasonably fast, but cannot quite reach the performance of the algorithm for p=2 in Mérigot (2011). This is again mainly due to the nice decomposition property of the cost function for p=2 or, more blatantly, the fact that we minimize for p=2 over partitions formed by line rather than hyperbola segments. From an intuitive point of view p=1 and p=2 have both nice interpretations and depending on the application setting either the one or the other may be more justified. The difference is between thinking in terms of transportation logistics or in terms of fluid mechanics. If p=1, the optimal transport plan minimizes the cumulative distance by which mass is transported. This is (up to a factor that would not change the transport plan) the natural cost in the absence of fixed costs or any other savings on long-distance transportation. If p=2, the optimal transport plan is determined by a pressureless potential flow from μto νas seen from the kinetic energy minimization formulation of Benamou and Brenier (2000), Villani (2009, Chapter 7). The different behaviors in the two cases can be illustrated by the discrete toy example in Fig. 1. Each point along the incomplete circle denotes the location of one unit of mass of μ(blue x-points) and/or ν(red o-points). The unique solution for p=1 moves one unit of mass from one end of the circular structur to the other. This is how we would go about carrying boxes around to get from the blue scenario to the red scenario. The unique solution for p=2 on the other hand is to transport each unit a tiny bit further to the next one, corresponding to a (discretized) flow along the circle. It is straightforward to adapt this toy example for the semi-discrete or the continuous setting. A more complex semi-discrete example is given in Sect. 6.1. 123
Semi-discrete optimal transport: the unsquared Euclidean distance case 139 Fig. 1 Optimal transport maps from blue-x to red-o measure with unit mass at each point. Left: the transportation logistics solution (p=1); right: the fluid mechanics solution (p=2). (Color figure online) One argument in favour of the metric W1is its nice invariant properties that are not shared by the other Wp. In particular, considering finite measures μ, ν, α on Rd satisfying μ(Rd)=ν(Rd),p≥1 and c>0, we have W1(α +μ, α +ν) =W1(μ, ν), (4) W1(cμ, cν) =cW 1(μ, ν). (5) The first result is in general not true for any other p, the second result holds with a factor c1/pon the right hand side. We prove these statements in the appendix. These invariance properties have important implications for image analysis, where it is quite common to adjust for differring levels of brightness (in grayscale images) by affine transformations. While the above equalities show that it is safe to do so for p=1, it may change the resulting Wasserstein distance and the optimal transport plan dramatically for other p; see Appendix and Sect. 6.1. It is sometimes considered problematic that optimal transport plans for p=1are in general not unique. But this is not so in the semi-discrete case, as we will see in Sect. 2: the minimal transport cost in (1) is realized by a unique coupling π, which is always deterministic. The same is true for p=2. A major difference in the case p=1 is that for d>1 each cell of the optimal transport partition contains the support point of the target measure νthat it assigns its mass to. This can be seen as a consequence of cyclical monotonicity (Villani 2009, beginning of Chapter 8). In contrast, for p=2, optimal transport cells can be separated by many other cells from their support points, which can make the resulting partition hard to interpret without drawing corresponding arrows for the assignment; see the bottom panels of Fig. 5. For this reason we prefer to use p=1 for the goodness-of-fit partitions considered in Sect. 6.3. 2 Semi-discrete optimal transport We first concretize the semi-discrete setting and introduce some additional notation. Let now Xand Ybe Borel subsets of Rdand let μand νbe probability measures on 123
140 V. Hartmann, D. Schuhmacher Xand Y, respectively. This is just for notational convenience and does not change the set of admissible measures in an essential way: we may always set X=Y=Rd and any statement about μand νwe make can be easily recovered for cμand cνfor arbitrary c>0. For the rest of the article it is tacitly assumed that d≥2 to avoid certain pathologies of the one-dimensional case that would lead to a somewhat tedious distinction of cases in various results for a case that is well-understood anyway. Moreover, we always require μto be absolutely continuous with density with respect to d-dimensional Lebesgue measure Lebdand to satisfy Xxμ(dx)<∞.(6) We assume further that ν=n j=1νjδyj, where n∈N,y1,...,yn∈Yand ν1,...ν n∈(0,1]. Condition (6) guarantees that W1(μ, ν) ≤Xxμ(dx)+Yyν(dy)=: C<∞,(7) which simplifies certain arguments. The set of Borel subsets of Xis denoted by BX. Lebesgue mass is denoted by absolute value bars, i.e. |A|=Lebd(A)for every A∈BX. We call a partition C=(Cj)1≤j≤nof Xinto Borel sets satisfying μ(Cj)=νjfor every jatransport partition from μto ν. Any such partition characterizes a transport map Tfrom μto ν, where we set TC(x)=n j=1yj1{x∈Cj}for a given transport partition C=(Cj)1≤j≤nand CT=(T−1(yj))1≤j≤nfor a given transport map T. Monge’s problem for p=1 can then be equivalently formulated as finding inf CXx−TC(x)μ(dx)=inf C n j=1Cjx−yjμ(dx), (8) where the infima are taken over all transport partitions C=(Cj)1≤j≤nfrom μto ν. Contrary to the difficulties encountered for more general measures μand νwhen considering Monge’s problem with Euclidean costs, we can give a clear-cut existence and uniqueness theorem in the semi-discrete case, without any further restrictions. Theorem 1 In the semi-discrete setting with Euclidean costs (always including d ≥2 and (6))thereisaμ-a.e. unique solution T∗to Monge’s problem. The induced coupling πT∗is the unique solution to the Monge–Kantorovich problem, yielding W1(μ, ν) =Xx−T∗(x)μ(dx). (9) Proof The part concerning Monge’s problem is a consequence of the concrete construction in Sect. 3; see Theorem 2. 123
Semi-discrete optimal transport: the unsquared Euclidean distance case 147 4 The algorithm The previous section provides the theory needed to compute the optimal transport partition. It is sufficient to find a vector w∗where is locally optimal. By convexity, w∗is then a global minimizer of and Remark 1identifies the μ-a.e. unique optimal transport partition as (Vor w∗(j))1≤j≤n. For the optimization process we can choose from a variety of methods thanks to knowing the gradient ∇of analytically from Theorem 3. We consider iterative methods that start at an initial weight vector w(0)and apply steps of the form w(k+1)=w(k)+tkw(k),k≥0, where w(k)denotes the search direction and tkthe step size. Newton’s method would use w(k)=− D2(w(k))−1∇(w(k)), but the Hessian matrix D2(w(k))is not available to us. We therefore use a quasi-Newton method that makes use of the gradient. Just like Mérigot (2011) for the case p=2, we have obtained many good results using L-BFGS (Nocedal 1980), the limited-memory variant of the Broyden–Fletcher–Goldfarb–Shanno algorithm, which uses the value of the gradient at the current as well as at preceding steps for approximating the Hessian. The limited-memory variant works without storing the whole Hessian of size n×n, which is important since in applications our nis typically large. To determine a suitable step size tkfor L-BFGS, we use the Armijo rule (Armijo 1966), which has proven to be well suited for our problem. It considers different values for tkuntil it arrives at one that sufficiently decreases (w(k)): the step size tkneeds to fulfill (w(k)+tkw(k))≤(w(k))+ct k∇(w(k))Tw(k)for a small fixed cwith 0<c<1. We use the default value c=10−4of the L-BFGS library (Okazaki and Nocedal 2010) employed by our implementation, which is also given as an example by Nocedal and Wright (1999). An alternative that could be investigated is to use a non-monotone line search such as the one proposed in Grippo et al. (1986). There the above condition is relaxed by admitting a step whenever it sufficiently decreases a function value from one of the previous Kiterations, for some K≥1. This might lead to fewer function evaluations and also to convergence in fewer steps. We also considered replacing the Armijo rule with the strong Wolfe conditions (1969,1971) as done in Mérigot (2011), which contain an additional decrease requirement on the gradient. In our case, however, this requirement could often not be fulfilled because of the pixel splitting method used for computing the gradient (cf. Sect. 4.2), which made it less suited. 4.1 Multiscale approach to determine starting value To find a good starting value w(0)we use a multiscale method similar to the one proposed in Mérigot (2011). We first create a decomposition of ν, i.e. a sequence ν=ν(0),...,ν(L)of measures with decreasing cardinality of the support. Here ν(l) is obtained as a coarsening of ν(l−1)by merging the masses of several points into one point. 123
148 V. Hartmann, D. Schuhmacher It seems intuitively reasonable to choose ν(l)in such a way that W1(ν(l),ν(l−1)) is as small as possible, since the latter bounds |W1(μ, ν(l))−W1(μ, ν(l−1))|.This comes down to a capacitated location-allocation problem, which is NP-hard even in the one-dimensional case; see Sherali and Nordai (1988). Out of speed concerns and since we only need a reasonably good starting value for our algorithm, we decided to content ourselves with the same weighted K-means clustering algorithm used by Mérigot (2011) (referred to as Lloyd’s algorithm), which iteratively improves an initial aggregation of the support of ν(l−1)into |supp(ν(l))|clusters towards local optimality with respect to the squared Euclidean distance. The resulting ν(l)is then the discrete measure with the cluster centers as its support points and as weights the summed up weights of the points of ν(l−1)contained in each cluster; see Algorithm 3 in Hartmann (2016). The corresponding weighted K-median clustering algorithm, based on alternating between assignment of points to clusters and recomputation of cluster centers as the median of all weighted points in the cluster, should intuitively give a ν(l)based on which we obtain a better starting solution. This may sometimes compensate for the much longer time needed for performing K-median clustering. Having created the decomposition ν=ν(0),...,ν(L), we minimize along the sequence of these coarsened measures, beginning at ν(L)with the initial weight vector w(L,0)=0∈R|supp(ν(L))|and computing the optimal weight vector w(L,∗)for the transport from μto ν(L). Every time we pass to a finer measure ν(l−1)from the coarser measure ν(l), we generate the initial weight vector w(l−1,0)from the last optimal weight vector w(l,∗)by assigning the weight of each support point of ν(l)to all the support points of ν(l−1)from whose merging the point of ν(l)originated; see also Algorithm 2 in Hartmann (2016). 4.2 Numerical computation of 8and ∇8 For practical computation we assume here that Xis a bounded rectangle in R2and that the density of the measure μis of the form (x)= i∈I ai1Qi(x) for x∈X, where we assume that Iis a finite index set and (Qi)i∈Iis a partition of the domain Xinto (small) squares, called pixels, of equal side length. This is natural if is given as a grayscale image and we would then typically index the pixels Qi by their centers i∈I⊂Z2. It may also serve as an approximation for arbitrary . It is however easy enough to adapt the following considerations to more general (not necessarily congruent) tiles and to obtain better approximations if the function is specified more generally than piecewise constant. The optimization procedure requires the non-trivial evaluation of at a given weight vector w. This includes the integration over Voronoi cells and therefore the construction of a weighted Voronoi diagram. The latter task is solved by the package 2D Apollonius Graphs as part of the Computational Geometry Algorithms Library (CGAL 2015). The integrals we need to compute are 123
Semi-discrete optimal transport: the unsquared Euclidean distance case 149 Vor w(j) ρ(x)dx and Vor w(j)x−yjρ(x)dx. By definition the boundary of a Voronoi cell Vorw(j)is made up of hyperbola segments, each between yjand one of the other support points of ν. The integration could be performed by drawing lines from yjto the end points of those segments and integrating over the resulting triangle-shaped areas separately. This would be executed by applying an affinely-linear transformation that moves the hyperbola segment onto the hyperbola y=1/xto both the area and the function we want to integrate. The required transformation can be found in Hartmann (2016, Sect. 5.6). However, we take a somewhat more crude, but also more efficient path here, because it is a quite time-consuming task to decide which pixels intersect which weighted Voronoi cells and then to compute the (areas of the) intersections. We therefore approximate the intersections by splitting the pixels into a quadratic number of subpixels (unless the former are already very small) and assuming that each of them is completely contained in the Voronoi cell in which its center lies. This reduces the problem from computing intersections to determining the corresponding cell for each center, which the data structure used for storing the Voronoi diagram enables us to do in roughly O(log n)time; see Karavelas and Yvinec (2002). The operation can be performed even more efficiently: when considering a subpixel other than the very first one, we already know the cell that one of the neighboring subpixel’s center belongs to. Hence, we can begin our search at this cell, which is either already the cell we are looking for or lies very close to it. The downside of this approximation is that it can make the L-BFGS algorithm follow search directions along which the value of cannot be sufficiently decreased even though there exist different directions that allow a decrease. This usually only happens near a minimizing weight vector and can therefore be controlled by choosing a not too strict stopping criterion for a given subpixel resolution, see the next subsection. 4.3 Our implementation Implementing the algorithm described in this section requires two technical choices: the number of subpixels every pixel is being split into and the stopping criterion for the minimization of . We found that choosing the number of subpixels to be the smallest square number such that their total number is larger than or equal to 1000n gives a good compromise between performance and precision. The stopping criterion is implemented as follows: we terminate the optimization process once ∇(w)1/2≤εfor some ε>0. Due to Theorem 3b) this criterion yields an intuitive interpretation: ∇(w)1/2 is the amount of mass that is being mistransported, i.e. the total amount of mass missing or being in surplus at some νlocation yiwhen transporting according to the current tessellation. In our experience this mass is typically rather proportionally distributed among the different cells and tends to be assigned in a close neighborhood of the correct cell rather than far away. So even with somewhat larger ε, the computed Wasserstein distance and the overall visual impression of the optimal transport partition remain mostly the same. In the numerical examples in Sects. 5and 6we choose the value of ε=0.05. 123
150 V. Hartmann, D. Schuhmacher Fig. 3 Realizations of the measure μfor all six parameter combinations in Sect. 5. First row: smoothness s=0.5; second row: smoothness s=2.5. The correlation scale γis 0.05, 0.15 and 0.5 (from left to right) We implemented the algorithm in C++ and make it available on GitHub2under the MIT license. Our implementation uses libLBFGS (Okazaki and Nocedal 2010)forthe L-BFGS procedure and the geometry library CGAL (CGAL 2015) for the construction and querying of weighted Voronoi tessellations. The repository also contains a Matlab script to visualize such tessellations. Our implementation is also included in the latest version of the transport-package (Schuhmacher et al. 2019) for the statistical computing environment R(R Core Team 2017). 5 Performance evaluation We evaluated the performance of our algorithm by randomly generating measures μ and νwith varying features and computing the optimal transport partitions between them. The measure μwas generated by simulating its density as a Gaussian random field with Matérn covariance function on the rectangle [0,1]×[0,0.75], applying a quadratic function and normalizing the result to a probability density. Corresponding images were produced at resolution 256 ×196 pixels and were further divided into 25 subpixels each to compute integrals over Voronoi cells. In addition to a variance parameter, which we kept fixed, the Matérn covariance function has parameters for the scale γof the correlations, which we varied among 0.05, 0.15 and 0.5, and the smoothness sof the generated surface, which we varied between 0.5 and 2.5 corresponding to a continuous surface and a C2-surface, respectively. The simulation mechanism is similar to the ones for classes 2–5 in the benchmark DOTmark proposed in Schrieber et al. (2017), but allows to investigate the influence of individual parameters more directly. Figure 3shows one realization for each parameter combination. For the performance evaluation we generated 10 realizations each. 2https://github.com/valentin-hartmann-research/semi-discrete-transport. 123
Semi-discrete optimal transport: the unsquared Euclidean distance case 151 (a) n=250,s=0.5(b) n=250,s=2.5 (c) n=1000,s=0.5(d) n=1000,s=2.5 Fig. 4 Runtimes of the experiments of Sect. 5. Bars and lines indicate means and standard deviations over 200 experiments, combining 10 realizations of μwith 20 realizations of ν. The measures μare based on Gaussian random fields with Matérn covariance function; see Fig. 3. The measures νare based on support points picked uniformly at random with unit masses (blue) or masses picked from the corresponding μ (red). Rows: νwith n=250 versus n=1000 support points. Columns: smoothness parameter s=0.5 versus s=2.5. Note the different scaling. (Color figure online) The measures νhave nsupport points generated uniformly at random on [0,1]× [0,0.75],whereweusedn=250 and n=1000. We then assigned either mass 1 or mass (x)to each point xand normalized to obtain probability measures. We generated 20 independent ν-measures of the first kind (unit mass) and computed the optimal transport from each of the 10×6μ-measures for each of the 6 parameter combinations. We further generated for each of the 10 ×6μ-measures 20 corresponding ν-measures of the second kind (masses from μ) and computed again the corresponding optimal transports. The stopping criterion for the optimization was an amount of ≤0.05 of mistransported mass. The results for n=250 support points of νare shown in Fig. 4a, b, those for n=1000 support points in Fig. 4c, d. Each bar shows the mean of the runtimes on one core of a mobile Intel Core i7 across the 200 experiments for the respective parameter combination; the blue bars are for the νmeasures with uniform masses, the red bars for the measures with masses selected from the corresponding μmeasure. The lines indicate the standard deviations. 123
152 V. Hartmann, D. Schuhmacher We observe that computation times stay more or less the same between parameter choices (with some sampling variation) if the ν-masses are taken from the corresponding μ-measure. In this case mass can typically be assigned (very) locally, and slightly more so if ρhas fewer local fluctuations (higher γand/or s). This seems a plausible explanation for the relatively small computation times. In contrast, if all ν-masses are the same, the computation times are considerably higher and increase substantially with increasing γand somewhat with increasing smoothness. This seems consistent with the hypothesis that the more the optimal transport problem can be solved by assigning mass locally the lower the computation times. For larger scales many of the support points of νcompete strongly for the assignment of mass and a solution can only be found globally. A lower smoothness may alleviate the problem somewhat, because it creates locally more variation in the available mass. In addition to the runtimes, we also recorded how many update steps for the weight vector wwere performed until convergence. We only investigate the update steps for the transport to the original measure ν, not the coarsenings νl,l>0, because the former dominates the runtime, and also has a different dimensionality than the coarsenings. We have computed the Pearson and Spearman correlation coefficients between the numbers of update steps and the runtimes. Both for n=250 and n=1000 support points of ν, these correlation coefficients are larger than 0.99, indicating very high correlation. This strongly suggests that the differences in runtimes are not due to intricacies of the line search procedure or Voronoi cell computations, but rather due to differences in the structures of the simulated problem instances. We would like to note that to the best of our knowledge the present implementation is the first one for computing the optimal transport in the semi-discrete setting for the case p=1, which means that fair performance comparisons with other algorithms are not easily possible. 6 Applications We investigate three concrete problem settings in order to better understand the workings and performance of our algorithm as well as to illustrate various theoretical and practical aspects pointed out in the paper. 6.1 Optimal transport between two normal distributions We consider the two bivariate normal distributions μ=MVN2(a,σ2I2)and ν= MVN2(b,σ2I2), where a=0.8·1,b=2.2·1and σ2=0.1, i.e. they both have the same spherical covariance matrix such that one distribution is just a displacement of the other. For computations we have truncated both measures to the set X=[0,3]2. By discretization (quantization) a measure ˜νis obtained from ν. We then compute the optimal transport partition and the Wasserstein distances between μand ˜νfor both p=1 and p=2. Computations and plots for p=2 are obtained with the package transport (Schuhmacher et al. 2019) for the statistical computing environment R 123
Semi-discrete optimal transport: the unsquared Euclidean distance case 153 Table 1 Theoretical continuous and computed semi-discrete Wasserstein distances, together with the discretization error MVN versus MVN MVN +Leb versus MVN +Leb Theoretical Computed Discr. error Theoretical Computed Discr. error p=1 1.979899 1.965988 0.030962 1.979899 2.164697 0.653370 p=2 1.979899 1.965753 0.034454 Unknown 0.827809 0.220176 (R Core Team 2017). For p=1 we use our implementation presented in the previous section. Note that for the original problem of optimal transport from μto νthe solution is known exactly, so we can use this example to investigate the correct working of our implementation. In fact, for any probability measure μon Rdand its displacement ν=T#μ, where T:R2→R2,x→ x+(b−a)for some vector b−a∈Rd,it is immediately clear that the translation Tinduces an optimal transport plan for (1) and that Wp(μ,ν)=b−afor arbitrary p≥1. This holds because we obtain by Jensen’s inequality (EX−Yp)1/p≥E(X−Y)=b−afor X∼μ,Y∼ν; therefore Wp(μ,ν)≥b−aand Tis clearly a transport map from μto νthat achieves this lower bound. For p=2 Theorem 9.4 in Villani (2009) yields that Tis the unique optimal transport map and the induced plan πTis the unique optimal transport plan. In the case p=1 neither of these objects is unique due to the possibility to rearrange mass transported within the same line segment at no extra cost. Discretization was performed by applying the weighted K-means algorithm based on the discretization of μto a fine grid and an initial configuration of cluster centers drawn independently from distribution νand equipped with the corresponding density values of νas weights. The number of cluster centers was kept to n=300 for better visibility in the plots below. We write ˜ν=n i=1δyifor the discretized measure. The discretization error can be computed numerically by solving another semi-discrete transport problem, see the third column of Table 1below. The first column of Fig. 5depicts the measures μand ˜νand the resulting optimal transport partitions for p=1 and p=2. In the case p=1 the nuclei of the weighted Voronoi tessellation are always contained in their cells, whereas for p=2 this need not be the case. We therefore indicate the relation by a gray arrow pointing from the centroid of the cell to its nucleus whenever the nucleus is outside the cell. The theory for the case p=2, see e.g. Merigot (2011, Sect. 2), identifies the tessellation as a Laguerre tessellation (or power diagram), which consists of convex polygons. The partitions obtained for p=1 and p=2 look very different, but they both capture optimal transports along the direction b−avery closely. For p=2 we clearly see a close approximation of the optimal transport map Tintroduced above. For p=1 we see an approximation of an optimal transport plan πthat collects the mass for any y∈Ysomewhere along the way in the direction b−a. The second column of Table 1gives the Wasserstein distances computed numerically based on these partitions. Both of them are very close to the theoretical value of 123
154 V. Hartmann, D. Schuhmacher ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● μand ˜ ν (a) Measures (b) Measures μ and ˜ ν (c) p=1(d) p=1 (e) p=2(f) p=2 Fig. 5 Left column: semi-discrete transport between a bivariate normal distribution μand a discretized bivariate normal distribution ˜ν. Right column: same with Lebesgue measures added to both distributions (before discretization). Panels aand billustrate the measures. The densities of the continuous measures μ and μare displayed as gray level images, the point masses of the discrete measures νand νare shown as small discs with areas proportional to the masses placed there. Panels cto fshow the optimal transport partitions 123
Semi-discrete optimal transport: the unsquared Euclidean distance case 155 b−a=√2·1.4≈1.979899, and in particular they are well inside the boundaries set by the approximation error. We also investigate the effect of adding a common measure to both μand ν:let α=Lebd|Xand proceed in the same way as above for the measures μ=μ+αand ν=ν+α, calling the discretized measure ˜ν. Note that the discretization error (sixth column of Table 1) is considerably higher, on the one hand due to the fact that the n=300 support points of ˜νhave to be spread far wider, on the other hand because the total mass of each measure is 10 now compared to 1 before. The second column of Fig. 5depicts the measures μand ˜νand the resulting optimal transport partitions for p=1 and p=2. Both partitions look very different from their counterparts when no αis added. However the partition for p=1 clearly approximates a transport plan along the direction of b−aagain. Note that the movement of mass is much more local now, meaning the approximated optimal transport plan is not just obtained by keeping measure αin place and moving the remaining measure μ according to the optimal transport plan πapproximated in Fig. 5c, but a substantial amount of mass available from αis moved as well. Furthermore, Fig. 5dgivesthe impression of a slightly curved movement of mass. We attribute this to a combination of a boundary effect from trimming the Lebesgue measure to Xand numerical error based on the coarse discretization and a small amount of mistransported mass. The computed W1-value for this new example (last column of Table 1) lies in the vicinity of the theoretical value again if one allows for the rather large discretization error. The case p=2 exhibits the distinctive curved behavior that goes with the fluid mechanics interpretation discussed in Sect. 1.2, see also Fig. 1. Various of the other points mentioned in Sect. 1.2 can be observed as well, e.g. the numerically computed Wasserstein distance is much smaller than for p=1, which illustrates the lack of invariance and seems plausible in view of the example in Remark 2in the appendix. 6.2 A practical resource allocation problem We revisit the delivery problem mentioned in the introduction. A fast-food delivery service has 32 branches throughout a city area, depicted by the black dots on the map in Fig. 6. For simplicity of representation we assume that most branches have the same fixed production capacity and a few individual ones (marked by an extra circle around the dot) have twice that capacity. We assume further that the expected orders at peak times have a spatial distribution as indicated by the heatmap (where yellow to white means higher number of orders) and a total volume that matches the total capacity of the branches. The task of the fast-food chain is to partition the map into 32 delivery zones, matching expected orders in each zone with the capacity of the branches, in such a way that the expected cost in form of the travel distance between branch and customer is minimal. We assume here the Euclidean distance, either because of a street layout that comes close to it, see e.g. Boscoe et al. (2012), or because the deliveries are performed by drones. The desired partition, computed by our implementation described in Sect. 4.3, is also displayed in Fig. 6. A number 123
156 V. Hartmann, D. Schuhmacher Fig. 6 The optimal partition of the city area for the delivery example of elongated cells in the western and central parts of the city area suggest that future expansions of the fast-food chain should concentrate on the city center in the north. 6.3 A visual tool for detecting deviations from a density map Very recently, asymptotic theory has been developed that allows, among other things, to test based on the Wasserstein metric Wpwhether a sample in Rdcomes from a given multivariate probability distribution Q. More precisely, assuming independent and identically distributed random vectors X1,...,Xnwith distribution P, limiting distributions have been derived for suitable standardizations of Wp(1 nn i=1δXi,Q) both if P=Qand if P= Q. Based on an observed value Wp(1 nn i=1δxi,Q), where x1,...,xn∈Rd, these distributions allow to assign a level of statistical certainty (p-value) to statements of P=Qand P= Q, respectively. See Sommerfeld and Munk (2018), which uses general p≥1, but requires discrete distributions Pand Q; and del Barrio and Loubes (2018), which is constraint to p=2, but allows for quite general distributions (Pand Qnot both discrete). We propose here the optimal transport partition between an absolutely continuous Qand 1 nn i=1δxias a simple but useful tool for assessing the hypothesis P=Q.We refer to this tool as goodness-of-fit (GOF) partition.Ifd=2, relevant information may be gained from a simple plot of this partition in a similar way as residual plots are used for assessing the fit of linear models. As a general rule of thumb the partition is consistent with the hypothesis P=Qif it consists of many “round” cells that contain their respective P-points roughly in their middle. The size of cells may vary according to local densities and there are bound to be some elongated cells due to sampling error (i.e. the fact that we can only sample from Pand do not know it exactly), but a local accumulation of many elongated cells should give rise to the suspicion that P=Q may be violated in a specific way. Thus GOF partitions provide the data scientist both with a global impression for the plausibility of P=Qand with detailed local 123
Semi-discrete optimal transport: the unsquared Euclidean distance case 163 Mérigot Q (2011) A multiscale approach to optimal transport. Comput Graph. Forum 30(5):1583–1592 Monge G (1781) Mémoire sur la théorie des déblais et des remblais. In: Histoire de l’Académie Royale des Sciences de Paris, avec les Mémoires de Mathématique et de Physique pour la même année, pp 666–704 Nicolas P (2016) Optimal transport for image processing. Habilitation thesis, Signal and Image Processing, Université de Bordeaux. https://hal.archives-ouvertes.fr/tel-01246096v6 Nocedal J (1980) Updating quasi-Newton matrices with limited storage. Math Comput 35(151):773–782 Nocedal J, Wright S (1999) Numerical optimization. Springer Sci 35(67–68):7 Núñez M, Scarsini M (2016) Competing over a finite number of locations. Econ Theory Bull 4(2):125–136 Okazaki N, Nocedal J (2010) libLBFGS (Version 1.10). http://www.chokkan.org/software/liblbfgs/ Peyré G, Cuturi M (2018) Computational optimal transport. now Publishers. arXiv:1803.00567 Pratelli A (2007) On the equality between Monge’s infimum and Kantorovich’s minimum in optimal mass transportation. Ann Inst H Poincaré Probab Stat 43(1):1–13 R Core Team (2017) R: a Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. Version 3.3.0. https://www.R-project.org/ Rippl T, Munk A, Sturm A (2016) Limit laws of the empirical Wasserstein distance: Gaussian distributions. J Multivar Anal 151:90–109 Santambrogio F (2015) Optimal transport for applied mathematicians, Progress in nonlinear differential equations and their applications, vol 87. Birkhäuser/Springer, Cham Schmitz MA, Heitz M, Bonneel N, Ngolè F, Coeurjolly D, Cuturi M, Peyré G, Starck JL (2018) Wasserstein dictionary learning: optimal transport-based unsupervised nonlinear dictionary learning. SIAM J Imaging Sci 11(1):643–678 Schmitzer B (2016) A sparse multiscale algorithm for dense optimal transport. J Math Imaging Vis 56(2):238–259 Schmitzer B (2019) Stabilized sparse scaling algorithms for entropy regularized transport problems. SIAM J Sci Comput 41(3):A1443–A1481 Schmitzer B, Wirth B (2019) A framework for Wasserstein-1-type metrics. J Convex Anal 26(2):353–396 Schrieber J, Schuhmacher D, Gottschlich C (2017) DOTmark: a benchmark for discrete optimal transport. IEEE Access, 5 Schuhmacher D, Bähre B, Gottschlich C, Hartmann V, Heinemann F, Schmitzer B, Schrieber J (2019) Transport: computation of optimal transport plans and Wasserstein distances. R package version 0.11- 1. https://cran.r-project.org/package=transport Sherali HD, Nordai FL (1988) NP-hard, capacitated, balanced p-median problems on a chain graph with a continuum of link demands. Math Oper Res 13(1):32–49 Solomon J, de Goes F, Peyré G, Cuturi M, Butscher A, Nguyen A, Du T, Guibas L (2015) Convolutional Wasserstein distances: efficient optimal transportation on geometric domains. ACM Trans Graph 34(4): 66:1–66:11 Solomon J, Rustamov R, Guibas L, Butscher A (2014) Earth mover’s distances on discrete surfaces. ACM Trans Graph 33(4): 67:1–67:12 Sommerfeld M, Munk A (2018) Inference for empirical Wasserstein distances on finite spaces. J R Stat Soc: Ser B (Statistical Methodology) 80(1):219–238 Villani C (2009) Optimal transport, old and new, Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], vol 338. Springer, Berlin Wolansky G (2015) Semi-discrete approximation of optimal mass transport. Preprint. arXiv:1502.04309v1 Wolfe P (1969) Convergence conditions for ascent methods. SIAM Rev 11:226–235 Wolfe P (1971) Convergence conditions for ascent methods. II. Some corrections. SIAM Rev 13:185–188 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123