Pith. sign in

REVIEW 2 major objections 5 minor 31 references

Closing the Gap: Efficient Algorithms for Discrete Wasserstein Barycenters

T0 review · 2 major / 5 minor · reviewed 2026-08-04 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read The paper establishes the first polynomial-time approximation scheme for the discrete Wasserstein barycenter problem, achieving a (1+α)-approximation in time polynomial in (nk)^{1/α} and d.

desk verdict Solid core proof for a PTAS for discrete Wasserstein barycenters, but the runtime claim has an unproven ellipsoid step that needs fixing before the result is fully rigorous. read the letter →

arxiv 2511.04607 v2 pith:K6MDJ447 submitted 2025-11-06 math.OC

Jiaqi Wang , Weijun Xie This is my paper
classification math.OC MSC 49Q2268W2590C0890C25
keywords Wassersteinbarycenterdiscreteoptimaltransportmulti-marginalpolynomial-timeapproximationschemecandidatesupportellipsoidmethodalgorithmssample-meanvariancebound
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Computing the Wasserstein barycenter of k discrete measures—the probability measure minimizing a weighted sum of squared Wasserstein distances to all of them—is NP-hard, and the best previous polynomial-time guarantee was a factor of 2. This paper breaks that barrier by building a small candidate support set: sample t input measures according to the weights, take one support point from each sampled measure, average them, and solve a restricted multi-marginal optimal transport problem over those averages. Because each optimal barycenter atom is itself a weighted average of k input points, the random average of t independent draws has expected squared error at most (1/t) of the optimal fitting cost; this yields an expected (1+1/t)-approximation, and choosing t=⌈1/α⌉ gives a PTAS. Deterministic enumeration of all t-way averages gives the same ratio, and sampling without replacement improves it to 1+(k−t)/(t(k−1)) when the measures are equally weighted. The restricted problem is solved in polynomial time via a separation oracle, and experiments show the t=2 variant achieves near-optimal objectives on standard image datasets.

What carries the argument

Key machinery: candidate-support reduction powered by the sample-mean variance bound. The optimal barycenter support S* consists of weighted averages of every k-tuple of input points—exponentially many. The algorithms instead solve the restricted multi-marginal optimal transport dual over a candidate set S of t-way averages of sampled input points. Lemma 2 makes this tractable: feasibility of the dual is checked in O(kn|S|d) time by decomposing the minimum over index tuples into independent per-measure minimizations for each candidate support point, so the ellipsoid method runs in O~(n³k³d|S|). Lemma 3 recovers an explicit barycenter measure on S, and Theorem 1 converts pointwise support app

What would settle it

A direct check: for a small instance (e.g., k=3 measures with two support points each), enumerate all 2^3=8 candidate atoms to compute the exact optimum v*, then run the deterministic t=2 enumeration and test whether the restricted objective v(S_2^1) ≤ 1.5 v*; any violation refutes Theorem 3. Separately, instrument the ellipsoid solver on random instances of increasing n,k to verify that separation-call counts and bit-lengths match the claimed soft-O polynomial bound.

Watch

Extended reading notes

Core claim

The central discovery is a support-reduction inequality: for any k-tuple of input support points with weighted center c, if T is an i.i.d. t-sample from the weight vector, the expected squared distance from c to the best candidate average in S_T is at most (1/t) Σ_i λ_i‖x_i−c‖². Feeding this inequality into the multi-marginal optimal transport dual shows that the restricted barycenter over S_T has expected cost at most (1+1/t) v*. Deterministic enumeration gives the same ratio, and sampling without repetition under equal weights improves it to 1+(k−t)/(t(k−1)). Choosing t=⌈1/α⌉ returns a (1+α)-approximate barycenter in time polynomial in (nk)^{1/α} and d.

Load-bearing premise

The load-bearing premise is that the restricted multi-marginal optimal transport dual over the candidate support can be solved in polynomial time by the ellipsoid method with the closed-form separation oracle—that is, the bounding box for the dual variables, the number of separation calls, and the arithmetic precision all stay within the claimed polynomial bounds, which the paper invokes but does not explicitly construct.

Editorial extensions

If this is right

  • For any prescribed α∈(0,1], a (1+α)-approximate discrete Wasserstein barycenter can be found in time polynomial in (nk)^{1/α} and d—the first PTAS for this problem in general dimension.
  • The previous factor-2 approximation is recovered as the t=1 special case, so the method is a tunable accuracy–runtime tradeoff rather than a separate algorithm.
  • Under equal weights, the without-replacement variant attains the strictly better ratio 1+(k−t)/(t(k−1)), which for t>1 is smaller than 1+1/t.
  • Because the guarantee holds for every fixed input tuple, the deterministic enumeration version is fully deterministic and gives the same approximation ratio with support size O((nk)^t).
  • Numerical experiments on ellipse, digit, and sign-language benchmarks show the t=2 deterministic variant reaching near-optimal objectives in instances where exact multi-marginal formulations would require on the order of 10^22 to 10^120 operations.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The paper leaves implicit that the t-averaging construction is a general variance-reduction device: any optimization over probability measures whose objective is a convex average of squared distances could carry the same (1+1/t) guarantee whenever optimal atoms are convex combinations of input points—e.g., Wasserstein regression or distributionally robust location.
  • A testable extension is to replace i.i.d. sampling with stratified or low-discrepancy sampling of indices; the variance bound suggests the effective constant could move from (1+1/t) toward 1 at lower enumeration cost than (nk)^t.
  • For non-uniform weights, a without-replacement scheme drawn from the weight distribution should yield an intermediate ratio between 1+1/t and 1+(k−t)/(t(k−1)); the paper does not analyze this hybrid.
  • The hybrid candidate-support strategy used in the large-support experiments suggests a practical route to near-exact barycenters: use the PTAS support to warm-start column generation or exact solvers, potentially inheriting both the theoretical guarantee and practical scalability.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

2 major / 5 minor

Summary. The paper studies the discrete Wasserstein barycenter problem and proposes approximation algorithms that improve the known factor-2 guarantee. For general weights, the authors construct a candidate barycenter support by averaging support points from t measures sampled i.i.d. according to the weights (randomized) or by enumerating all such averages (deterministic), and prove that the restricted barycenter value is at most (1+1/t) v* by an elementary variance decomposition. For equally weighted measures, sampling without replacement yields the improved factor 1+(k-t)/(t(k-1)). With t = ceil(1/alpha), this gives a (1+alpha)-approximation. The restricted MOT dual is solved via the ellipsoid method with a separation oracle (Lemma 2). Numerical experiments on synthetic and real data illustrate the trade-off. The paper's central approximation-ratio arguments are correct and parameter-free, but the computational engine—Lemma 2—has a significant omitted justification.

Significance. If the claimed PTAS is fully established, it is a genuine advance: it breaks the factor-2 barrier for discrete Wasserstein barycenters in general dimension, matching the NP-hardness result with a polynomial-time approximation scheme whose exponent depends on 1/alpha. The proof of the approximation guarantee is simple, transparent, and does not rely on any fitted parameters or circular assumptions. The deterministic variants are a plus, and the equal-weight improvement is a nice refinement. However, the polynomial-time claim rests on Lemma 2, whose ellipsoid-method invocation omits standard but necessary boundedness and bit-length conditions. This is a fixable gap, but it is load-bearing for the computational half of the main claim.

major comments (2)
  1. [Section 2.2, Lemma 2] The lemma asserts that the restricted MOT dual (3) can be solved in O~(n^3 k^3 d |S|) time via the ellipsoid method. The separation oracle itself is correctly described, but the ellipsoid invocation is incomplete. The dual feasible region is unbounded: adding constants a_i to all variables gamma_{i,j} of marginal i leaves feasibility and objective unchanged whenever sum_i a_i = 0. No gauge fixing, bounding box, or initial ellipsoid is provided. Moreover, the ellipsoid iteration count depends on log(R/r) and on the bit-length of the data; these are suppressed but not bounded. Without these ingredients, the claimed O~(n^3 k^3 d |S|) bound is not justified. Since the PTAS runtime in the abstract and in Theorem 3 relies on Lemma 2, this is a load-bearing gap. A direct fix is to solve the fixed-support barycenter LP with O(k n |S|) variables, which is polynomial size and avoids ellipsoid subt
  2. [Section 2.2 and Algorithm 1, step 4] The algorithms state 'Solve the restricted MOT dual (3) and obtain its corresponding primal optimal solution Pi*(S)' and then invoke Lemma 3 to output a barycenter. Lemma 2, however, only yields an optimal dual solution (or value) via the ellipsoid method; it does not explain how to extract an optimal primal solution of the exponentially sized MOT primal (5). Since the output is an explicit probability measure, this recovery step is part of the algorithm. The paper should either solve the equivalent polynomial-size fixed-support LP directly, or describe a valid primal-recovery procedure from the dual solution (e.g., via complementary slackness with a guarantee of polynomial complexity).
minor comments (5)
  1. [Algorithm 1, line 2] The line reads 'Uniformly sample a size-t multi-subset T={T_i}_{i∈[k]} from subset [k]', but the construction in Section 3.1 requires i.i.d. draws T_1,...,T_t with Prob[T_j=i]=lambda_i, and the index set should be i∈[t], not i∈[k]. This is a typo, but it is confusing because the algorithm as printed is not the one analyzed in Theorem 2.
  2. [Theorem 5 and Section 4.2] Theorem 5 references 'Algorithm 3' but the deterministic counterpart is Algorithm 4. Its statement also writes 'min_{T:|T|=t} v(S_t^2)', which is not meaningful because S_t^2 is already the union over all subsets; it should simply be v(S_t^2) ≤ ... . These are presentation issues, not mathematical errors.
  3. [Theorem 3 proof] In the first displayed inequality of the proof, the support is written as S_t^2, but this theorem concerns S_t^1. The subsequent logic is correct, so this is only a typo.
  4. [Theorem 4 proof] The final displayed equality after 'Plugging (7d) to (7c)' drops the factor 1/k on the right-hand side: the result should be (1 + (k-t)/(t(k-1))) * (1/k) sum_j ||c-x_j||^2. The derivation immediately before it is correct, so the missing factor appears to be a typographical error.
  5. [Section 5] The numerical section alternates between 'Algorithm 2' and 'Algorithm 4' (e.g., Section 5.2 says Algorithm 2, while Section 5.3 says Algorithm 4) without stating the weights used. Since the equal-weight and general-weight guarantees differ, the experimental setup should state which algorithm and weight regime applies to each dataset.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: the PTAS guarantee is derived from an elementary variance identity, with no fitted inputs or load-bearing self-citations.

full rationale

This paper's central derivation chain is self-contained. Theorem 2's (1+1/t) guarantee follows from condition (4), which is proven directly in equations (6a)-(6e) using the variance identity E||(1/t)Σ x_{T_j} - c||² = (1/t)Σλ_i||x_i-c||² for i.i.d. samples with probabilities λ_i; the target approximation factor is the output of this computation, not an input. Theorem 1 then transfers the support-approximation condition to the restricted MOT value via the inequality E_S[v(˜S)] ≤ E_S⟨C(˜S),Π*⟩ ≤ (1+α)v(S*), with Π* an optimal full-MOT solution feasible for the restricted problem; this is a standard domination argument, not a definitional identity. Theorems 4-5 use the same argument with sampling without replacement and a correct covariance computation. There are no fitted parameters, no data-dependent tuning invoked in the proof, and no load-bearing self-citations; cited external results (MOT equivalence, ellipsoid method, hardness) are used as background and are not substitutes for the main argument. The only notable weakness is Lemma 2's ellipsoid invocation: the feasible set of dual variables is defined only up to per-marginal constants and no explicit bounding box/R ratio is supplied, so the claimed O~(n^3k^3d|S|) runtime is not fully justified as written. That is a correctness/complexity gap, not circularity, because the approximation-ratio theorems do not depend on that solver bound. Accordingly no circular step is established; score 0.

Assumptions & free parameters 0 free parameters · 4 assumptions · 0 invented entities

No numbers are fitted: t is a user-selected accuracy parameter, not a fitted value. The proof relies on standard LP duality, the MOT equivalence, and ellipsoid separation complexity. No new entities are introduced.

assumptions (4)
  • standard math Equivalence between (WBCenter) and the multi-marginal optimal transport LP, and the characterization of optimal support S* as weighted least-squares minimizers over k-tuples (Lemma 1).
    Imported from Altschuler & Boix-Adsera (2022, Prop 2.1) and Carlier & Ekeland (2010, Sec 6); used to motivate support reduction.
  • standard math LP strong duality for the MOT primal/dual pair (2)-(3).
    Used to equate v(S) with the restricted primal optimum and to recover a barycenter in Lemma 3.
  • standard math Ellipsoid method solves LPs with O(nk) variables and a polynomial-time separation oracle in O~((nk)^2) iterations.
    Used in Lemma 2 to claim polynomial solvability of the restricted dual; assumes rational data and exact arithmetic.
  • domain assumption Problem setting: finite-support measures with at most n points, weights λ, objective is weighted sum of squared type-2 Wasserstein distances.
    This is the problem definition from Section 1; the approximation ratio and runtime are relative to this objective.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Closing the Gap: Efficient Algorithms for Discrete Wasserstein Barycenters." pith.science (2026). https://pith.science/paper/K6MDJ447

@misc{pith2026251104607,
  author       = {Pith},
  title        = {Pith review of: Closing the Gap: Efficient Algorithms for Discrete Wasserstein Barycenters},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/K6MDJ447}},
  note         = {Machine review of arXiv:2511.04607}
}
abstract

Given a collection of probability measures, the Wasserstein barycenter problem seeks a probability measure that minimizes the weighted sum of their squared type-$2$ Wasserstein distances. We study the discrete setting, in which every input measure has finite support; this setting arises frequently in machine learning and operations research. Since computing a discrete Wasserstein barycenter is NP-hard, we develop approximation algorithms with provable performance guarantees. The best previously known polynomial-time guarantee is a factor of $2$. We improve this bound by proposing a polynomial-time approximation scheme (PTAS) that generalizes the existing $2$-approximation: for any prescribed $\alpha\in(0,1]$, the proposed algorithms return a $(1+\alpha)$-approximate barycenter in time polynomial in $(nk)^{1/\alpha}$ and $d$. We provide both randomized and deterministic constructions and derive a sharper guarantee when the input measures are equally weighted. Numerical experiments on synthetic and real datasets show that the algorithms are computationally practical and produce near-optimal barycenters.

Figures

Figures reproduced from arXiv: 2511.04607 by the authors.

Figure 1
Figure 1. Ten images from the nested ellipses dataset. [PITH_FULL_IMAGE:figures/full_fig_p016_1.png] view at source ↗
Figure 2
Figure 2. Barycenters produced by different algorithms for nested ellipse. [PITH_FULL_IMAGE:figures/full_fig_p017_2.png] view at source ↗
Figure 3
Figure 3. Barycenters produced by different algorithms for MNIST. [PITH_FULL_IMAGE:figures/full_fig_p018_3.png] view at source ↗
Figures from the paper (3 more)
Figure 4
Figure 4. Figure 4: Barycenters produced by different algorithms for gesture 1. [PITH_FULL_IMAGE:figures/full_fig_p019_4.png]
Figure 5
Figure 5. Figure 5: Barycenters produced by different algorithms for gesture 2. [PITH_FULL_IMAGE:figures/full_fig_p019_5.png]
Figure 6
Figure 6. Figure 6: Barycenters produced by different algorithms for gesture 3. [PITH_FULL_IMAGE:figures/full_fig_p020_6.png]

Discussion (0). Sign in to comment.

Reference graph

Works this paper leans on

31 extracted references · 2 linked inside Pith

  1. [1]

    & Carlier, G

    Agueh, M. & Carlier, G. (2011), ‘Barycenters in the wasserstein space’,SIAM Journal on Mathe- matical Analysis43(2), 904–924

    2011

  2. [2]

    Altschuler, J. M. &Boix-Adsera, E. (2021), ‘Wasserstein barycenters can be computedin polynomial 20 time in fixed dimension’,Journal of Machine Learning Research22(44), 1–19

    2021

  3. [3]

    Altschuler, J. M. & Boix-Adsera, E. (2022), ‘Wasserstein barycenters are np-hard to compute’, SIAM Journal on Mathematics of Data Science4(1), 179–203

    2022

  4. [4]

    & Miller, J

    Anderes, E., Borgwardt, S. & Miller, J. (2016), ‘Discrete wasserstein barycenters: Optimal transport for discrete data’,Mathematical Methods of Operations Research84(2), 389–409

    2016

  5. [5]

    & Peyré, G

    Benamou, J.-D., Carlier, G., Cuturi, M., Nenna, L. & Peyré, G. (2015), ‘Iterative bregman projections for regularized transportation problems’,SIAM Journal on Scientific Computing 37(2), A1111–A1138

    2015

  6. [6]

    Bigot, J., Klein, T. et al. (2012), ‘Consistent estimation of a population barycenter in the wasserstein space’,ArXiv e-prints

    2012

  7. [7]

    & Cuturi, M

    Bonneel, N., Peyré, G. & Cuturi, M. (2016), ‘Wasserstein barycentric coordinates: histogram re- gression using optimal transport.’,ACM Trans. Graph.35(4), 71–1

    2016

  8. [8]

    & Pfister, H

    Bonneel, N., Rabin, J., Peyré, G. & Pfister, H. (2015), ‘Sliced and radon wasserstein barycenters of measures’,Journal of Mathematical Imaging and Vision51(1), 22–45

    2015

Show all 31 references
  1. [9]

    (2022), ‘An lp-based, strongly-polynomial 2-approximation algorithm for sparse wasserstein barycenters’,Operational Research22(2), 1511–1551

    Borgwardt, S. (2022), ‘An lp-based, strongly-polynomial 2-approximation algorithm for sparse wasserstein barycenters’,Operational Research22(2), 1511–1551

    2022

  2. [10]

    & Patterson, S

    Borgwardt, S. & Patterson, S. (2020), ‘Improved linear programs for discrete barycenters’,Informs Journal on Optimization2(1), 14–33

    2020

  3. [11]

    & Patterson, S

    Borgwardt, S. & Patterson, S. (2021), ‘On the computational complexity of finding a sparse wasser- stein barycenter’,Journal of Combinatorial Optimization41(3), 736–761

    2021

  4. [12]

    & Patterson, S

    Borgwardt, S. & Patterson, S. (2022), ‘A column generation approach to the discrete barycenter problem’,Discrete Optimization43, 100674

    2022

  5. [13]

    & Ekeland, I

    Carlier, G. & Ekeland, I. (2010), ‘Matching for teams’,Economic Theory42(2), 397–418

    2010

  6. [14]

    Cheng, K., Aeron, S., Hughes, M. C. & Miller, E. L. (2021), ‘Dynamical wasserstein barycenters for time-series modeling’,Advances in Neural Information Processing Systems34, 27991–28003

    2021

  7. [15]

    & Solomon, J

    Claici, S., Chien, E. & Solomon, J. (2018), Stochastic wasserstein barycenters,in‘International Conference on Machine Learning’, PMLR, pp. 999–1008

    2018

  8. [16]

    & Doucet, A

    Cuturi, M. & Doucet, A. (2014), Fast computation of wasserstein barycenters,in‘International conference on machine learning’, PMLR, pp. 685–693

    2014

  9. [17]

    Ge, D., Wang, H., Xiong, Z. & Ye, Y. (2019), ‘Interior-point methods strike back: Solving the wasserstein barycenter problem’,Advances in neural information processing systems. 21 Grötschel, M., Lovász, L. & Schrijver, A. (2012),Geometric algorithms and combinatorial optimiza-...

    2019

  10. [18]

    & Lai, L

    Huang, M., Ma, S. & Lai, L. (2021), Projection robust wasserstein barycenters,in‘International Conference on Machine Learning’, PMLR, pp. 4456–4465

    2021

  11. [19]

    & Uribe, C

    Kroshnin, A., Tupitsa, N., Dvinskikh, D., Dvurechensky, P., Gasnikov, A. & Uribe, C. (2019), On the complexity of approximating wasserstein barycenters,in‘International conference on machine learning’, PMLR, pp. 3530–3540

    2019

  12. [20]

    Lau, T. T.-K. & Liu, H. (2022), ‘Wasserstein distributionally robust optimization with wasserstein barycenters’,arXiv preprint arXiv:2203.12136

    2022 arXiv

  13. [21]

    & Haffner, P

    LeCun, Y., Bottou, L., Bengio, Y. & Haffner, P. (2002), ‘Gradient-based learning applied to docu- ment recognition’,Proceedings of the IEEE86(11), 2278–2324

    2002

  14. [22]

    & Jordan, M

    Lin, T., Ho, N., Chen, X., Cuturi, M. & Jordan, M. (2020), ‘Fixed-support wasserstein barycenters: Computational hardness and fast algorithm’,Advances in neural information processing systems 33, 5368–5380

    2020

  15. [23]

    & Ruszczyński, A

    Lin, Z. & Ruszczyński, A. (2025), ‘Federated calculation of the free-support transportation barycen- ter by single-loop dual decomposition’,arXiv preprint arXiv:2507.19627

    2025 arXiv

  16. [24]

    Lindheim, J. v. (2023), ‘Simple approximative algorithms for free-support wasserstein barycenters’, Computational Optimization and Applications85(1), 213–246

    2023

  17. [25]

    & Ciliberto, C

    Luise, G., Salzo, S., Pontil, M. & Ciliberto, C. (2019), ‘Sinkhorn barycenters with free support via frank-wolfe algorithm’,Advances in neural information processing systems

    2019

  18. [26]

    A., Heitz, M., Bonneel, N., Ngole, F., Coeurjolly, D., Cuturi, M., Peyré, G

    Schmitz, M. A., Heitz, M., Bonneel, N., Ngole, F., Coeurjolly, D., Cuturi, M., Peyré, G. & Starck, J.-L. (2018), ‘Wasserstein dictionary learning: Optimal transport-based unsupervised nonlinear dictionary learning’,SIAM Journal on Imaging Sciences11(1), 643–678

    2018

  19. [27]

    & Aberdam, A

    Simon, D. & Aberdam, A. (2020), Barycenters of natural images constrained wasserstein barycenters for image morphing,in‘Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition’, pp. 7910–7919

    2020

  20. [28]

    A., Dvinskikh, D., Dvurechensky, P., Gasnikov, A

    Uribe, C. A., Dvinskikh, D., Dvurechensky, P., Gasnikov, A. & Nedić, A. (2018), Distributed computation of wasserstein barycenters over networks,in‘2018 IEEE Conference on Decision and Control (CDC)’, IEEE, pp. 6544–6549. 22

    2018

  21. [29]

    & Chateau, T

    Yan, Y., Duffner, S., Phutane, P., Berthelier, A., Blanc, C., Garcia, C. & Chateau, T. (2021), ‘2d wasserstein loss for robust facial landmark detection’,Pattern Recognition116, 107945

    2021

  22. [30]

    & Toh, K.-C

    Yang, L., Li, J., Sun, D. & Toh, K.-C. (2021), ‘A fast globally linearly convergent algorithm for the computation of wasserstein barycenters’,Journal of Machine Learning Research22(21), 1–37

    2021

  23. [31]

    Ye, J., Wu, P., Wang, J. Z. & Li, J. (2017), ‘Fast discrete distribution clustering using wasserstein barycenter with sparse support’,IEEE Transactions on Signal Processing65(9), 2317–2332. 23

    2017

Pith tools

Reviewed August 4, 2026 · model on record in the stance chip above.