Semi-Discrete Optimal Transport:  Nearly Minimax Estimation With Stochastic Gradient Descent and Adaptive Entropic Regularization

Authors: ferdinand  genans
antoine  godichon-baggioni
françois-xavier  vialard
olivier  wintenberger
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 10:54:47

Abstract: Optimal Transport (OT) based distances are powerful tools for machine learning to compare probability measures and manipulate them using OT maps. In this field, a setting of interest is semi-discrete OT, where the source measure $\mu$ is continuous, while the target $\nu$ is discrete. Recent works have shown that the minimax rate for the OT map is $\mathcal{O}(t^{-1/2})$ when using $t$ i.i.d. subsamples from each measure (two-sample setting). An open question is whether a better convergence rate can be achieved when the full information of the discrete measure $\nu$ is known (one-sample setting). In this work, we 
    answer positively to this question by (i) proving an $\mathcal{O}(t^{-1})$ lower bound rate for the OT map, using the similarity between Laguerre cells estimation and density support estimation, and (ii) proposing a Stochastic Gradient Descent (SGD) algorithm with adaptive entropic regularization and averaging acceleration. To nearly achieve the desired fast rate, characteristic of non-regular parametric problems, we design an entropic regularization scheme decreasing with the number of samples. 
    Another key step in our algorithm consists of using a projection step that permits to leverage the local strong convexity of the regularized OT problem. Our convergence analysis integrates online convex optimization and stochastic gradient techniques, complemented by the specificities of the OT semi-dual. Moreover, while being as computationally and memory efficient as vanilla SGD, our algorithm achieves the unusual fast rates of our theory in numerical experiments.

Main paper: 
 introduction: 
 Optimal transport (OT) is now a widely used tool to compare probability distributions in different areas of data science such as machine learning [12, 25, 5] , computational biology [44] , imaging [20, 6] , even economics [22] or material sciences [8] . The computational and statistical efficiency of OT solvers is the key to facilitating their use in practical applications. Therefore, both computational methods and the statistical bottleneck in optimal transport (OT), often referred to as the curse of dimensionality, have received significant attention over the past decade [38, 53] . Regularization such as Entropic OT (EOT) [14] is a popular method to alleviate these two issues. It consists of adding an entropic regularization term to the objective function. Annealing schemes on the regularization parameter to approximate the true solution of OT by its entropic approximation are efficient, as shown in [31, 45, 19] . Still largely open is the theoretical understanding of these methods [48, 45] , which can shed light on the design of the annealing scheme, also called ε-scaling.OT and its entropic regularization apply to different contexts of interest. The most general context is when the two distributions are accessed via samples and one wants to estimate the OT distance and the correspondence plan or map. Another context of interest in some applications is the case of semi-discrete OT, as in [29] , in which one of the two distributions is discrete and the other continuous. This setting is slightly simpler than the general case since (i) the OT problem reduces to the estimation of Laguerre cells and (ii) the curse of dimensionality is alleviated [41] .Related works. In many applications of OT, one or both of the measures are accessed via i.i.d. samples. The goal then becomes to construct estimators of the OT map and/or cost. It is known that without any assumptions on the measures, the estimation of OT quantities suffers from the curse of dimensionality. For instance, estimating the Wasserstein distance from t samples achieves a rate of O(t -1 d ) for d ≥ 3. Despite the curse of dimensionality, estimating OT quantities attracts a lot of interest. Relevant works include [21, 53, 10, 42] for the OT cost, and [16, 27, 50, 41] for the OT map.We study here the estimation of OT quantities in the semi-discrete setting, where the continuous distribution is accessed through sampling, similarly to [33, 41] , but we assume full access to the discrete target measure, as in [24, 4] . This setting is of interest since, as recently shown in [41] , the OT map estimation escapes the curse of dimensionality, even without assuming the map to be smooth or continuous. Indeed, they showed that a rate of O(t -12 ) is achievable in the "one-sample" and "two-sample" settings (sampling only from the source measure or from both measures). To do so, their work uses the EOT map estimator [47, 40] with a regularization ε ≍ t -1 2 , as well as results on the convergence rate of the entropic optimal potential to the Kantorovich potential in the semi-discrete setting proved in [1, 17] . Moreover, they showed that the rate O(t -12 ) is minimax for the estimation of the OT map in the two-sample setting.Beyond the statistical challenges, building efficient solvers for semi-discrete OT is also a considerable challenge. Many solvers of (E)OT in this setting are based on optimizing the semi-dual, which is a finite-dimensional convex optimization problem. In particular, efficient Newton and quasi-Newton methods [34, 32, 29] are proposed for low dimensions, employing meshes when the source density is known. For arbitrary dimensions, or when the source measure is only accessible via samples, [24] propose using semi-dual EOT and Stochastic Gradient Descent (SGD) based solvers as proxies for OT. The study of SGD and Averaged SGD (ASGD) for EOT was further investigated by [4] , which notably demonstrated that the objective function is self-concordant and benefits from enhanced strong convexity near an optimum. Using these facts, [4] showed that a convergence rate of O(t -1 ) can be achieved for the squared Euclidean distance estimation of the discrete entropic optimal potential. However, terms in ε -1 were considered negligible in their study, thus excluding small regularization.Contributions. Our main contribution is twofold. First, we introduce an SGD-based algorithm to solve the semi-dual formulation of OT. This algorithm incorporates a projection step and an entropic regularization scheme that decreases with the number of samples. While being as computationally and memory efficient as vanilla SGD, our algorithm achieves enhanced convergence rates, thanks to the decreasing regularization. Specifically, given t i.i.d. samples of the source measure, it achieves a convergence rate of O(t -2b ) with b ∈ (1/2, 1) for both the discrete Kantorovich potential and OT cost estimation. We then construct an OT map estimator based on our discrete potential estimator and the closed form of the gradient of Fenchel transforms. By studying the difference between the Laguerre cells formed by the Kantorovich potential and our estimator, we retrieve a O(t -b ) convergence rate for the OT map, for b ∈ (1/2, 1).Second, building upon the parallel between measure support estimation and Laguerre cell estimations, we derive two new minimax lower bounds, characteristic of the fast rates of non-regular models: a O(t -2 ) rate for the Kantorovich potential and a O(t -1 ) rate for the OT map (compared to O(t -1/2 ) in the two-sample setting [41] ). These lower bounds are nearly achieved by our estimators since b < 1. Finally, we numerically showcase the convergence rates of our algorithm for the OT potential, map, and cost estimators. EQUATIONwhere Π(µ, ν) is the set of joints probability measures on R d × R d with marginals µ and ν. Mild conditions on µ, ν and the cost can be made so that this problem is well-posed, see [51] . When ε = 0, Problem (1) recovers the Kantorovich formulation of OT. In this article, we focus on the quadratic cost c(x, y) = 1 2 ∥x -y∥ 2 , although some of our results can be extended to more general costs. Our analysis relies on the semi-dual formulation of the convex problem (1) given byEQUATIONwhere for all y ∈ R d ,f c,ε (y) : = min x∈R d c(x, y) -f (x) if ε = 0, -ε log R d exp f (x)-c(x,y) ε dµ(x) if ε > 0.Under mild conditions on the cost or densities, a positive ε makes the semi-dual formulation ε -1smooth [15] . The key property of this semi-dual formulation of (E)OT is to retain more convexity than the standard dual of (1) (see [27, 49] ).Optimal maps and Brenier's theorem. We consider the quadratic cost, ε = 0 and µ, ν having second-order moments. Under the additional assumption that the measure µ is absolutely continuous, the optimal potential f * , called Kantorovich potential, is (locally) Lipschitz and the mapEQUATIONpushes forward µ onto ν (see [7] ). In addition, T µ,ν is the gradient of a convex function. This optimal map has more importance than the OT cost in subfields of machine learning such as generative modeling [28] or domain adaptation [13] . 
semi-discrete ot: 
 Semi-discrete (E)OT is when the source measure µ is absolutely continuous and the target measure ν = M j=1 w j δ yj is a finite sum of M ≥ 1 Dirac masses with weights w j > 0. In this case, the semi-dual formulation reduces to a finite-dimensional convex optimization problem onEQUATIONwhere we also writeH ε (g) = R d h ε (x, g)dµ(x), and where for all x ∈ R d , g c,ε (x) is a (vectorial) (c, ε)-transform with respect to a vector g = (g 1 , . . . , g M ) ∈ R M , defined byg c,ε (x) = min j∈ 1,M 1 2 ∥x -y j ∥ 2 -g j if ε = 0, -ε ln M j=1 exp -1 2 ∥x-yj ∥ 2 +gj ε w j if ε > 0.The vector g corresponds to the value of the potential function at the points y j . For all g ∈ R M and given X ∼ µ, an unbiased estimator of the gradient is given by∇ g h ε (X, g) j = -w j + χ ε j (X, g), 1 ≤ j ≤ M ,where for x ∈ R d , g ∈ R M , we haveχ ε j (x, g) = exp -1 2 ∥x-yj ∥ 2 +gj ε w j M k=1 exp -1 2 ∥x-y k ∥ 2 +g k ε w k • For ε = 0, χ j (x, g) = 1 Lj (g) (x) is an indicator function and we have a partition R d = M j=1 L j (g), where for all j ∈ 1, M ,L j (g) := x ∈ R d ; g c (x) = 1 2 ∥x -y j ∥ 2 -g j .The convex sets L j (g) are called power or Laguerre cells and µ(L i (g) ∩ L j (g)) = 0 when i ̸ = j. By the first-order optimality condition, solving semi-discrete OT amounts to finding g such that for all j ∈ 1, M , µ(L j (g)) = w j . Semi-discrete OT is a case of application of Brenier's theorem. Given the optimal potential g * , the OT map is T µ,ν (x) = x -∇(g * ) c (x) = xy j for x inside L j (g * ).Exploiting its finite-dimensional nature, optimizing the OT semi-dual H 0 has become a popular approach. Notably, Newton and quasi-Newton methods are highly effective in scenarios with low dimensions and known source densities, utilizing meshes to approximate the source density [34, 32, 29] . In scenarios involving arbitrary dimensions or when only sample-based access to the source measure is available, EOT emerges as a favored strategy. Notably, to avoid working with a discretized version of the source measure, such as with the Sinkhorn Algorithm, [24] recommend employing stochastic optimization to solve (4). Indeed, the semi-dual EOT problem has a convex objective of the formH ε (g) = E X∼µ [h ε (X, g)],with X as a random variable under µ. As noted in [24] , the main advantage of stochastic optimization algorithms is that they are suited for large-scale problems, keeping in memory only the discrete measure ν. Moreover, not relying on discretization permits an unbiased approach to solving the semi-discrete EOT problem.For a given fixed regularization parameter ε > 0, stochastic first-order methods are predominantly employed to solve (4). Starting with an initial value g 0 ∈ R M , these algorithms consider at each iteration one or many samples X t ∼ µ and rely on an update of the formg t = g t-1 -γ t ∇ g h ε (X t , g t-1).At time t, the Averaged Stochastic Gradient Descent (ASGD) returns the averaged estimate g t = 1 t+1 t k=0 g k , while Stochastic Gradient Descent (SGD) returns g t . ASGD, as an acceleration of SGD, has been widely studied in the literature (see [39, 37, 3] , and [4] for the specific case of EOT).Choosing the regularization parameter ε for EOT. Approximating the EOT problem rather than the OT one benefits from an enhanced convergence rate, especially in the discrete setting. The introduction of the Sinkhorn Algorithm for solving the EOT problem, as highlighted by [14] , has led to a resurgence of interest in OT within the machine learning community.The choice of the regularization parameter ε then becomes a practical and/or statistical problem:1. In the discrete case, selecting the regularization parameter is a practical issue that aims to strike an optimal balance between convergence speed and accuracy [14, 18] . To address this trade-off, some heuristics, such as ε-scaling [46] , which involves a decreasing regularization scheme, are employed, although they lack strong theoretical guarantees. 2. In the semi-discrete and continuous settings, the initial statistical problem is to determine the number of samples needed to accurately approximate the OT quantities. In this line of work, the use of EOT to construct estimators has also been proven to be satisfactory. In this case, studies show that regularization must decrease as the number of samples increases [40, 41] . However, discrete solvers do not adjust to the number of drawn points, as the solver is initiated once the points to approximate the measures have been sampled.3 DRAG: Decreasing Regularization Averaged Gradient 
setting: 
 We focus here on the one-sample setting of semi-discrete OT. Specifically, we sample from the source measure µ and leverage the full information of the discrete measure ν. Furthermore, fixing R > 0 and α ∈ (0, 1], we make the following mild assumption, already present in [17, 41] .Assumption 1. Let µ ∈ P(R d ), such that Supp(µ) ⊂ B(0, R) is convex and its density dµ is α-Hölder with, 0 < dµ < ∞ on its support. We note P α (B(0, R)) the set of these measures.The target measure ν is discrete, of the form ν = M j=1 w j δ yj , with w = (w 1 , . . . , w M ) its probability weights and (y 1 , . . . , y M ) ∈ B(0, R) M its support.Remark that with Assumption 1, we have that w min > 0 is lower-bounded.To accurately estimate the true OT cost and map, it is crucial to use a regularization parameter ε that decreases as the number of drawn samples increases. However, no existing algorithm in the OT literature simultaneously adapts to both entropic regularization and sample size. Inspired by the decreasing regularization scheme from the discrete OT setting, such as ε-annealing [46] , which is known for accelerating the convergence of the Sinkhorn algorithm in practice, and considering that SGD algorithms are inherently adaptive to the number of samples, we introduce the Decreasing Regularization projected Averaged Stochastic Gradient Descent (DRAG) to solve the semi-dual (2). Our algorithm employs a decreasing regularization sequence (ε t ) t and replaces the usual gradient step in ASGD with a projected step using adaptive regularizationg t = Proj C g t-1 -γ t ∇ g h εt-1 (X t , g t-1 ) ,where for U ⊂ R M convex, we define the projector as Proj U (g) := arg min{∥gg ′ ∥, g ′ ∈ U }. This method can be interpreted as a decreasing bias SGD scheme. For such a method, employing a projection step can be highly effective in ensuring convergence [11, 23] . In the context of EOT, it is well established that the (c, ε)-transform enables the localization of a ∥.∥ ∞ -ball, where a minimum of the semi-dual problem lies when the cost is bounded [36] . Specifically, since sup{c(x, y j ); x ∈ Supp(µ), j ∈ 1, M } < 2R 2 by Assumption 1, a preliminary projection set can be expressed as C ∞ := [0, 2R 2 ] M and we know that we can search for a minimum in this set. Nonetheless, leveraging the regularity of the cost function, we can have a projection set with a unique optimizer, as described in the following Lemma. Lemma 1. (Proof in Appendix B.7) Under Assumption 1, for all ε ≥ 0, there exists a unique solutiong * ε to (4) in C u := {g ∈ R M ; g 1 = 0 and |g j | ≤ R∥y 1 -y j ∥, j ∈ 1, M }. Algorithm 1 DRAG Parameters: (γ 1 , a, b, C) Initialize g 0 ∈ C, g 0 = g 0 , ε 0 = 1. for k = 1 to t do γ k = γ 1 k -b X k ∼ µ g k = Proj C g k-1 -γ k ∇ g h ε k-1 (X k , g k-1)g k = 1 k+1 g k + k k+1 g k-1 ε k = k -aend for return g t Note that the choice g 1 = 0 is arbitrary. In what follows, we refer to C = C ∞ or C = C u as our projection set. Note that for both sets, the projection's computational complexity is only O(M ), as it involves merely clipping each coordinate of our vector.Finally, in order to accelerate the convergence, we consider the Decreasing Regularization projected Averaged stochastic Gradient descent (DRAG) defined byg t = 1 t + 1 g t + t t + 1 g t-1 ,with g 0 = g 0 . The pseudo-code of our algorithm is given in Algorithm 1. A main advantage of DRAG is that it has a O(dtM ) computational complexity and O(dM ) spatial complexity.As a key step to the convergence rate of DRAG, we will provide the convergence rate of the non averaged estimate g t to g * εt , solving (4) with regularization ε t . Note that, up to a transformation of the form g * εt + a1 M , where a ∈ R * , the minimizer of the semi-dual is unique. Consequently, no matter the set C chosen, we focus our analysis on the orthogonal complement of the subspace spanned by 1 M , denoted as Vect(1 M ) ⊥ . For simplicity, for g, g ′ ∈ R M , we denote forp ∈ [1, ∞] ∥g -g ′ ∥ p := ∥g -g ′ ∥ p Vect(1 M ) ⊥ , ⟨g, g ′ ⟩ := ⟨g, g ′ ⟩ Vect(1 M ) ⊥ .Our analysis is greatly influenced by the findings in Corollary 2.2 from [17] , which states that for 0 ≤ ε ′ ≤ ε, under Assumption 1 with µ ∈ P α (B(0, R)), for any α ′ ∈ (0, α), there exists a constant K 0 , notably depending on the characteristics of ν (see [17] ), such thatEQUATIONIn addition, the convergence rates of our algorithm take advantage of the two following properties of the entropic semi-dual. For any ε > 0, noting w min := min j∈ 1,M w j ,• H ε is locally strongly convex on Vect(1 M ) ⊥ and the smallest eigenvalue of its Hessian atg * ε on Vect(1 M ) ⊥ is greater than w min ε -1 ([4], Lemma A.1). • H ε is 1 ε -self concordant ([4], Lemma A.2).Let us emphasize that, surprisingly, the first point reveals that the strong convexity at the optimum increases as we decrease the parameter ε. By combining these two points and benefiting from our projection step, we derive the following lemma.Lemma 2. (Proof in Appendix 2) For all regularization ε > 0 and for all g ∈ C, we haveEQUATIONMoreover, definingEQUATION)While technical, this lemma is a key step for our convergence guarantees and thus warrants further discussion. Note that A g,ε /K w can be interpreted as a form of local strong convexity coefficient of H ε . However, if ∥gg * ε ∥/ε is small, the term A g,ε tends to 0, and we would not be able to exploit more local strong convexity. This situation is unavoidable with any fixed regularization ε, if convergence to g * ε is desired. The use of a decreasing regularization scheme helps to avoid this problem. Indeed, if the term A gt,εt remains small for any t and ε t tends to 0, then ∥g tg * εt ∥ also tends to 0. However, if at time t, A gt,εt is close to 1, we can exploit strong convexity. Thus, a decreasing regularization scheme ensures good convergence behavior, regardless of A gt,εt . Building on these essential properties, we obtain the convergence rate for the non-averaged iterates of DRAG.Theorem 1. (Proof in Appendix B.1) Under Assumption 1 with µ ∈ P α (B(0, R)), taking the parameters (γ 1 , a, b) of DRAG such that γ 1 > 0, 1 + a + aα > 2b, a ≥ b 2 and b ∈ 1 2 , 1 , we have E ∥g t -g * εt ∥ 2p ≲ 1 w p min t bp , t ≥ 1, p ∈ {1, 2}.Remarkably, we achieve a convergence rate without any undesirable dependence on regularization. Our projection step and the improvements in Lemma 2, compared to Lemma A.1 in [4] , were crucial for this achievement. In contrast, [4] derived a convergence rate of the form O(ε -c t -b ) for a fixed regularization, with c at least equal to 1. Note that having no adverse dependence on the regularization parameter is essential for our algorithm, as it (i) employs a decreasing regularization scheme and (ii) aims to leverage the increased strong convexity at the optimum as ε t decreases. This last point will be further discussed in the next section.In convex stochastic optimization, it is known that averaging SGD iterations can lead to acceleration. More precisely, ASGD can adapt to the possibly unknown local strong convexity of the objective function at the optimizer [2] . As we saw previously, the strong convexity of H ε increases as the regularization parameter ε decreases. Despite the fact that our objective function changes at each time t, Theorem 2 (Proof in Appendix B.2) shows that DRAG fully exploits the increase in local strong convexity.Theorem 2. (Proof in Appendix B.2) Under the same assumptions as in Theorem 1, taking a ≥ b,E[∥g t -g * ∥ 2 ] ≲ 1 w 4 min t 2b , t ≥ 1 .Note that as b tends to 1, we achieve a quadratic convergence rate. We emphasize that this convergence rate is surprising, since for a general strongly convex function, the expected convergence rate would typically be linear. This difference comes from the fact that we face a Laguerre cells support problem.In parametric statistics, support problems are known to often be non regular and can yield an enhanced quadratic convergence rate (see, for instance, [52] , Chapter 15). In the next theorem, we show that our convergence rate to g * is nearly minimax.Theorem 3. (Proof in Appendix B.5) Let ν ∈ P(R) be a fixed discrete measure of M points. Then, inf g (t) sup µ∈Pα(B(0,R)) E ∥g (t) -g * ∥ 2 ≳ M t 2 ,where g * is the discrete optimal vector, solving the non regularized semi-dual in (4). The infimum is taken over all vectors g (t) ∈ R M constructed using t ∈ N * i.i.d samples of µ.Remark: While the dependence on w min (or M ) may seem minor in our context since it is a constant, we have included it in our analysis. This is pertinent, especially when applying DRAG to a discretized version of a continuous measure, which could result in a large M . We highlight that such results, demonstrating explicit dependence on the weights or number of points, are novel in the semi-discrete optimal transport (OT) literature. Additionally, when the weights of the discrete measure are uniform, our analysis achieves a convergence rate closer to O(M 2 t -2b ) (refer to the proof of Theorem 2 for further details). We believe that a theoretical convergence rate of O(w -2 min t -2b ) is achievable for DRAG. Indeed, a quadratic dependence on the strong-convexity coefficient is commonly observed in ASGD [2] . This dependence is illustrated in Figure 4 . 
optimal transport cost and brenier map estimation rate with drag 4.1 ot and eot cost estimation: 
 In this part, we derive convergence rates of the (E)OT costs using DRAG.Corollary 1. (Proof in Appendix B.3) Taking the same assumptions as Theorem 2, with 0 < ε < 1 and 0 < α ′ < 1, we have the following convergence rate for the approximation of the (E)OT costsEQUATIONEQUATIONOnce again, we achieve a superior rate compared to the typical O(t -1 ) observed in strongly convex and/or smooth scenarios, highlighting that semi-discrete OT deviates from conventional problems. While one could use a regularization parameter ε > 0, such as ε t , to approximate the OT cost, the best estimation is achieved when ε = 0. In this case, using H 0 (g t ) yields an error of the order O(t -2b ), a rate that is not achievable using (8) and ε > 0. This convergence rate is possible since H 0 is smooth, as noted in Theorem 4.1 of [30] .When employing entropic regularization, a popular choice involves using the estimator of the entropic Brenier map [40]EQUATIONIndeed, for ĝ ∈ R M , T ε µ,ν (ĝ)(x) could serve as an estimator. The objective is then to find an accurate estimator, ĝ, close to g * ε , and to analyze its performance based on the bias-variance decomposition 41] , Theorem 3.4). However, the mapping g → T ε µ,ν (g) is ε -1 -Lipschitz, complicating the bias-variance trade-off given that ε t = t -b . Instead, we rely on the gradient computed thanks to the c-transform of the estimator g t of DRAG. In fact, for any x ∈ R d , if there exists j ∈ 1, M such that x is in the interior of L j (g * ) ∩ L j (g t ), we have∥T µ,ν -T ε µ,ν (ĝ)∥ 2 L 2 (µ) ≲ ∥T ε µ,ν (ĝ) -T ε µ,ν (g * ε )∥ 2 L 2 (µ) + ε, using the fact that ∥T µ,ν -T ε µ,ν (g * ε )∥ 2 L 2 (µ) ≲ ε ([T µ,ν (x) = x -∇(g t ) c (x).Indeed, no matter g, as soon as x ∈ R d is in the interior of L j (g), the gradient of g c is given byEQUATIONBy analyzing the differences of Laguerre cells partitions between L(g t ) and L(g * ), we derive the following theorem. Theorem 4. (Proof in Appendix B.4) Under the same assumptions as Theorem 1, defining for all x ∈ R d and time t ≥ 0 T (g t )(x) = x -∇g c t , we have for all 1 ≤ p < ∞ the convergence rateE ∥T µ,ν -T µ,ν (g t )∥ p L p (µ) ≲ 1 t b .Minimax estimation. In the two-sample setting, where we subsample from both µ and ν, [41] shows that a convergence rate of O(t -1/2 ) is minimax for the squared L 2 error of the Brenier map estimation. As we see in Theorem 4, this rate can be improved to O(t -1 ) in the one-sample setting, as b tends to 1. In the following theorem, we prove that this rate is minimax. T (t) sup µ∈Pα(B(0,R)) E ∥T µ,ν -T (t) ∥ p L p (µ) ≳ 1 t ,where the infimum is taken over all maps T (t) constructed using t ∈ N * iid samples of µ. 
numerical experiments: 
 In this section, we numerically verify our convergence rate guarantees through various examples. For each example, we know the theoretical OT map, cost, and discrete potential. The first two examples are similar to those in [41] . In all figures, we fixed the parameters of DRAG to (γ 1 = √ w min , a = b = 0.75). While increasing b leads asymptotically to a better convergence rate, it decreases the step size of our gradient descent. Therefore, we need to wait longer to observe the acceleration from averaging. Our numerical investigation found that our parameter selection achieves a good compromise between convergence rate and the time before acceleration and is robust without further hypertuning.Examples settings: (1) µ ∼ U([0, 1] 10 ), Supp(ν) = {y j = ( j-1/2 J , 1 2 , ..., 1 2 ), j ∈ 1, 100 }, w = 1 100 1 100 . (2) µ ∼ U([0, 1] 10), M = 30 and y 1 , ..., y M randomly generated in [0, 1] 10 . We then also randomly generate g * ∈ R 30 and approximate w with Monte Carlo (MC) , such that g * is the discrete optimum potential. This setting led to w min = 0.00103.(3) µ ∼ U([δ, 1 + δ]), δ = 0.5, Supp(ν) = k M ; k ∈ 1, M , w = 1 M 1 M , M = 1000.While in dimension 1, this example is interesting since it appears in the proofs of Theorem 3 and 5. OT cost, map and potential convergence.Error ∝ t -b ∝ t -2b |H 0 (g t ) -H 0 (g * )| 1 |H 0 (g t ) -H 0 (g * )| 2 |H 0 (g t ) -H 0 (g * )| 3In Figures 1 and 2 , we show the convergence rates of the OT cost, map, and discrete potential. As we can see, we match our theoretical rates perfectly, except for the OT cost, where the rates are slightly slower. This discrepancy could be due to (i) our results for the OT cost estimations being asymptotic and (ii) our OT cost estimation already being extremely precise, with 10 7 MC samples proving insufficient to achieve precision around 5 • 10 -4 . Example 2∝ 1/t b ∝ 1/t 2b gt -g * 2 g t -g * 2Tµ,ν -T (ḡt) 2 Visualisation of the OT map estimators with DRAG.We visualize our OT map estimator T (g t ) =x -∇(g t ) c on a concrete example of Monge-Kantorovich (MK) quantiles [9] . In this context, having a target measure ν to investigate, the source measure is set to be the uniform measure on the unit Euclidean ball µ ∼ U(B(0, 1)). The goal is then to visualize the destinations through the OT map of points in regions B(0, (k + 1)/10) \ B(0, k/10) for k ∈ 0, 9 , which define MK quantile regions. We used M = 10 5 points to approximate ν, a discrete version of a boomerang-shaped measure. Finally, we launched DRAG with t iterations. In Figure 3 , we present the estimated MK quantiles regions of ν, where each color represents a region, starting from B(0, 0.1) in the center. In this example, taking t = 10 7 samples was sufficient and produced a similar result to when more samples were used. Error∝ 1/t g t -g * 2 1 gt -g * 2 1 g t -g * 2 3gtg * 2 3Figure 4 : Error evolution of DRAG for g t and g t as M grows, for Examples 1 and 3, t = M 2 iterations.Dependence of our convergence rate as M grows.As discussed in Section 3.4, our theoretical analysis indicates a dependence on w 4 min . In Examples 1 and 3, where similar problems arise with increased point counts, we run our algorithm with progressively larger M and M 2 iterations. Our theory predicts that the error of the estimator g t should decrease linearly. However, if the dependence of DRAG is indeed on w 4 min , the error of the estimator g t would remain constant or decrease linearly if the actual dependence is w 2 min . As illustrated in Figure 4 , our theoretical bound accurately matches the behavior of ∥g tg * ∥ 2 . Moreover, the behavior of ∥g tg * ∥ 2 suggests that our theoretical bound may not be sharp, as discussed after Theorem 2, since we observe a linear decrease.In the appendix, we present additional experiments that, while not altering our theoretical findings, could be highly beneficial for practitioners. Specifically, we provide evidence that mini-batching with GPU computation and weighted averaging of the iterates g t can significantly accelerate the algorithm. We also briefly discuss the choice of the parameters a and b and compare DRAG with SGD and ASGD. 
conclusion: 
 In EOT, a decreasing regularization parameter naturally appeals to practitioners who use annealing schemes to speed up Sinkhorn-like algorithms. Similarly, in the statistical community, a regularization that decreases with the number of points is preferred to more accurately approximate true OT quantities. With our algorithm, DRAG, we show that these two motivations for decreasing regularization can coexist successfully. To the best of our knowledge, this is the first algorithm in the OT literature that adapts to both the regularization and the sample size. Additionally, we derive two new minimax lower bounds for approximating OT quantities in the one-sample semi-discrete OT setting and show that DRAG nearly achieves these bounds.Our algorithm nearly achieves the minimax rate when b is close to 1. However, the closer b is to 1, the higher the constants in the rates. In practice, the choice a = b ≈ 0.75 gives robust practical results, as shown in Figure 7 in the appendix. An open direction is to design an improvement of our DRAG algorithm that achieves the minimax lower bound, while not suffering from large multiplicative constants, and remaining as computationally and memory efficient as our algorithm.Our results can also motivate further investigation into different lines of work: (i) Studying the convergence of the entropic discrete potential in semi-discrete OT for various costs. The key challenge in extending our algorithm's convergence proof is achieving results analogous to those in ([17], Corollary 2.2) for different cost functions. Alternatively, one could design an algorithm matching our bounds without entropic regularization, avoiding the need for such results. (ii) Developing decreasing regularization algorithms in the continuous case to efficiently approximate OT distances and maps. (iii) Adapting our approach to demonstrate or improve the acceleration of entropic annealing schemes for EOT solvers in the discrete case.Supplementary Material for "Semi-Discrete Optimal Transport: Nearly Minimax Estimation With Stochastic Gradient Descent and Adaptive Entropic Regularization"  
 
 Appendix: Weighted Averaging: Maintaining a better trade-off between averaged and non-averaged iterations. Since the dependence of DRAG iterates g t on the number of points M is at least quadratic, whereas for the non-averaged iterates g t it is only linear, when the total number of iterations t is insufficient (i.e., t ≤ M 2 ), g t can outperform g t as an estimator. One strategy to try to consistently achieve the best estimator regardless of the time t is through weighted averaging [35] .Namely, we replace the averaged estimator g t = 1 t+1 t k=0 g t , byg (ω) t := 1 t k=0 ln(k + 1) ω t k=0 ln(k + 1) ω g k ,with a parameter ω > 0. The parameter ω balances the weights assigned to the estimators g k . As ω increases, greater importance is given to the more recent estimates, while we retrieve g t when ω goes to 0. As for the usual averaged estimators, we can perform the weighted average online, without having to store all the iterates, with the recursiong (ω) t+1 = 1 - ln(t + 1) ω t k=0 ln(k + 1) ω g (ω) t + ln(t + 1) ω t k=0 ln(k + 1) ω g t+1 .It is important to note that g (ω) twill have the same asymptotic convergence guarantees as g t .In the following experiments, we operate under conditions where the number of iterations t is insufficient for the estimator g t to outperform g t . We set M = 1000 in Examples 1 and 3, select ω = 2 for the weighted average parameter, and fix t at 10 5 . Example 3 As illustrated in Figure 5 , the estimator g t begins to converge after approximately M iterations and remains superior to g t throughout the figure, since we are still within the regime where t ≤ M 2 . However, we see that the weighted average estimator g (ω) t consistently outperforms g t and already surpasses g t in performance after 10 5 iterations in Example 1.∝ 1/t b ∝ 1/t 2b gt -g * 2 g t -g * 2 g (ω) t -g * 2Mini-batch DRAG. As for Vanilla SGD, we can take advantage of GPU parallelization and replace the gradient estimator using one sample X ∼ µ ∇ g h ε (X, g) by a mini-batch estimator, using n b ≥ 1 i.i.d samples X 1 , ..., X n b samples of the source measure at onceEQUATIONOf course, no matter the choice n b , ( 12) defines an unbiased estimator of ∇H ε (g).Using a mini-batch of size n b , we suggest multiplying γ 1 by √ n b , as is usual with mini-batch SGD.The following figure shows the acceleration due to mini-batching in Example 2, while maintaining the same computational time when using a GPU. Indeed, each mini-batch estimator has an error an order of magnitude lower than the non-batched ones, even with a small mini-batch size of n b = 16. Error Influence of the parameter a and b. In Figure 7 , we illustrate the behavior of DRAG when changing the parameters a and b, on Example 2. As we can see, the choice a = b = 0.75 seems to be a good compromise on this experiments. We also see on Figure 7b that the non-averaged estimates with the best convergence rate is when b = 0.9. This behabiour is concordant with our theory. However, as we can see, the parameter b = 0.9 does not yet benefits from the acceleration thanks to averaging.gt -g * 2 ḡt -g * 2 ḡ(ω) t -g * 2 gt -g * 2 (Mini-batch) ḡt -g * 2 (Mini-batch) ḡ(ω) t -g * 2 (Mini-batch)DRAG compared with SGD and ASGD. We compare here the performance of our algorithm DRAG compared to the vanilla SGD and ASGD, introduced in [24] for EOT, on Example 2. For our comparison, since we fixed the parameters of DRAG to ( √ M , 3/4, 3/4) and ran the algorithm for t = 10 5 iterations, we have ε t = 10 -15/4 ≃ 10 -4 . We thus set ε = ε t to run SGD and ASGD. As we can see in Figure 8 , DRAG clearly outperforms SGD and ASGD. We note that the poor convergence of SGD and ASGD is not surprising with a small regularization parameter, as already observed, for instance, in [47] . Error gt -g * 2 , DRG (Non-Averaged) ḡt -g * 2 , DRAG gt -g * 2 , SGD, γ1 = √ M ḡt -g * 2 , ASGD, γ1 = √ M gt -g * 2 , SGD, γ1 = 1 ḡt -g * 2 , ASGD, γ1 = 1Additionnal notations for the proofs.For any c > 0 we define the function t → Ψ c (t) such thatEQUATIONFor a sequence (u t ) t∈N , if t 2 / ∈ N, u t 2 must be understood as u ⌈ t 2 ⌉ .B.1 Proof of Theorem 1: Convergence rate of the non averaged iterates.In all the sequel, we note∆ t = ∥g t -g * εt ∥ 2 .Remark that the dependence in t is both in the estimator g t and the optimizer g * εt . We also recall that we noteD C := sup g,g ′ ∈C ∥g -g ′ ∥ < ∞ .We will divide the proof into two parts.p = 1.Proof. By definition of the gradient step at time t + 1 and since g * εt+1 ∈ C, we have∆ t+1 = ∥g t+1 -g * εt+1 ∥ 2 = ∥Proj C (g t -γ t+1 ∇ g h εt (g t , X t+1 )) -g * εt+1 ∥ 2 ≤ ∥g t -γ t+1 ∇ g h εt (g t , X t+1 ) -g * εt+1 ∥ 2 .Then, incorporating the change of optimum between time t and t + 1, we get∆ t+1 ≤ ∥g t -γ t+1 ∇ g h εt (g t , X t+1 ) -g * εt + g * εt -g * εt+1 ∥ 2 ≤ ∥g t -γ t+1 ∇ g h εt (g t , X t+1 ) -g * εt ∥ 2 + 2 g t -γ t+1 ∇ g h εt (g t , X t+1 ) -g * εt , g * εt -g * εt+1 + ∥g * εt -g * εt+1 ∥ 2 .Using Corollary 2.2 in [17] , see (5) , there exists K 0 > 0 such that for any α ′ ∈]0, α[EQUATION)For clarity, we define r t := aK 0 t -(1+a+aα ′ ) and R t := (2D C + 2γ t+1 + r t )r t .Using that for all t, g t ∈ C, and that for allx ∈ R d , g ∈ R M , ∥∇ g h εt (g, x)∥ ≤ 2, we obtain ∆ t+1 ≤ ∥g t -γ t+1 ∇ g h εt (g t , X t+1 ) -g * εt ∥ 2 + (2D C + 2γ t+1 )∥g * εt -g * εt+1 ∥ + ∥g * εt -g * εt+1 ∥ 2 ≤ ∥g t -γ t+1 ∇ g h εt (g t , X t+1 ) -g * εt ∥ 2 + R t ≤ ∥g t -g * εt ∥ 2 -2γ t+1 ∇ g h εt (g t , X t+1 ), g t -g * εt + γ 2 t+1 ∥∇ g h εt (g t , X t+1 )∥ 2 + R t ≤ ∆ t -2γ t+1 ∇ g h εt (g t , X t+1 ), g t -g * εt + 4γ 2 t+1 + R t .Note that, since we have 1 + a + aα > 2b, we can also take α ′ ∈]0, α[ such that 1 + a + aα ′ > 2b.Therefore, the sequence R t /γ 2 t is decreasing and tends to 0. For conciseness, we noteEQUATION)For any t ≥ t a,α , we then obtain the following upper bound of ∆ t+1 in terms of ∆ t and the gradient direction:EQUATION)Noting F t the filtration generated by the samples X 1 , . . . , X t iid ∼ µ, that is F t = σ (X 1 , . . . , X t ) and taking the conditional expectation, we haveEQUATION)Using Lemma 2 and denoting A gt,εt = 1e -2 ε t [1∧∥gt-g * ε t ∥] , one has for all tEQUATIONThen, it comesEQUATION)We noteEQUATIONand note that λ t ≥ λ if ∥g tg * εt ∥ ∞ ≥ ε t . Therefore, we haveE [∆ t+1 | F t ] ≤ (1 -2λγ t+1 ) ∆ t + 2(λ -λ t )1 ∥gt-g * ε t ∥∞≤εt γ t+1 ∆ t + 5γ 2 t+1 . Moreover, ∥g t -g * εt ∥ ∞ ≤ ε t implies that ∆ t ≤ M ε 2 t . Therefore, E [∆ t+1 | F t ] ≤ (1 -2λγ t+1 ) ∆ t + 2(λ -λ t )1 ∥gt-g * ε t ∥∞≤εt γ t+1 M ε 2 t + 5γ 2 t+1 .Using that (λλ t )1 ∥gt-g * ε t ∥∞ ≤ λ and taking the expectation, we obtainE [∆ t+1 ] ≤ (1 -2λγ t+1 ) E[∆ t ] + 2λM ε 2 t γ t+1 + 5γ 2 t+1 .Noting t γ := min {t, 2λγ t+1 ≤ 1} and t 0 := max{t a,α , t γ } , we use Proposition 1 to obtainEQUATIONApplying Corollary 2, the exponential product converges exponentially to 0. Therefore, using the value of λ defined in (20) , an asymptotic comparison givesE[∆ t ] ≤ 5[2R 2 ∨ 1] w min (1 -e -2 ) γ t 2 -1 + M ε 2 t 2 -1 + O(γ 2 t ) .In the usual case where the discrete measure ν has uniform weights equal to 1 M , we deduce from the relation 2a ≥ b, by the assumption of the theorem, thatE[∆ t ] = O (M γ t ) . B.1.2 Part 2: proof for p = 2.Proof. Building on the proof of the case p = 1, we start by squaring equation (16) . For t ≥ t a,α , where t a,α is defined in (15) , we have∆ 2 t+1 ≤ ∆ t -2γ t+1 ∇ g h εt (g t , X t+1 ), g t -g * εt + 5γ 2 t+1 2 ≤ ∆ 2 t + 4γ 2 t+1 ∇ g h εt (g t , X t+1 ), g t -g * εt 2 + 25γ 4 t+1 -2∆ t γ t+1 ∇ g h εt (g t , X t+1 ), g t -g * εt =:A +5∆ t γ 2 t+1 -10γ 3 t+1 ∇ g h εt (g t , X t+1 ), g t -g * εt =:B .Taking the conditional expectation, recalling that λ t is defined in 18, we obtain thanks to Lemma 2E[A | F t ] ≥ 2∆ 2 t λ t γ t+1 .We also use the simple boundE[B | F t ] ≥ 0.These two inequalities lead toEQUATION)Using that the gradient norm is bounded by two, we use Cauchy-Schwarz inequality to obtain4γ 2 t+1 ∇ g h εt (g t , X t+1 ), g t -g * εt 2 ≤ 16γ 2 t+1 ∥g t -g * εt ∥ 2 ≤ 16∆ t γ 2 t+1 .Recalling the value of λ defined in (20) :EQUATIONwe use Hölder's inequality to obtain21∆ t γ 2 t+1 ≤ ∆ t √ 2λ 1 √ 2λ 21γ t+1 γ t+1 ≤ γ t+1 ∆ 2 t λ + 21 2 4λ γ 3 t+1 .Summing up these inequalities, we obtainE[∆ 2 t+1 | F t ] ≤ (1 -2λ t γ t+1 + λγ t+1 )∆ 2 t + 21 2 4λ γ 3 t+1 + 25γ 4 t+1 .Similarly to the case p = 1, we haveE[∆ 2 t+1 | F t ] ≤ (1 -2λγ t+1 + λγ t+1 )∆ 2 t + 2(λ -λ t )1 ∥gt-g * ε t ∥∞≤εt ∆ 2 t γ t+1 + 21 2 4λ γ 3 t+1 + 25γ 4 t+1 ≤ (1 -λγ t+1 )∆ 2 t + 2λM 2 ε 4 t γ t+1 + 21 2 4λ γ 3 t+1 + 25γ 4 t+1 .Taking the expectation, we obtainE[∆ 2 t+1 ] ≤ (1 -λγ t+1 )E[∆ 2 t ] + 2λM 2 ε 4 t γ t+1 + 21 2 4λ γ 3 t+1 + 25γ 4 t+1 .Proceeding as for the case p = 1, that is, applying Proposition 1 and Corollary 2 concludes the proof.Proof. We start by a decomposition of the gradient step, already present in [26] . By abuse of notation, we note∇ 2 k := ∇ 2 H ε k (g * ε k) and define the following differences:p k := Proj C (g k -γ k+1 ∇ g h ε k (g k , X k+1 )) -(g k -γ k+1 ∇ g h ε k (g k , X k+1 )) , ξ k+1 := ∇H ε k (g k ) -∇ g h ε k (g k , X k+1 ) , δ k := ∇H ε k (g k ) -∇ 2 k (g k -g * k ) .The term p k represents the difference between the projected and non-projected steps. Remark thatp k = 0 if g k -γ k+1 ∇ g h ε k (g k , X k+1 ) ∈ C.The difference of martingale ξ k represents the difference between the gradient and its unbiased version. Finally, δ k represents the difference between the gradient at g k with the linearized Hessian at the optimum.Noting I M the identity matrix of M M (R), observe that for any k ∈ Ng k+1 -g * ε k = Proj C (g k -γ k+1 ∇ g h ε k (g k , X k+1 )) -g * ε k = g k -γ k+1 ∇ g h ε k (g k , X k+1 ) -g * ε k -p k = g k -γ k+1 ∇H ε k (g k , X k+1 ) -g * ε k + γ k+1 ξ k+1 -p k = I M -γ k+1 ∇ 2 k (g k -g * k ) -γ k+1 δ k + γ k+1 ξ k+1 + p k .Thus, we have that∇ 2 k g k -g * ε k = g k -g k+1 γ k+1 -δ k + ξ k+1 + p k γ k .Observe that there is an orthogonal matrixU k such that ∇ 2 k = U k diag (λ k,1 , . . . , λ k,M -1 , 0) U ⊤ k .Therefore, in the following, we denote∇ 2 k -1 = U k diag λ -1 k,1 , . . . , λ -1 k,M -1 , 0 U ⊤ k the inverse of ∇ 2 k in the space Vect(1 M ) ⊥ . Note that we have ([4], Lemma A.1, equation (A.4)) min j∈ 1,M -1 λ k,j ≥ w min ε k , k ≥ 0.Taking all the equalities in Vect(1 M ) ⊥ , that is, considering all our vectors in the subspace Vect(1 M ) ⊥ , we haveg t -g * εt = 1 t + 1 t k=0 ∇ 2 k -1 g k -g k+1 γ k+1 :=L1,t - 1 t + 1 t k=0 ∇ 2 k -1 δ k :=L2,t + 1 t + 1 t k=0 ∇ 2 k -1 ξ k+1 :=Mt + 1 t + 1 t k=0 ∇ 2 k -1 p k γ k+1 :=L3,t + 1 t + 1 t k=0 (g * k -g * t ) =:Dt .Remark that the term 1 t+1 D t comes from the difference between 1 t+1 t k=0 g * ε k and g * εt . We will now bound the convergence rate for each of the sums in our decomposition. Note that the terms L 1,t , L 2,t and L 3,t will be, surprisingly, the limiting terms. Indeed, in stochastic optimization, M t is usually the main term. Nevertheless, the presence of the inverse of the Hessian ∇ 2 k -1 , whose largest eigenvalues is of order ε k , decreasing with k ≥ 1, makes it negligible.• Convergence rate for L 1,t . By the definition of our gradient step, we have1 t + 1 L 1,t = 1 t + 1 t k=0 ∇ 2 k -1 γ k+1 ∇ g h ε k (g k , x k+1 ) γ k+1 .Then, using that the gradient norm is bounded by 2, we obtain1 t + 1 E ∥L 1,t ∥ 2 1 2 ≤ 2 t + 1 t k=0 ∥(∇ 2 k ) -1 ∥ ≤ 2w -1 min t + 1 Ψ a (t).• Convergence rate for L 2,t . Using Lemma 2, for all k ≥ 0, we have∥δ k ∥ = ∇H ε k (g k ) -∇ 2 H ε k g * ε k g k -g * ε k ≤ 4 ε k ∥g k -g * ε k ∥ 2 ∞ .In addition, thanks to Theorem 1,E ∆ 2 k = O(w -2 min γ 2 k ). That is, there is a positive constant C 2 such that for all k ≥ 1, we have E ∆ 2 k ≤ C 2 w -2 min k -2b. Therefore,1 t + 1 E ∥L 2,t ∥ 2 1 2 ≤ 4w -1 min (t + 1) t k=0 E g k -g * ε k 4 ∞ ≤ 4w -1 min (t + 1) t k=0 E [∆ 2 k ] ≤ 4w -2 min √ C 2 (t + 1) Ψ b (t) = O(w -2 min t -b ) .Remark: When the weights are uniform, i.e., w min = 1/M , the bound can be of the order of Msmaller since ∥.∥ ∞ ≤ ∥.∥ ≤ √ M -1∥.∥ ∞ .Therefore, the bound can be closer to1 t + 1 E ∥L 2,t ∥ 2 1 2 = O(M t -b ).To emphasis this, we can fix β ∈ [0, 1] such that, when w = 1 M 1 M , we have1 t + 1 E ∥L 2,t ∥ 2 1 2 = O(M 1+β t -b ).• Convergence rate for L 3,t . In the same way as for L 1,t , we have that for any k∥p k ∥ ≤ 2γ k , such that 1 t + 1 E ∥L 3,t ∥ 2 1 2 ≤ t k=0 2∥(∇ 2 k ) -1 ∥ ≤ 2w -1 min Ψ a (t).However, we can retrieve a better convergence rate for this term.• Convergence rate for M t . Observe thatE[∥M t ∥ 2 ] = E ∥M t-1 ∥ 2 + 2 (∇ 2 t ) -1⊤ M t-1 , ξ t + ∥(∇ 2 t ) -1 ∥ 2 ∥ξ t ∥ 2 , with E (∇ 2 t ) -1⊤ M t-1 , ξ t = 0.Moreover, we have ∥ξ∥ ≤ 4, such that1 t + 1 E[∥M t ∥ 2 ] 1/2 = 4w -1 min t + 1 Ψ 2a (t) ≤ 2a 2a -1 4w -1 min t + 1 .• Convergence rate for D t . Thanks to (5), one as for all 0 < α ′ < α,1 t + 1 D t ≤ K 0 t + 1 t k=0 ε α ′ k (ε k -ε t ) ≤ K 0 t + 1 t k=0 ε 1+α ′ k ≤ K 0 t + 1 Ψ a+aα ′ (t),and this term is negligible since a + α ′ > b.• Conclusion. Taking a ≥ b as in the Theorem's assumption and summing up the inequalities, we obtainE ∥g t -g * εt ∥ 2 1 2 ≤ O(w 2 min t -b ) + o(t -b ). When w = 1 M 1 M , we obtain E ∥g t -g * εt ∥ 2 1 2 ≤ O(M 1+β t -b ) + o(t -b ).Using (5) , for any α ′ < α we have∥g * εt -g * ∥ ≤ K 0 ε 1+α ′ t ≤ K 0 t a+aα ′ = o(t -b ).Finally, we haveE ∥g t -g * ∥ 2 1 2 ≤ O(w 2 min t -b ) + o(t -b ),and whenw = 1 M 1 M , E ∥g t -g * ∥ 2 1 2 ≤ O(M 1+β t -b ) + o(t -b ).Remark: The main theorem considers a ≥ b to have the best convergence rate. However, note that from the proof, we can read the result when b/2 ≤ a < b. In this case, the limiting terms are only L 1,t and L 3,t .Proof. EOT cost estimation.For any ε > 0, the function H ε is 1 ε -smooth. Therefore, for any g ∈ R M , we haveH ε (g) -H ε (g * ε ) ≤ 1 2ε ∥g * ε -g∥ 2 .Using our estimator g t and (5), we obtainH ε (g t ) -H ε (g * ε ) ≤ 1 ε ∥g * ε -g * t ∥ 2 + ∥g t -g * t ∥ 2 ≲ ε 2α ′ -1 (ε -ε t ) 2 + 1 εt 2b .Remark. Using the triangular inequality and Theorem 2.3 in [17] , we also have|H 0 (g * ) -H ε (g t )| ≲ ε 2 + ε 2α ′ -1 (ε -ε t ) 2 + 1 εt 2b .OT cost estimation. For any vector g ∈ R M , we recall the definition ofL(g) = M j=1 L i (g) : for all j ∈ 1, M , L j (g) := x ∈ R d ; g c (x) = 1 2 ∥x -y j ∥ 2 2 -g j .Note that L(g) defines a partition, i.e. µ (L i (g) ∩ L j (g)) = 0 when i ̸ = j, and the convex sets L j (g) are called power or Laguerre cells. We define the setK δ := g : R M → R | ∀i ∈ 1, M , µ (L i (g)) > δ .Using Theorem 4.1 in [30] , under Assumption 1, H 0 is uniformly C 2,α on K δ . That is, there exists a constant L such that H 0 is L-smooth on K δ . Note that the constant L depends on µ min , δ, R. We refer to [30] , Remark 4.1 for more details.By the first order condition, as soon as δ ≤ w min , we have g * ∈ K δ . Indeed, at the optimum, we have for all i ∈ 1, M , L i (g * ) = w i . We fix here δ = 1 10 w min . Thanks to the L-smoothness, for any g ∈ K δ , we have|H 0 (g) -H 0 (g * )| ≤ L 2 ∥g -g * ∥ 2 .Note that, for any g ∈ R M and i ∈ 1, M , the difference of measure of the Laguerre cells L i (g) and L i (g * ) is at most linear with respect to ∥gg * ∥ ∞ . We refer to Theorem 4 or Section 6.4.2 in [43] for more details.Therefore, there exists a constant C L such that, as soon as ∥gg * ∥ 2 ≤ C L , we have that g ∈ K δ . This constant depends on δ, µ max , R and d as in Theorem 4. Using Theorem 2,E[∥g t -g * ∥ 2 ] = O(t -2b ). Then E [|H 0 (g t ) -H 0 (g * )|] = E |H 0 (g t ) -H 0 (g * )|1 g t ∈K δ + E |H 0 (g t ) -H 0 (g * )|1 g t / ∈K δ ≤ L 2 E[∥g t -g * ∥ 2 ] + max g∈C |H 0 (g) -H 0 (g * )|E[1 g t / ∈K δ ] ≤ L 2 E[∥g t -g * ∥ 2 ] + max g∈C |H 0 (g) -H 0 (g * )|E[1 ∥g t -g * ∥ 2 >C L ] = O(t -2b) , where the Markov inequality of order 1 was used on E[1 ∥g t -g * ∥ 2 >C L ].Proof. We will show here that a rate of convergence of g t to g * 0 gives a convergence rate for the map estimation. The Brenier map is equal to T µ,ν (x) = x -∇(g * 0 ) c (x); see for instance [43] , Theorem 1.17. We will thus focus on the convergence of ∇g c t to ∇(g * 0 ) c . For all j ∈ 1, M , if x is the interior of L j (g), we have ∇g c (x) = xy j .Therefore, given g, g ′ ∈ R M , if there exists a j ∈ 1, M such that x is the interior of L j (g) ∩ L j (g ′ ) we have ∇g c (x) = ∇(g ′ ) c (x). We will now follow arguments from [43] , Section 6.4.2. Fix j, j ′ ∈ 1, M such that j ̸ = j ′ and x is in the interior of L j (g) ∩ L j ′ (g ′ ). By definition of the c-transform, we can see that L j (g) is defined by M -1 linear inequalities of the form ⟨x, y j ′y j ⟩ ≤ a g (j, j ′ ) := g jg j ′ + 1 2∥y j ′ ∥ 2 2 - 1 2 ∥y j ∥ 2 2 .Similarly, interchanging the role of g, g ′ and j, j ′ we have ⟨x, y j -EQUATION) We have µ A := x ∈ R d , ∇g c (x) ̸ = ∇(g ′ ) c (x) = µ   j<j ′ L j (g) ∩ L j ′ (g ′ )   ≤ j<j ′ µ (L j (g) ∩ L j ′ (g ′ )) ≤ j<j ′ µ {x ∈ R d : -a g ′ (j ′ , j) ≤ ⟨x, y j ′ -y j ⟩ ≤ a g (j, j ′ )} .Under Assumption 1 µ is a measure such that Supp(µ) ⊂ B(0, R) and it admits a density dµ bounded by dµ max . Thus,µ(A) ≤ dµ max j<j ′ λ R d ({x ∈ B(0, R) : -a g ′ (j ′ , j) ≤ ⟨x, y j ′ -y j ⟩ ≤ a g (j, j ′ )}) ≤ dµ max j<j ′ λ R d x ∈ B(0, R) : - a g ′ (j ′ , j) ∥y j ′ -y j ∥ 2 ≤ x, y j ′ -y j ∥y j ′ -y j ∥ 2 ≤ a g (j, j ′ ) ∥y j ′ -y j ∥ 2 ≤ dµ max j<j ′ λ R d x ∈ B(0, R) : - a g ′ (j ′ , j) ∥y j ′ -y j ∥ 2 ≤ x 1 ≤ a g (j, j ′ ) ∥y j ′ -y j ∥ 2 ,by isotropy of the Lebesgue measure. Combining this remark with (24) yieldsµ(A) ≤ dµ max R d-1 j<j ′ |h j ′ -h j | ∥y j ′ -y j ∥ 2 .Similarly∥∥ (∇g c (•) -∇(g ′ ) c (•)∥ q ∥ p L p (µ) ≤ j<j ′ Lj (g)∩L j ′ (g ′ ) ∥ (∇g c (•) -∇(g ′ ) c (•)∥ q dµ(x) ≤ j<j ′ ∥y j ′ -y j ∥ q µ (L j (g) ∩ L j ′ (g ′ )) ≤ dµ max R d-1 j<j ′ ∥y j ′ -y j ∥ q |h j ′ -h j | ∥y j ′ -y j ∥ 2 ≤ dµ max M (2-q)+/2q R d-1 2M ∥h∥ 1 .So, in particular, there exists C ∆ > 0 independent of the location of the points y j but growing at least linearly in M such that∥∥ (∇g c (•) -∇(g ′ ) c (•)∥ q ∥ p L p (µ) ≤ C ∆ ∥g -g ′ ∥ 1 ≤ C ∆ √ M ∥g -g ′ ∥.Plugging the convergence rate of g t to g * concludes the proof.Proof. Let Θ ⊆ θ = (θ 1 , ..., θ M ) ∈ R M ; θ 1 = 0 and ν be a fixed discrete measure. For each θ ∈ Θ, consider ρ θ ∈ P α (B(0, R)) such that θ is the only vector in Θ for which the couple (θ c , θ) is solution of the dual of OT(ρ θ , ν).In our class of probabilities, the minimax estimation of the optimal transport potential θ, given t > 0 i.i.d samples of the source measure, can be written asR Θ M,t := inf θ(t) sup θ∈Θ E ρ θ ∥ θ(t) -θ∥ 2 ,where θ(t) is constructed with the t iid samples from the source measure µ. Note thatR Θ M,t ≤ inf g (t)sup µ∈Pα(B(0,R))EQUATION)where the infimum is taken over all vectors g (t) constructed with the t iid samples of µ.Let M ≥ 2 and take ν M the uniform measure on the points k M ; k ∈ 1, M . For δ ≥ 0, we note ρ θ δ ∼ U([δ, δ + 1]). Note that since d = 1, the optimal transport map is monotone non-decreasing (see, for instance, Chapter 2 in [43] ). Thus, for all k ∈ 1, M , we must have the identityT ρ θ δ ,ν (x) = k/M, x ∈ [δ + (k -1)/M ; δ + k/M ].Using the above information, the vector θ δ ∈ Θ is optimal for the semi-dual problem if and only if it satisfies the following inequalities for all k ∈ 1, M -1 ∀x ∈ [δ + (k -1)1/M, δ + k/M ] : θ δ,k+1 -θ δ,k ≤ - 1 M x + (2k + 1) 2 2M 2 , ∀x ∈ [δ + k/M, δ + (k + 1)/M ] : θ δ,k+1 -θ δ,k ≥ - 1 M x + (2k + 1) 2 2M 2 . For all k ∈ 1, M -1 , we thus obtain that θ δ,k+1 -θ δ,k = 1 2M 2 -δ 1 M .In particular, for any δ ≥ 0, we have∥θ 0 -θ δ ∥ 2 2 = M -1 k=1 kδ 1 M 2 = 1 M 2 δ 2 6 [M (M -1)(2M -1)] ≥ 1 6 δ 2 (M -1).Taking P, Q ∈ P(R d ) with densities ρ P and ρ Q , we recall that the Hellinger distance is defined byEQUATION.In particular, we have d H (ρ θ0 , ρ θ δ ) = √ 2δ. Applying Le Cam's Lemma (see, for instance [52] , Chapter 15) with δ = 1 8t givesR Θ M,t ≥ 1 4 1 - √ td H (ρ θ0 , ρ θ δ ) ∥θ 0 -θ δ ∥ 2 2 ≥ (M -1) 3072t 2 .Using the inequality (25) concludes the the proof.Proof. We fix the source measure ν =1 2 δ 1 + 1 2 δ 2 . For p ∈ [1, ∞[, we define Q 2,t := inf T (t) sup µ∈Pα(B(0,R)) E µ ∥∥ T (t) -T µ,ν ∥∥ p L p (U ([0,1])) ,where T (t) is constructed with t i.i.d samples from the source measure µ. Note that we have infEQUATIONWe define the family of source measures ρ δ = U([δ, 1 + δ]). Since the Brenier map is monotone increasing on the support of the source measure, we haveT δ (x) = 1, ∀x ∈ δ, 1 2 + δ , T δ (x) = 2, ∀x ∈ 1 2 + δ, 1 + δ .Fixing δ > 0, we see that∥T ρ θ 0 ,ν -T ρ θ 1 ,ν ∥ p L p (U ([0,1])) = δ.Using Le Cam's Lemma with δ = 1 8t , as in the proof of Theorem 3, we obtainEQUATIONUsing ( 27) concludes the proof.Proof. Following [36] , we know that an optimal couple of functions (f ε , g ε ) optimizing the dual formulation of EOT with regularization ε ≥ 0 satisfies the Schrödinger equations. That is, we can take for all y ∈ R d , g ε (y) = f c,ε ε (y). Moreover, 1 2 ∥x -y∥ 2 is R-Lipschitz on B(0, R). Therefore, since by Assumption 1, we have Supp(µ) ⊂ B(0, R) and Supp(ν) ⊂ B(0, R), we can exploit the Lipschitz property of our cost function on B(0, R). Using that the (c, ε)-transform has the same modulus of continuity as c (see Lemma 3.1 in [36] ), we get, for all y, y ′ ∈ R d :|f c,ε ε (y) -f c,ε ε (y ′ )| ≤ R∥y -y ′ ∥.That is, coming back to the function g, we have for all j, j ′ ∈ {1, ..., M } :ε (y j ) -g ε (y j ′ )| ≤ R ∥y j -y j ′ ∥ .By writing back our dual potential as a vector, that is g * = (g * 1 , . . . , g * M ), where for all j ∈ 1, M , g * j = g ε (y j ), we have|g * j -g * j ′ | ≤ R∥y j -y j ′ ∥.Moreover, if g * optimizes the semi-dual H ε , then for any β ∈ R, the vector g * + β1 M optimizes H ε .In particular, g *g ε (y 1 )1 M , which we rename g * , optimizes the semi-dual, with g * 1 = 0. Hence, for all j ∈ 1, ..., M g * y1g * yj = g * yj ≤ R∥y 1y j ∥. That is, there exists an optimizer in the desired closed convex set.Remark: Note that for other costs such as c(x, y) = ∥x -y∥ which defines the 1-Wasserstein distance, this projection set can be more relevant. Indeed, in this case, the cost is 1-Lipschitz and the projection set depends only on the target measure ν and no assumption of bounded cost is needed. In this case, the practitioner could choose the index k such that g k = 0, minimizing for instance the Euclidean diameter of the corresponding set.Proof. Since this proof heavily relies on Lemma A.2 in [4] , we will begin by rewriting the essential elements of their proof, using our notations, to derive our lemma. Note that in their proof, they study the concave problem -H ε (which they refer to as H ε ).Fix ε > 0 and g ∈ R M . Note g * ε ∈ Vect(1 M ) ⊥ such that min g∈R M H ε (g) = H ε (g * ε ). For any s ∈ [0, 1], denote g s = g * ε + s(g -g * ε ) and define the function φ : s ∈ [0, 1] → H ε (g s ).Following equation (A.21, [4] ), we haveEQUATIONwhere for all x ∈ R d and any s ∈ [0, 1], we define m(x, g s ) bym(x, g s ) := M j=1 χ ε j (x, g s )(g s -g * ε ).Instead of using Cauchy-Schwarz inequality as in [4] , we use Hölder's inequality with the Hölder conjugates p = 1, q = +∞ to obtainEQUATIONPlugging ( 29) in (28) givesEQUATIONThen, following from equation (A. 23, [4] ) to (A.27, [4] ) with our new inequality (29) leads to∇H ε (g) -∇ 2 H ε (g * ε ) (g -g * ε ) ≤ 2 ε ∥g -g * ε ∥ ∞ (φ(1) -φ(0)) ,where φ(0) = H ε (g * ε ) and φ(1) = H ε (g). Remark that for all ε > 0, H ε is 2-Lipschitz for the ∥.∥ ∞ norm. That is, we haveφ(1) -φ(0) = H ε (g) -H ε (g * ε ) ≤ 2∥g -g * ε ∥ ∞ .Therefore, we have the desired first bound in (6) .For the second bound in our proof, we still follow Lemma A.2 in [4] starting from line (A.28), with our new valueδ = 2 ε ∥g -g * ε ∥ ∞ ,such that using (30), we haveφ ′′′ (s) φ ′′ (s) ≥ -δ.Integrating between 0 and givesEQUATIONUsing that φ′′ (s) = (g -g * ε ) T ∇ 2 H ε (g s ) (g -g * ε ) and that the smallest eigenvalue of ∇ 2 H ε (g * ε ) is greater than w min /ε (Lemma A.1,[4]) implies that φ ′′ (0) ≥ w min ε ∥g -g * ε ∥ 2 .Then, using that φ ′ (s) = ⟨∇H ε (g s ), gg * ⟩ and integrating 31 between 0 and 1 gives⟨∇H ε (g), g -g * ε ⟩ ≥ w min ε 1 δ (1 -exp(-δ)) ∥g -g * ε ∥ 2 .Using a disjunction of cases, we obtain⟨∇H ε (g), g -g * ε ⟩ ≥    wmin ε ε 2 1 -exp -2∥g-g * ε ∥ ∞ ε ∥g -g * ε ∥ 2 ∞ if ∥g -g * ε ∥ ∞ ≤ 1, wmin ∥g-g * ε ∥ ∞ ε ε 2 1 -exp -2 ε ∥g -g * ε ∥ 2 ∞ if ∥g -g * ε ∥ ∞ ≥ 1.Then, using the projection step, no matter ifC = C ∞ or C = C u , we have sup x,y∈C {∥x -y∥ ∞ } ≤ 2R 2 .We thus have ∥g -g * ε ∥ ∞ ≤ 2R 2 , which leads to min w min ε ε 2 ; w min ∥g -g * ε ∥ ∞ ε ε 2 ≤ w min [2R 2 ∨ 1] 1 2 .Finally, noticing thatexp -2[∥g -g * ε ∥ ∞ ∧ 1] ε = exp -2∥g-g * ε ∥ ∞ ε if ∥g -g * ε ∥ ∞ ≤ 1, exp -2 ε if ∥g -g * ε ∥ ∞ ≥ 1,we obtain the desired bound.We can also derive a convergence rate without evaluating the regularized semi-dual nor using the unknown fixed smoothness of H 0 , noticing that the vectorial c-transform is non-expansive. That is, considering g 1 , g 2 ∈ R M , we have ∥g c 1g c 2 ∥ ∞ ≤ ∥g 1g 2 ∥ L ∞ (µ) . Theorem 6. Under the same assumptions as Theorem 2, we haveE∥f * -g c t ∥ 2 ∞ ≲ 1 t 2b ,which leads toE OT c (µ, ν) -g c t dν - M i=1 w i g i 2 ≲ 1 t 2b .Proof. By definition of the c-transform, for all x, y ∈ R d and all function g : R d → Rg c (x) + g j ≤ 1 2 ∥x -y j ∥ 2 .That is, for any f ∈ R M , we havef c (x) = inf j∈ 1,M 1 2 ∥x -y j ∥ 2 -f j ≥ g c (x) + inf j∈ 1,M [g j -f j ] ,such that we obtainEQUATIONTherefore, changing the role of f and g in (32) Proposition 1. Let (γ t ) t≥0 and (ν t ) t≥0 be some positive and decreasing sequences and let (δ t ) t≥0 , satisfying the following:• The sequence δ t follows the recursive relation:EQUATIONwith δ 0 ≥ 0 and ω > 0.• Let γ t converge to 0.• Let t 0 = inf {t ≥ 1 : ωγ t+1 ≤ 1}.Then, for all t ≥ t 0 , we have the upper bound:δ t ≤ exp -ω t i=t0+1 γ i t k=t0 γ k ν k + δ t0 + 1 ω ν ⌈ t 2 ⌉-1Proof. For all t ≥ t 0 , one hasδ t ≤ n i=t0+1(1ωγ i ) δ t0 =:U1,t+ t k=t0+1 t i=k+1 (1 -ωγ i ) γ k ν k =:U2,tOne can consider two cases: ⌈t/2⌉ -1 ≤ t 0 and ⌈t/2⌉ -1 > t 0 .Case where ⌈t/2⌉ -1 ≤ t 0 < t: Since ν k is decreasing,U 2,t ≤ ν t0+1 t k=t0+1 t i=k+1 (1 -ωγ i ) γ k = 1 ω ν t0+1 t k=t0+1 t i=k+1 (1 -ωγ i ) - t i=k (1 -ωγ i ) = 1 ω ν t0+1 1 - t i=t0+1 (1 -ωγ i ) ≤ 1 ω ν t0+1Since ν k is decreasing, it comes U 2,t ≤ 1 ω ν ⌈t/2⌉ . Case where ⌈t/2⌉ -1 > t 0 : As in [2] , for all m = t 0 + 1, . . . , t, one has Corollary 2. Let (γ t ) t≥0 and (ν t ) t≥0 be some positive and decreasing sequences and let (δ t ) t≥0 , satisfying the following:• The sequence δ t follows the recursive relation:EQUATIONwith δ 0 ≥ 0 and ω > 0.• Let γ t = c γ t -α with α ∈ (0, 1).• Let t 0 = inf {t ≥ 1 : ωγ t+1 ≤ 1}.Then, for all t ∈ N, we have the upper bound: δ t ≤ 


Decision: The paper deals with the problem of statistical estimation of the optimal transport (OT) map in the semi-discrete setting (i.e., the source measure is continuous and the target is discrete) with the target measure fully specified. When both the source and target measures are known only through samples, the minimax optimal rate of estimating the OT map is known to n^{-1/2}, where n is the number of samples from the source measure (as well as from the target measure). However, in the case where the target measure is fully known but the source measure is known from n iid samples, the paper shows that the OT map can be estimated at a minimax rate of n^{-1}. This is achieved by proposing an SGD-based algorithm to solve the semi-dual formulation of OT, which involves an entropic regularization scheme with a decreasing regularization parameter. All the reviewers are positive about the paper's novelty and contributions, which I agree with. Based on the discussion with the reviewers, I think the paper is a nice contribution to the statistical OT literature and therefore I am recommending acceptance. I strongly urge the authors to incorporate all the reviewer comments along with any clarifications they provided in the rebuttal into the camera-ready version.

Reviews: review: paper_summary: paper_summary: Optimal transport (OT) distances are a powerful computational tool to compare probability distributions and they play an increasingly preponderant role in machine learning, statistics, and computer vision. The comparison includes the cases of continuous-OT when the source and target probability distributions are continuous, discrete-OT for discrete distributions, and semi-discrete-OT for continuous source and discrete target. 

This paper studies the discrete-OT problem and provides a minimax lower bound for the OT map through the semi-dual formulation and using Laguerre cells and density support estimations. The authors propose the DRAG algorithm to solve the regularized semi-discrete OT that consists of an average stochastic gradient descent scheme coupled with a decreasing path on the regularization parameter $\varepsilon$.  Several theoretical guarantees are derived, specifically the Kantorovich potential that enjoys a convergence of $O(t^{-1})$ (see Theorem 1). Numerical experiments are conducted to corroborate the theoretical results on simulated data.
strength_weakness: strength_weakness: **Weakness and Questions**

- In Corollary 1, I am wondering why the convergence rate of the approximation cost is not with respect to the path $(\varepsilon_t)_t$. In the statement, it was for any $\varepsilon \in (0,1)$. Recall the DRAG algorithm is based on the $\varepsilon$-annealing scheme, then I think one can derive a bound for each $t$ as it was done in Theorem 1.
questions: questions: The notation $t$ is used for the sample size and the maximal number of iterations (see Algorithm 1 and the optimization convergence results).  I find some difficult to understand if the convergence bound is in term of sample size (that we have used to denote it by $n$ in OT-community) or the maximal number of iterations.
limitations: limitations: NA
paper_summary: paper_summary: The paper provides (almost tight) upper and lower bounds for computational optimal transport in the semi-discrete, one sample setting. The algorithm which is used for the upper bound is a stochastic gradient scheme, which solves the EOT estimation problem with a decreasing scheme for the regularization parameter.
strength_weakness: strength_weakness: The rates obtained are not exactly tight (a cost is paid as $b \to 1$, so there will always be a polynomial gap between upper and lower bounds in sample complexity. Do the authors conjecture that this can be removed, for instance by allowing an algorithm with exponential cost in computational complexity?
questions: questions: How flexible is this approach in handling different regularization schemes, or various gradient descent schemes which may have faster computational complexity? In particular, do the authors conjecture that the computational/memory costs in the paper are close to optimal in asymptotic order?

Typos:

82: "euclidean" -> "Euclidean"

158 "Setting." -> "Setting"

162 "H\"olderian" -> "H\"older"

206 "caracterstics" -> "characteristics"

777: "unknowned" -> "unknown"
limitations: limitations: All relevant limitations have been addressed.
paper_summary: paper_summary: This paper enhances our understanding of the semi-discrete optimal transport problem. Assuming full access to the discrete measure and only samples from the source, the authors provide a novel SGD-based algorithm based on entropic optimal transport. They go on to show that their estimator has more favorable properties than the naive estimator based on Sinkhorn's algorithm for estimating the potentials, and thus the map. They also prove matching lower-bounds that demonstrate that, in the one-sample setting, their result is tight. Experiments support their novel findings.
strength_weakness: strength_weakness: N/A but many small questions/comments/issues below
questions: questions: - Line 318 --- what exactly is $\nu$ here? 
- At the start of the appendix, $\ln$ and $\log$ are both used (maybe pick one).
- Typo in notations section "Euclidean" not "euclidean"
- Between lines 113 and 114 --- you could put $h_\varepsilon(x,g)$ in (4)
- Line 133: "really" is a tad informal
- Line 145: "Machine Learning community" should not be capitalized
- I do not understand what lines 172-173 mean in the larger context of the paragraph
- Lines 178-180: What does "localization of the potentials" mean?
- Line 188: "both set*s*"
- Line 189: Is the projection cost not $O(M)$?
- If you were the first to propose this DRAG algorithm, this should be highlighted more somehow (as written, it is not clear if you came up with it or if it already exists)
- Line 201: I presume $\mathbb{R}^*$ is $\mathbb{R}\setminus\{0\}$ but this should be mentioned somewhere
- Maybe add $w_{\min} = \min_{j \in [M]} w_j$ as an assumption, since this needs to always be lower-bounded, right?
- Maybe add a ref for (10)?
- Is it possible to include the dependence on "M" in the lower bound of Theorem 5?

Minor things I caught in the appendix:
- In Section B.1.1 --- this is incredibly minor, but I think it should be $R_t = (2D_{\mathcal{C}} + 4 \gamma_{t+1} + r_t)r_t$ 
- In B.2: The algorithm title seems to be different from DRAG
- Is $g_k^* = g_{\varepsilon_k}^*$? Otherwise I'm a confused
limitations: limitations: N/A

Comments: comment: I thank the authors for their prompt and thorough response, as well as for the suggested reference. I personally think this paper is strong and I am happy to recommend acceptance. I will elect to keep my score as it presently stands.