Non-asymptotic convergence bounds for Wasserstein approximation using point clouds
Authors: Quentin Mérigot
Filippo Santambrogio
Clément SARRAZIN
Venue: NeurIPS 2021 Poster (2021)
Keywords: Optimal transport
Wasserstein distance
Non-convex optimization
Polyak-Lojasiewicz inequality
Lloyd's algorithm
Point clouds

Abstract: Several issues in machine learning and inverse problems require to generate discrete data, as if sampled from a model probability
distribution. A common way to do so relies on the construction of a uniform probability distribution over a set of $N$ points which
minimizes the Wasserstein distance to the model distribution. This minimization problem, where the unknowns are the positions of the atoms, is non-convex. Yet, in most cases, a suitably adjusted version of Lloyd's algorithm in which Voronoi cells are replaced by Power cells, leads to configurations with small Wasserstein error. This is surprising because, again, of the non-convex nature of the problem, which moreover admits spurious critical points. We provide explicit upper bounds for the convergence speed of this Lloyd-type algorithm, starting from a cloud of points sufficiently far from each other. This already works after one step of the iteration procedure, and similar bounds can be deduced, for the corresponding gradient descent. These bounds naturally lead to a sort of Poliak-Łojasiewicz inequality for the Wasserstein distance cost, with an error term depending on the distances between Dirac masses in the discrete distribution.

-as-asympio convergence bounds for Wasserstein approximation using point clouds
Quentin Mérigot
Filippo Santambrogio
Université Paris-Saclay, CNRS, Laboratoire de mathématiques d' d'Orsay, 91405, Orsay, France Institut Universitaire de France
Univ Lyon,
Université Claude Bernard Lyon 1, CNRS
UMR 5208, Institut Camille Jordan
F-69622 Villeurbanne
Institut Universitaire de France
Clément Sarrazin
Université Paris- Sacia-Sclly, CNRS, Laboratoire de mathématiques d' O'Orsay 91 91405, Orsay, France
Abstract
Several issues in machine learning and inverse problems require to generate discrete data, if sampled from model probability distribution. common way do so relies on the construction of uniform probability distribution over set of N points which minimizes the Wasserstein distance to the model distribution. This minimization problem, where the unknowns are the positions of the atoms, is non-convex. Yet, most cases, suitably adjusted version of Lloyd's algorithm in which Voronoi cells are replaced by Power cells - leads to configurations with small Wasserstein error. This surprising because, again, of the non-convex nature of the problem, as well as the existence of spurious critical points. We provide explicit estimates for the convergence of this Lloyd-type algorithm, starting from cloud of points that are sufficiently far from each other. Our estimates are tight when the algorithm is initialized from an point cloud that is evenly distributed the ambient space. Similar bounds can be deduced for the corresponding gradient descent. These bounds naturally lead to modified Poliak Foliak-Lojasion Lojssiewicz inequality for the Wasserstein distance cost, with error term depending the distances between Dirac masses in the discrete distribution.
Introduction
In recent years, the theory optimal transport has been the source of stimulating ideas in machine learning and in inverse problems. Optimal transport can be used to define distances, called Wasserstein or earth-mover distances, between probability distributions over a metric space. These distances allows one to measure the closeness between a generated distribution and a model distribution, and they have been used with success as data attachment terms in inverse problems. Practically, it has been observed for several different inverse problems that replacing usual loss functions with Wasserstein distances tend to increase the basin of convergence of the methods towards good solution of problem, even to convexify the landscape of the minimized energy [8] 77. This good behaviour is not fully understood, but one may attribute it partly to the fact that the Wasserstein distances encodes the geometry of the underlying space. A notable use of Wasserstein distances in machine learning is the field of generative adversarial networks, where one seeks to design neural network able to produce random examples whose distribution is close prescribed model distribution (2).
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

Wasserstein distance and Wasserstein regression Given two probability distributions p. l on R the Wasserstein distance of exponent p between p and μ a way to measure the total cost of moving mass distribution described by / to j,, knowing that moving unit mass from x to y costs || yll".
Formally, is defined as the value of an optimal transport problem between p and j.:
W Wp(p = (TEl(pp))) / ||= v1l'dT(z,)))
1/r
(1)
where we minimize over the set [p, 14) of transport plans between p and μ,, i.e. probability distributions over IR Rd with marginals p and μ. Standard references on the theory of optimal
transport include books by Villani and by Santambrogio [21 [2 (0], while the computational and
statistical aspects are discussed in survey of Cuturi and Peyré [17].
In this article, we consider regression problems with respect to the Wasserstein metric, which can be put in the following form
min WP WP (T) EEB
where j is the reference distribution, a probability measure on [0, 1J", is the model distribution, a probability measure on R and where T. : [0, 11 - IR a family of maps indexed by a parameter 0 E θ. In the previous formula, we also denoted # the image of the measure l under the map To. also called pushforward of μ under Tg. This image measure is defined by (B) := (TT (B)) for any measurable set B in the codomain of Tg. In this work, we will concentrate on the quadratic Wasserstein distance W2. Several problems related to the design of generative models can be put under the form (2) see for instance 9 2 Alternatively, in some occasions (4) and initially [88), extracting uniform quantization can be used an intermediate step in k- k-means clustering. Let note briefly here that our estimates do not make assumptions on the absolute continuity of the sampled measure and are therefore valid in this fully discrete case. In any case, solving 2 numerically is challenging for several reasons, but in this article we will concentrate on one of them: the non-convexity of the Wasserstein distance under displacement of the measures.
Non-convexity of the Wasserstein distance under displacements. It is well known that the Wasserstein distance is convex for the standard (linear) structure of the space of probability measures, meaning that if vo and 11 are two probability measures and V4 = (1 1)1o + (v,, then the map t e [0,1] + WZ p) is convex. Using a terminology from physics, we may say that the Wasserstein distance is convex for the Eulerian structure of the space of probability measures, e.g. when one interpolates linearly between the densities. However, in the regression problem (2), the perturbations are Lagrangian rather than Fulerian, in the sense that modifications of the parameter 0 leads to displacement of the support of the measure Ta This appears very clearly in particular when μ is the uniform measure over a set X (Z1,.... ENN of N point in [0, 1]d, ie. μ = 6x with
N
8x def
8Z1-
N
In this case To p is the uniform measure over the set Tx(X) (To(x))..... To(xn)), i.e. To (x). In this article, we will therefore be interested by the function
1 FN : Y E (R + W (0,6y).
(3)
2
This function Fx is not convex, and actually exhibits beemi-)ooncavity properties. This has been observed first in [1] (Theorem 7.3.2), and is related to the positive curvature of the Wasserstein space. A precise statement in the context considered here may also be found Proposition 21 in [5]. A practical consequence of the lack of convexity of F'y is that critical points of FN are not necessarily global minimizers. It is actually easy construct examples of families of critical points YN of FN such that Fv (N)) is bounded from below by a positive constant, while lim +00 min FN = 0, so that the ratio between Fv (Yx) and min EN arbitrarily large as N + +o0. This can be done by concentrating the points YN on lower-dimensional subspaces of RRd as in RemarksBandd
When applying gradient descent to the nonconvex optimization problem (2). is in principle possible to end up on local minima corresponding to high energy critical points of the Wasserstein distance,

---

Figure I: From left to right, point cloud Y in the square 52 = [0, 1] x [0, 1/. the associated power cells P,(Y) in the optimal transport to the Lebesgue measure on Ωt the vectors -NVEX(XY)) = BY - y followed during the Lloyd step and the positions of the barycenters Y = (Y).
regardless of the non- con-lincay of the map 0 + tu Our main theorem, or rather its Corollary 6 shows that if the points of Y are at distance least from one another, then
1 Nga-I SNIVE
(x C
In the previous inequality, /V (Y) denotes the Euclidean norm of the vector in RN obtained by putting one after the other the gradients of F'y w.r.t. the positions of the atoms y:- We note that due to the weights 1 1/N in the atomic measure ōy, the components of this vector are in general of the order of 1/N, see Proposition I. This inequality resembles the Malyaa---jjaieiiir inequality, and shows in particular that if the quantization error FN (Y) (2 8y) is large, i.e. larger than (2-d/N. then the point cloud Y is not critical for Fx. From this, we deduce in Theorem 7 that if the points in the initial cloud are not too close to each other at the initialization, then the iterates of fixed step gradient descent converge to points with low energy FN, despite the non-convexity of FN.
Relation to optimal quantization Our main result also has implications in terms of the uniform optimal guantization problem, where one seeks point cloud Y (y1, ................................ yN) in (IRd) such that the uniform measure supported over Y, denoted ōy, is as close as possible to the model distribution with respect to the 2- Wasserstein distance:
min (Y). YERR
4
The uniform optimal quantization problem 44) is a very natural variant of the (standard) optimal quantization problem, where one does not impose that the measure supported on Y is uniform:
min GN(Y), where GN E (IR))N + min Y Y
W? HEAN
P.
(5)
and where AN C R is the probability simplex. This standard optimal quantization problem is cornerstone of sampling theory, and refer the reader to the book of Graf and Luschgy [1] and to the survey by Pagès [16]. The uniform quantization problem (4) less common, but also very natural. It has been used in imaging to produce stipplings of an image 15 3 or for meshing purposes [10] A common difficulty for solving (5) and (4) numerically is that the minimized functionals FN and GN are non-convex and have many critical points with high energy. However, in practice, simple fixed-point or gradient descent stategies behave well when the initial point cloud is not chosen adversely. Our second contribution a quantitative explanation for this good behaviour in the case of the uniform optimal quantization problem.
Lloyd's algorithm [13 is fixed point algorithm for solving approximately the standard optimal quantization problem (5). Starting from a point cloud y = (vf,-. v E (RRd)N with distinct points, one defines the next iterate y* in two steps. First, one computes the Voronoi diagram of Y, tessellation of the space into convex polyhedra (V.(Y* )))ssN)),,, where
Vi(Y) (x E Ω Vj E (1, | vill s ||x -
In the second step, moves every point yy towards the barycenter, with respect p, of the corresponding cell V,(K)). This algorithm can also be interpreted as a fixed point algorithm for
3

---

solving the first-order optimality condition for (5), i.c. (Y) = 0. One can show that the energy (Gx(Y")) decreases k. The convergence of y towards a critical point of Fy as - +o0 has been studied in [6]. but the energy of this limit critical point not guaranteed be small.
In the case the uniform quantization problem one can try to minimize the energy PN by gradient descent, defining the iterates
yk+1 y TNVFF(YY".
where is the time step. The factor N in front FN is set as a compensation for the fact that we have, in general, VFN O1/NN). When = 1, one recovers version of Lloyd's ss algorithm for the uniform quantization problem, involving barycenters By (Y) of Power cells, rather than Voronoi cells, associated to Y. More precisely, Proposition proves that FF) = (Y (Y))))/ so that y = Bx(Y when = 1. Quite surprisingly. we prove in Corollary 4 that if the points in the initial cloud Y are not too close each other, then the uniform measure over the point cloud Y = yu NV Fx (Y obtained after only one step of Lloyd's algorithm is close to This is illustrated in Figure M We prove in particular the following statement.
Theorem (Particular case of Corollary 4). Let p be a probability density over a compact convex set Ω S RR, let Y (yf.... v S2d and assume that the points lie at some positive distance from one another: for some constant
Vi j. v11 cN-1/4
corresponding for instance to point cloud sampled a regular grid. Then, the point cloud -NV obtained after one step of Lloyd's 's algorithm satisfies
W p) S Ce
where Ce.d,1 is a constant depending on c, d and diam([).
Outline Section 2 we start by a short review of background material on optimal transport and optimal uniform quantization. We then establish our main result (Theorem B the approximation of a measure p by barycenters of Power cells. This theorem yields error estimates for one step of
Lloyd's algorithm in deterministic and probabilistic settings (Corollaries and 5). In Section 3
establish a Polyak-Lon Lojyak--jssiewiic-typee inequality (Corollary 6) for the function FN Y + 1 W3 (p. dy
introduced in (3), and study the convergence of gradient descent algorithm for Fv (Theorem 7 Finally, in Section 4 we report numerical results on optimal uniform quantization dimension d= d=2 2.
2 Lloyd's algorithm for optimal uniform quantization
Optimal transport and Kantorovich duality this section we briefly review Kantorovich duality and its relation to semidiscrete optimal transport. The cost is fixed to c(x, y) =| | assume that p is a probability density over a compact convex domain Ω. In this setting, Brenier's theorem implies that given any probability measure μ supported on Ω the optimal transport plan between p and j., i.e. the minimizer æ in the definition of the Wasserstein distance (T) with p = 2, is induced by transport map Tμ Ω Ω, meaning æ (Tμ, Id) p. One can derive alternative expression for the Wasserstein distance using Kantorovich duality, which leads to a more precise
y11² and
description of the optimal transport map (00 Theorem 1.39]:
/ odu
/ d'dp,
(2. H)
max o:Y +R gd
(
where (x) min, c(z, B1) o. When 8y is the uniform probability measure over point cloud = (y1,..., yN) containing N distinct points, we set o. = o(y.) and we define the ith Power cell associated to the couple (Y, v))
Pow,(Y,6) = (x E Rd
| - 11 di < || wzll² -z)-
Then, the Kantorovich dual 8 of the optimal transport problem between and 8y turns into a finite dimensoimensional concave maximization problem
W2(
max DER
from:(Yss)) (1 wil12 - 0i) d)(())

---

By Corollary 1.2 in [2], vector & IN is optimal for this maximization problem if and only if the potential is such that each Power cell contains the same amount of mass, i.c. if
1
Vi E Ny, p(Pow, N
(10)
From now on, we denote P.(Y) = PY, Φ) n where á E R satisfies (10) The optimal transport map Ty between P and dy sends every Power cell P.(Y) the point yt. i.c. is defined
p-almost everywhere by Ty P.(Y)
We refer again to the introduction of [2 for more details.
Optimal uniform quantization In this article, study the behaviour of the squared Wasserstein distance between the (fixed) probability density p and a uniform finitely supported measure 5y where
Y (11..... ......... .UN) is cloud of N points, in terms of variations of Y. As in equation (3). we denote Fx = (p..). Proposition 21 in [13] gives an expression for the gradient of F, and proves its semiconcavity. We recall that F is called a-seemiconcave, with o 2 0, if the function F - 2 concave. We denote y the generalized diagonal
DN = (Y E (R4) 3i #j yr = =y}}.
Proposition 1 (Gradient of Fx). The function FN is -eemicooncvee on (R4) and is of class C² on (IR N N In addition, one has
1 AY E (RA) A DN, VFN(Y) = N -
where Bx = =(b((Y)....b>(>))) (11)
and where bY) is the barycenter of the ith power cell, i.e. by(Y)) = N Sp,(r zdp(x).
It is not difficult to prove that FN admits at least one minimizer, and that this minimizer Y satisfies the first-order optimality condition Y Bx point cloud that satisfies this condition is called critical.
Remark 1 (Upper bound the minimum of FN). We note from 15 Proposition 12] that when p is supported on a compact subset of IR, then
-
if 2
min Fny
min YE(RA
v W3(p,or)
log
if [ 2 if d 1.
(12)
Remark (High energy critical points). On the other hand, since FN is not convex, this first-order condition is sufficient have minimizer Fy. For instance, if = on the unit square f? = [o, 1]² can check that the point cloud
((2))) YN (( 2N 2
3
2N - )
2N
2.N
a critical point of Fy but not minimizer of Fx. In fact, this critical point becomes arbitrarily bad as N + +oo in the sense that
Fx(YN) +00.
lim +00 min FN
On the other hand, note that the point cloud YN is highly concentrated, in the sense that the distance between consecutive points in YN is 2N. whereas in an evenly distributed point cloud, one would expect the minimum distance between points to be of order N-
Gradient descent and Lloyd's algorithm One can find critical point of Fy by following the discrete gradient flow of FN, defined in 77. starting from initial position y e (IRd) Dy. Thanks to the expression of FN given in Proposition [ the discrete gradient flow may be written as
k+1 yk + (T((BN(( y*),
(13)
5

---

where TN is fixed time step. For TN = 1, one recovers a variant of Lloyd's algorithm, where one moves every point to the barycenter of its Power cell P; (Y* at each iteration:
y = BN (Y**))
(14)
We can state the following result about Lloyd's algorithm for the uniform quantization problem, whose proof is postponed the appendix.
Proposition 2. Lel N be fixed integer and (Y*))20 be the iterates of (1). with yo & D Then, the energy k + (Y*) is decreasing, and limg->+oe [ = O. Moreover, the sequence (z2 belongs to a compact subset of (R4) \ Dy and every limit point of a converging subsequence of it is critical point for Fn.
Experiments suggest that following the discrete gradient flow of Fy does not bring us high energy critical points of FN. such as those described in Remark 2 unless we started from an adversely chosen point cloud. The following theorem and its corollaries, the main results of this article, backs up this experimental evidence. It shows that if the point cloud Y is not too concentrated, then the uniform measure over the barycenters of the power cells, (Y), ,, is a good quantization of probability density i.e. it bounds the quantization error after one step of Lloyd's algorithm (14
We will use the following notation for 0:
(VY (i (1...., N) Vj # 1, yy 2 e)
DN,= = {Y E (d 3i + j. 11:/ M2 E Note that D, is an E-neighborhood around the generalized diagonal DN-
Theorem 3 (Quantization by barycenters). Let Ω C Rd be compact convex set, p probability density on and consider a point cloud Y (U1,..., .NN 22 \ TTx. Then, for all <es1,
W3 (p. ()) & Ca. (
E (s1-d
+1 Card(4((X)))
(15)
N
N
where Ca.
224-1 Wa
(!)) + 1)d+1 and where wa-1 is the volume of the unit ball in IR-1
The proof relies on arguments from convex geometry. In particular, we denote A a B the Minkowski sum of sets: A + = {a + b (a, b) A x B}.
Proof. Let o' R be the solution to the dual Kantorovich problem (00) between p and ōy. We let o' = =o and we denote p = Pow,(X (Y, o') n L the ith Power cell intersected with the slightly enlarged convex set ! l2 + (00 1). This way, P. 2 P.YY) whereas P is in fact the intersection of the i-th Voronoi cell defined in (6) with Ω
We will now prove an upper bound on the sum of the diameters of the cells P.(Y) whose index lies in (Y)). First, we notice the following inclusion, which holds for any e [0, 1):
(1 (1- P) e LPP Ct,
(16)
Indeed, let x0 E pd and x E P!,, so that for all j E (1..... N) and e (0.1}.
||** 22 d; < ||* y
Expanding the squares and substracting || 1 on both sides these inequalities become linear in Φs. 2 and x² implying directly that x' (1 t)z² 0 + tx S P. as desired.
For any index E Ig. the point y: is at distance at least e from other points, implying that ( f) is contained in the Voronoi cell Vi(Y) with Ω Using that po = Vi(Y) n that Ω = Ω @ B 1) and that y/ E S2 we deduce that 1.00 contains the same ball. On the other hand, P! contains segment St of length diam(P!) and inclusion (6) with t gives
(B(ys, e/2) @ S.) C / 2
The Minkowski sum the left-hand side contains in particular the product of (d - -iddimensional ball of radius e/2 with an orthogonal segment with length diam(P!) 2 dY))( Thus.
1
d-1
diam diam((.()))) (((X))) < |
Wa Wa-1 2
2d

---

where I-| denotes the Lehesgue measure. Using that the Power cells P form a tesselation of the domain S2' we therefore obtain
22d-1
22d-1
diam (P(()) -
(di + 1)del-d
(17)
wa-1
Wd-1
We now estimate the transport cost between on and the density p. where B By (Y). The transport cost due to the points whose indices do not belong (Y) can be bounded in a crude way by
Σ Σ
|| b, (x S Card (Y) diam
ig1,0 P.0))
N
Note that we used p(P:(Y)) (Y) can be bounded using
On the other hand, the transport cost associated with indices in 17 and (iaam(:(()) I diam((?);
Σ
x 18-6,(3 b,)r)
Σ
diam(P.(Y))²
iEL(Y) PP)
iETL(Y)
(1)
diam (P((Y)
N
xET
2n 1
(diam(!)) + 1) E l-d
Wa-1
N
In conclusion, obtain the desired estimate:
22d-1
(diam(St) + 1) s--d W
N + diam(S))² (1- Card Ie
W2 (p.
N
This theorem could be extended mutatis mutandis to the case where & is a general probability measure (i.e. not density). However, this would imply some technical complications in the definition of the barycenters b; by introducing a disintegration of p with respect to the transport plan æ.
Consequence for the uniform Lloyd's algorithm (1) In the next corollary, we assume that any pair of distinet points in Yv E (R4y bounded from below by nN CN-B, implying that Iex (Yn) N. This corresponds to the value could expect for point set uniformly sampled from set with Minkowski dimension 3. When 3 > d - 1, the corollary asserts that one step of Lloyd's algorithm is enough to approximate p. in the sense that the uniform measure d (Yw) over the barycenters converges towards p as N +oc.
Corollary 4 (Quantization by barycenters, asymptotic case). Assume EN 2 C. 1-8 C-N-1/8 with C, O. Then, with 1 d-1
Ca.
VY E (R DN W2lp.dna(n) S -
(18)
and in particular, B >
1,
max +++ Y ((
lim
Wa(p, =.
(19)
Remark 3 (Optimality of the exponent when 3= d). There is no reason to believe that the exponent in the upper bound (18) is optimal in general. However, it seems be optimal in "worst-case sense" when 3 d. More precisely, we show the following result in Appendix E for any 8 e (0, 1). and for every N nd (n E N) there exists a sequence of separable probability densities aN over X [-1, 1] (pn is truncated Gaussian distributions, whose variance converges to zero slowly as N + +oo) such that if YN is a uniform grid of size n X X N' in X, then
W2 PN) 2 CN 12-4)
where independent of N. On the other hand, in this setting every point in YN is distance at least 11/ from any other point in Yx. Applying Corollary 4 with d, i.c. a J. we get
W (6 S C'N-à.
Comparing this upper bound on W (5 with the above lower bound, one sees that is not possible to improve the exponent.
7

---

Remark 4 (Optimality of (19). The assumption 8 > d - 1 for 19 is tight: if p is the Lebesgue measure on [0. 11 is possible to construct point cloud YN with N points on the (d 1) cuue (2) x [0, 1]0-1 such that distinct point in Yv are at distance at least EN 2 C..N CCN-I/(d-1) Then, the barycenters By (Yx) are also contained in the cube, so that (2 (Bw 12
The next corollary is probabilistic analogue of Corollary | assuming that the initial point cloud Y is drawn from probability density on f2. Note that a can be distinct from The proof of this corollary relies on McDiarmid's inequality to quantity the proportion of s-isolated points in a point cloud that is drawn randomly and independently from a. The proof of this result is in Appendix B
Corollary 5 (Quantization by barycenters, probabilistic case). Let L (SI) and let X1, ..., X be i..dd random variables with distribution E L (R*). Then, there exists constant C > 0 depending only on |o and d, such that for N large enough,
P
W2
1 N
Σ 0ogi0)
N
)
2 21 21-e-CNa CCN
24-3
3 Gradient flow and a
inequality
Theorem 3 can be interpreted as modified
(PL for short) inequality for the function FN. The usual PL inequality for a differentiable function F IR - E is of the form
VY E RR, F(Y) min F I C))VE(Y))??
where C a positive constant. This inequality has been originally used by Polyak to prove conver- gence of gradient descent towards the global minimum of F. Note in particular that such an inequality implies that any critical point of F is a global minimum of F. By Remark Fx has critical points that are not minimizers, so that we cannot expect the standard PE. inequality to hold. What we get a similar inequality relating FN (Y) and [/ (2)12 but with a term involving the minimimum distance between the points in place of min FN.
Corollary (Polyak-L.jjaa foolyak--ojasiewizz-yyy inequality). Let Y E (IRA)N \ NN.. Then,
1 FY Ca.a N
(
d-1
S (NIVFN((Y))
(20)
A proof of Corollary 6 can be found in Appendix C
We note that when & = N-1/d, the second term of the left-hand side of (20) has order N-1/2, On the other hand, as recalled in Remark | min FN S 2/d when d > 2. Thus, (20 is strictly weaker than the PL inequality, which would
involve the minimum of FN.
Convergence of discrete gradient flow The modified Lojyyak--ojsiiiiizzzz inequality (0) sug- gests that the discrete gradient flow [3 will bring us close point cloud with low Wasserstein distance to p. provided that we can guarantee that the the points clouds Y remain far from the generalized diagonal Dx during the iterations. We prove Lemma 3 in Appendix D that
ykk+1 = yk TN V Px (Y*) and TN E (0, 1), then
Vi j.
| 1
21 (1-TN))l 1/
(21)
We note that this inequality ensures that y never touches the generalized diagonal Dx, that the gradient FN (Y*) is well-defined at each step. Combining this inequality with Theorem 3 one can actually prove that if the points in the initial cloud Y are not too close to each other, then a few steps of gradient discrete gradient descent leads to discrete measure Y that is close to the target p. We stress here that the goal of this paragraph is to showcase these kinds of several-steps estimates, in the simple case of our gradient flow for Fx. However, in this case, it provides a worse convergence rate to (w.r.t. N) than that given by Corollary and the decrease of the energy along the iterations from Proposition
2 Precisely, we arrive at the following theorem, proved in Appendix D
Theorem 7. Let 0 < a < र J. eN 2 N - and Y N \ [BE Let &Y ) be the iterates
of 133) starting from Y with timestep 0 < TN 1. We assume that Iim IN and we set
1
KN =
In IF (X))
drN
8

---

Figure 2: Optimal quantization of density P corresponding gray-scale image (Wikimedia Commons, CC BY-SA 3.0). (Left) We display the point clouds obtained after one step of Lloyd's algorithm, starting from regular grid of size N E (3750, 7350, 15000, 43350). (Right) Quantization error W (p. dBw) as function of N the number of points, showing that W2 (p. dBw) = N-0.00
Then,
W (p. 8ybe
=N ( +
(22)
Remark 5. Note that the exponential behavior implied by 21 and Lemma 3 is coherent with the estimates that are known in the absolutely continuous setting for the continuous gradient flow. When transitioning from discrete measures probability densities, lower bounds on the distance between points become upper bounds on the density. The gradient flow jrę = 1/2 V W pe) has an explicit a constant-speed geodesic in the Wasserstein space with do = po and 01 = p. In this case, a simple adaptation of the estimates in Theorem 2 in [9] shows the bound |4 - I edd [/olll== Still in this absolutely continous setting, is possible to remove the exponential growth if the target density is also bounded, as a consequence of displacement convexity [4 Theorem 2.2]. There seems be no discrete counterpart this argument, explaining in part the discrepancy between the exponent of N in 22 with the one obtained in Corollary 4
solution [4 = 01 01-e-- 1, where
4 Numerical results
In this section, we report some experimental results in dimension d 2.
Gray-scale image As we mentioned in the introduction, uniform optimal quantization allows to sparsely represent (gray scale) image via points, clustered more closely in areas where the image is darker 3). On figure we ploted the point clouds obtained after single Lloyd step toward the density representing the image on the left (Puffin), starting from regular grids. The rate of convergence observed on the sanht side chart, namely N 1100, is coherent with the theoretical estimate log (NN of Remark
Gaussian density with small variance We now consider toy model where we approximate gaussian density truncated the unit square [0, 1)2, p(a, v) = e -8((=-(()))))) +v= 4)*) where Z a normalization constant. On the left column of Fig. the initial point clouds Y are randomly distributed in [0, 1ja. The three point clouds represented above are obtained after one step of Lloyd's algorithm (1).. The blue curve displays in a log-log scale the mean values of (B (Yn)) over a hundred random point clouds, for N E {400, 961, 1600, 2500). In this case, we observe a decrease rate N 0.95 with respect to the number of points, similar to the case of the gray scale picture.
However, an interesting phenomena occurs when the initial point cloud Y is aligned on a axis- aligned grid. The pictures the right column of Fig. 3 where computed starting from such grid with N E {400, 961, 1600, 2500} points. As in the randomly initialized case, we represented the values of Fx xBN (Yx)) )) in log- scale. The corresponding discrete probability measure (xx) seems to converge to p as N + oc, but with much worse rate for these "low" values of N: Fx (Yw))
- In this specific setting, with a separable density and an axis-aligned grid Yo the power cells are rectangles and a single Lloyd step brings to critical point of Fx. Thanks to this remark, it is possible to estimate the approximation error from the one-dimensional case. In fact, Appendix E shows that for any 5 E (0, 1), there exists variances ON x8) such that

---

Figure Optimal quantization of Gaussian truncated to the unit square. On the left, the initial point cloud YN is drawn randomly and uniformly from [0, 1 while on the right YN is on a regular grid. The top row displays the point clouds obtained after one step of Lloyd's algorithm. The bottom row displays the quantization error after one step of Lloyd's algorithm Fy (Bx (Yx)) as fuction of the number of points. We get Fn (By (Yn)) N-0.95 when Yy is a random uniform point cloud in [0, 1] and EN (N((N(YN) ~ N 0.8 when YN is regular grid.
the approximation error W (PaN. 5B is of order ~ On the other hand, for fixed a, the approximation enror is of order N to be compared with the bound log(N)) for general measures.
5
Discussion
We have studied the problem of minimizing the Wasserstein distance between a fixed probability measure p and uniform measure over N points öy, parametrized by the position of the points
== (31..... ...,, (N). The main difficulty is the nonconvexity of the Wasserstein distance FN Y E (p, 5y), which we tackled by means of a modified Polyak-toojccewiic inequality 20) One limitation of work is that the exponents in estimates depends the ambient dimension: would be more interesting (but quite nontrivial) to establish similar estimates depending on the dimension of the support of p. Another direction of research would be derive consequences for the algorithmic resolution of Wasserstein regression problems ming ( To starting with the case where 0 + To is linear.
[Rd
Acknowledgments and Disclosure of Funding
This work was supported by grant from the French ANR (MAGA, ANR-16-CE40000144)
10

---

References
[1] Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient flows: in metric spaces and in the space of probability measures. Springer Science & Business Media, 2008.
[2] Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In International conference on machine learning, pages 214-223. PMLR, 2017.
[3] Michael Balzer, Thomas Schlömer, and Oliver Deussen. Capacity-constraine point distribu- tions: a variant of lloyd's 's method. ACM Transactions on Graphics (TOG), 28(3):1 2009.
[4] Marco Cuturi and Arnaud Doucet. Fast computation of wasserstein barycenters. In Eric P. Xing and Tony Jebara, editors, Proceedings of the 31 International Conference Machine Learning, volume 32 of Proceedings of Machine Learning Research, pages 085-693, Bejing, China, 22-24 Jun 2014. PMLR.
[5] Fernando De Goes, Katherine Breeden, Victor Ostromoukhov, and Mathieu Desbrun. Blue noise through optimal transport. ACM Transactions on Graphics (TOG), (11: 2012. 2012.
[6] Qiang Du. Maria Emelianenko, and Lili Ju. Convergence of the lloyd algorithm for computing centroidal voronoi tessellations. SIAM journal on numerical analysis, 44(1):102-119; 2006.
[7] Bjorn Engquist, Brittany D Froese, and Yunan Yang. Optimal transport for seismic full waveform inversion. Communicationn in Mathematical Sciences, 14(8):2309-2300 23300,,,,,,, 2016.
[8] Jean Feydy, Benjamin Charlier, François-Xavis Vialard, and Gabriel Peyré. Optimal transport for diffeomorphic registration. In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages
Springer, 2017.
[9] Aude Genevay, Gabriel Peyré, and Marco Cuturi. Learning generative models with Sinkhorn divergences. In International Conference on Artificial Intelligence and Statistics, pages 1608- 1617. PMLR, 2018.
[10] Fernando de Goes, Pooran Memari, Patrick Mullen, and Mathieu Desbrun. Weighted tri- angulations for geometry processing. ACM Transactions on Graphics (TOG), 33(3);1-13, 2014.
[11] Siegfried Graf and Harald Luschgy. Foundations of quantizaation for probability distributions. Springer, 2007.
[12] Jun Kitagawa, Quentin Mérigot, and Boris Thibert. Convergence of newton algorithm for semi- discrete optimal transport. Journal of the European Mathematical Society, 21(9):2603-28) 2651, 2019.
[13] Stuart Lloyd. Least squares quantization in pem. IEEE transactions on information theory, 1892)11 1982.
[14] Robert J McCann. A convexity principle for interacting gases. Advances in mathematics, 128(1):153-179, 1997.
[15] Quentin Mérigot and Jean-Marie Mirebeau. Minimal geodesics along volume-preserving maps, through semidiscrete optimal transport. SIAM Journal on Numerical Analysis, 54(6):3465-33) 3492, 2016.
[16] Gilles Pagès. Introduction to vector quantization and its applications for numerics. ESAIM: proceedings and surveys, 48:29-79, 2015.
[17] Gabriel Peyré, Marco Cuturi, al. Computational optimal transport: With applications to data science. Foundations and Trends@ in Machine Learning. 115-6))355-660 607, 2019.
[18) Sebastian Reich. A nonparametric ensemble transform method for bayesian inference. SIAM Journal on Scientific Computing. 35(4):22013-2222,, 2013.
[19] Filippo Santambrogio. Absolute continuity and summability of transport densities: simpler proofs and new estimates. Calculus of variations and partial differential equations, ( 354, 2009.
[20] Filippo Santambrogio. Optimal transport for applied mathematicians, volume 55. Springer, 2015.
[21] Cédric Villani. Topics in optimal transportation. Number 58. American Mathematical Soc., 2003.

---

[22] Cédric Villani. Optimal transport: old and new, volume 338. Springer Science & Business Media, 2008.
12


Reviews:
[{"review_id": "jsWNyWWBa2", "reviewer": "Reviewer_t3bd", "summary": "This paper studies approximation of continuous measures in Wasserstein distance using point clouds, which define empirical distributions over their support. Optimizing over the locations of the points in these clouds results in a nonconvex optimization problem. Despite this, the authors develop a Polyak-Lojasiewicz-like inequality that states that, from a well-separated initialization, one can get to a point that is close to the underlying measure (on the order of $N^{-1/d}$). An iterative version of this method is shown to converge with a slower order.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Upsides:\n\n+ This paper tackles an important problem of approximation of measures in Wasserstein distance by discrete distributions and gives the first theoretical results for this problem.\n+ The PL-like inequality is new and interesting from a nonconvex optimization perspective. In particular, it gives an approximation error for a nonconvex setting that decreases with the sample size. This type of analysis may have broader applications.\n+ The main results and their proofs are clear.\n\nDownsides:\n\n- The experiments are not extensive, and rather serve as a proof of concept. Furthermore, they do not illustrate the curse of dimensionality encountered in the theorems.\n- How does one compute the barycenter of the Voronoi cells in practice? What is the general computational complexity and practicality of implementing this method?\n- Theorem 7 and Remark 5 are not sufficiently clear. In particular, the discrepancy in the rates of Theorem 7 and Corollary 4 is not clearly expounded on. In general, it would be helpful to collect the results for the various methods/theorems together to make them easier to interpret and compare.\n- In some sense, for the one-step algorithm, the method seems to be learning the barycentric projection optimal transport map, provided that the initial $Y^0$ is drawn from some density. This would then explain the types of rates that are seen in this setting, due to the rates of estimation for the optimal transport distance. It would be nice if the authors could comment on the connection here.\n- The guarantees do not give hints as to what types of measures are easier to approximate than others. Furthermore, the dimension dependence of the constants is not discussed.\n\n---\n\nI still find the result interesting and think that it is close to being ready for publication. However, due to some issues raised by other reviewers, as well as the limited practical experiments, I will leave my score at a 6.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations appear to be discussed in the conclusion. Societal impact/applications of the methodology are not adequately discussed. When does one which to quantize a distribution with an empirical distribution, and why does it make sense to do this with Wasserstein distance?", "review_text": "", "extra_scores": {"ethical_concerns": "N/A", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "dZUthAWixD5", "reviewer": "Reviewer_Xcrv", "summary": "This paper studies the optimal quantization problem, where one is tasked with fitting a discrete uniform measure supported on a fixed number of points to a continuous measure, measured by 2-Wasserstein loss. Their primary theoretical finding is that one step of gradient descent on this objective, makes progress proportional to the minimum distance between the starting uniform measure's support.", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The main finding in this work is that a well-separated discrete measure leads to progress in gradient descent for the optimal quantization problem. This result is achieved by a short argument which proceeds by showing that a well-separted discrete measure must induce well-conditioned power cells, and so their barycenters must have relatively little transport cost. The result is further developed by combination with a previous observation that their objective is semi-concave, which allows some understanding of multiple steps of gradient descent.\n\nWhile the primary result is interesting and novel, I feel that there is insufficient context and empirical or theoretical development to merit acceptance. In particular, I wonder whether this result is of questionable significance since in practice, to perform the gradient descent step one must estimate integrals of the density over the relevant power cells. If this is done via Monte-Carlo (which would be standard in a high-dimensional context), then the natural competitor would be the empirical measure on a sample of size $N$, which would achieve rates of the order $N^{-1/d}$ [1]. How does this compare to the estimator proposed here in terms of relative computational/statistical efficiency? If theoretical results are difficult, can an empirical comparison be made?\n\nAnother issue is that I feel the context for the obtained rates is not clearly explained. This is somewhat addressed in Remark 3, but I feel it should be more clearly explained that generically we expect for $N$ points and a \"nice\" density supported a bounded subset of $\\mathbb{R}^d$, we cannot do better than $N^{-1/d}$ (this being a strength of the obtained result that is not properly expounded). Finally, I feel the experimental results section could be improved by a much more thorough substantiation of the claim that Lloyd's algorithm typically finds good local minima. Or, if there is already such a work in the literature, it should be cited (apologies if this has been cited and I missed it).\n\n[1] Fournier and Guillin, \"On the rate of convergence in Wasserstein distance of the empirical measure\"", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors addressed limitations adequately other than this issue regarding the situation of applying Lloyd's algorithm in practice.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "JbstFZxjMPY", "reviewer": "Reviewer_fudr", "summary": "The authors study the semi-discrete optimal transport problem between the source measure, which is the distribution of a set of points after one iteration of the Lloyd's algorithm and the target continuous measure. In particular, they prove that the 2-Wasserstein distance between these two probability distribution can be approximation with a given accuracy, depending on $N$ the number of support points of the discrete measure and $d$ the dimension of the underlying Euclidean space, providing that the initial points in the initial point cloud are sufficiently spaced. A more general result on the convergence of Lloyd's algorithm towards a critical point of the Wasserstein distance with respect to the point cloud is also provided. These issues are related to uniform optimal quantization problems.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality : To my knowledge, this paper is the first to provide accuracy bounds for the Wasserstein distance between a continuous measure and a point cloud after the first step of Lloyd's algorithm, reducing the gap between experimental and theoretical results, the latter being mostly asymptotic (in particular in the optimal quantization community). In particular, this paper proves the advantage of Lloyd's algorithm, even if the functional $F_N$ to optimise wrt to the point cloud $Y$ is non-convex.\n\nQuality : This paper contains strong theoretically results, for which the assumptions and consequences are clearly stated and explained. Moreover, some unfamiliar behaviours, such as in Figure 2 (right) are still presented and commented.\n\nClarity : This paper well written and the problems are clearly stated. Moreover, natural results following the first part of the paper, such as Corollary 5 (probabilistic case of quantization) are also proven. However, in my opinion, the introduction of regression problems is not completely justified as it can be more confusing than anything else (it is neither used or leveraged in the applications for example).\n\nSignificance : Lloyd's algorithm is largely used in several communities, from optimal transport to optimal quantization problems, even if it only solves provide critical points of a non-convex problem. Therefore, such a study on the behaviour of Lloyd's algorithm, both in finite time and asymptotically, is of great interest.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations : As the authors discuss in the conclusion, the bound obtained for Lloyd's algorithm does not reach the theoretical rate of convergence, even if the initialisation points are sufficiently far apart.\n\nQuestions : - Considering the bounds on $W_2^2(\\rho, \\delta_{B_n(Y)})$, would the gain after $k$ iterations of Lloyd's algorithm be significantly better than after a single step?\n-  In practice, Lloyd's algorithm can produce solutions with some boundary issues on the compact domain $\\Omega$. From your results, I would say that this behaviour is avoided when the points in $Y$ are sufficiently far away. Could you comment on this?\n\nMinor comments : \n- line 60 : the pushforward operator should be above, i.e. $T_{\\theta}$ # $\\sigma$.\n- In Proposition 1, $Y$ is defined in both $DN$ and ${D}{N,\\varepsilon}$. Moreover, in the definition of the barycenter of the power cell, the variable $x$ is missing (it is written $b_i(Y) = N\\int_{P_i(Y)}d\\rho(x)$).\n- In Theorem 3, the point cloud $Y$ should belong to $\\mathbb{D}_{N,\\varepsilon}$.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "AYYuzwA3qf_", "reviewer": "Reviewer_AWCx", "summary": "This paper considers the problem of constructing a uniform probability distribution over a set of N points which minimizes the Wasserstein distance to a known reference distribution. This minimization problem, where the unknowns are the positions of the atoms, is non-convex. However, it can be shown that a suitably adjusted version of Lloyd's algorithm in which Voronoi cells are replaced by Power cells, leads to configurations with small Wasserstein error. Explicit upper bounds for the convergence speed of this Lloyd-type algorithm are presented in Theorem 7, starting from a cloud of points sufficiently far from each other. These bounds also lead to a sort of Poliak-Łojasiewicz inequality for the Wasserstein distance cost, with an error term depending on the distances between Dirac masses in the discrete distribution.", "questions": "", "limitations": "", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper considers an interesting problem, i.e., minimizing the Wasserstein distance between a known continuous distribution and a discrete distribution defined over a set of N points. Such a problem cannot be solved in KL divergence, as the KL divergence between any continuous distribution and discrete distribution is always infinite. \n\nHowever, I am a bit confused by the contribution of the paper; it seems that Lloyd’s algorithm based on power cell has been studied in [11] and [14] under similar setting, see, equation (10) and Proposition 1 in the current paper. The authors should make it clear whether the adjusted Lloyd’s algorithm is a main contribution of the paper or not. In addition, more details about how to efficiently compute the barycenter for the power cell should be provided.\n\nIn general, there are some merits in the paper, but it is not well-written in the current form, and a lot of details need to be clarified. \n\nMinor comments:\n1.\tIn Line 118, the notation $(T_\\mu, Id)$# is used without definition\n2.\tIn line 62, should be “In the previous inequality”.\n3.\tIn line 133, the definition of barycenter is not correct, the current definition gives $b_i(Y)=1$, should be $b_i(Y)=N\\int_{P_i(Y)}x d\\rho(x)$.\n4.\tThe order of Figure 2 and 3 should be swapped. \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "This work is mainly theoretical, its practical impact will depend on the specific applications.\n\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers recognize the theoretical interest of the authors' contribution, bounding in particular the error done after one step of gradient descent to minimize Wasserstein distance between a given density and a uniform discrete N-point distribution. They however raise valid concerns, particularly on the complexity needed to compute such gradients, on how it would rely on obtention of samples from the density being approximated. These are important aspects that are not treated in the paper. The reviewers also point to the lack of a discussion contrasting the guarantee they obtain to the one that results from taking as an approximating point cloud an N-sample from the target density. 

It therefore appears that the paper contains a very promising set of results, but is at this stage too preliminary and not ready yet for publication at NeurIPS.
Decision: Accept (Poster)