Wasserstein Coreset via Sinkhorn Loss
Authors: Haoyun Yin
Yixuan Qiu
Xiao Wang
Venue: Accepted by TMLR ()

Abstract: Coreset selection, a technique for compressing large datasets while preserving performance, is crucial for modern machine learning. This paper presents a novel method for generating high-quality Wasserstein coresets using the Sinkhorn loss, a powerful tool with computational advantages. However, existing approaches suffer from numerical instability in Sinkhorn's algorithm. We address this by proposing stable algorithms for the computation and differentiation of the Sinkhorn optimization problem, including an analytical formula for the derivative of the Sinkhorn loss and a rigorous stability analysis of our method. Extensive experiments demonstrate that our approach significantly outperforms existing methods in terms of sample selection quality, computational efficiency, and achieving a smaller Wasserstein distance.

Published in Transactions on Machine Learning Research (02/2025)
Wasserstein Coreset via Sinkhorn Loss
Haoyun Yin Department of Statistics
yin Spurdue. edu
Purdue University
Yixuan Qiu
Sunyinuan@sue..duu.com
School Statistics and Data Science Shanghai University of Finance and Economics
Xiao Wang
wangziaa Swtrruue edu
Department of Statistics
Purdue Universily
Reviewed OpenReview: https: oppnrevvepnn..v.eu.. net/ formm? id= DrKCDS881L
Abstract
Coreset selection, technique for compressing large datasets while preserving performance, is crucial for modern machine learning. This paper presents a novel method for generating high-quality Wasserstein coresets using the Sinkhorn loss, a powerful tool with computational advantages. However, existing approaches suffer from numerical instability in Sinkhorn's algorithm. We address this by proposing stable algorithms for the computation and differentiation of the Sinkhorn optimization problem, including an analytical formula for the derivative of the Sinkhorn loss and rigorous stability analysis of our method. Extensive experiments demonstrate that our approach signilicantly outperforms existing methods in terms of sample selection quality, computational efficiency, and achieving smaller Wasserstein distance.
1 Introduction
In the last decade, big data stored in databases have grown massively and become difficult capture, form. store, manage, share, analyze, and visualize using typical database software tools (Sagiroglu & Sinanc 2013 To address this challenge, the coreset technique, also known reprosentative sampling, has been developed to select a handful of representative samples that summarize the original dataset. Formalizing the notion of representative requires care, as representative sample for a clustering algorithm may differ from that for a classification algorithm (Claici et al. 2018). Consequently, coreset construction methods usually depend on specific tasks, such as clustering (Feldman 2020 , classification (Coleman et al. 2019). logistic regression (Wang et al. 2018), mixture models (Lucic al.) 2018), low-rank approximation (Cohen et al. 2016), matrix
sketching Drineas et al. 2012), and Bayesian inference (Campbell & Broderick 2018).
In this context, it is essential to establish clear framework for discussing coreset problems. Let B denote the number of samples in a coreset, and define the set Ag = {VB VB = 1 peeeserpoobbee & De X, X), which consists of all probability measures supported on B points, each having mass B. Given a measure & on the metric space (X'd) and a particular objective function f f X - IR, we say that un E An a d-coreset of g for the specific task defined by f if
[& [x(]( - Eug [S(X)] <d, with B N

---

Published in Transactions on Machine Learning Research (02/2025)
To extend this notion beyond single task, we consider family F of objective functions f: - R. We define the integral probability metric (TPM) associated with F as
dx(E,v) sup feF
- -E,
(1)
A measure uB E AB then called a lask-agnostic d-coreset of E dy(S, un) < 6, thereby controlling the approximation error uniformly over all functions in J.. This task-tgn formulation emphasizes the coreset's ability to perform well for wide class of objective functions, rather than being tailored to single specific task.
In particular, consider the following two classes of functions:
= = (f))))))))(1fff)))) f(u) - d(x, Vz,u e 2
and
F2 = (s 115l (l $)} where 2 is the Sobolev space ( 1:; dx.l E L2} endowned with the norm
I - (// [[x)] du(a) + J [VS((M))
For these function classes, with W1 and W2 being the 1- W-sseresttein distance and 2-Wasserstein distance, respectively, one can show that:
(a) dr, MB) = (1, (Santabbrgiio 2015
(b) d MB) < (C, HB), where C constant Peyréé et al. 2019
These results indicate that the Wasserstein distance provides a suitable upper bound for the relevant integral probability metrics. Consequently, using a Wasserstein- objective in coreset construction ensures that the resulting coreset is task-agnostic. By minimizing the Wasserstein distance, we inherently minimize an upper bound on the IPM for broad range of objective functions f Fi or F2, and approximates the target measure with respect to a large family of tasks.
These insights into lask-agnostic coresets rely on integral probability metrics and our specific focus on Wasserstein- based bounds. However, task-agnostic coreset constructions have also been explored through other frameworks. few examples of task-agnostic corescts are the mso-rep-points (Fang & Wang 1993)
energy-rep-points MMak & Joseph 2018), and mmd-ree-points (Dwivedi Mackey 2021
The mse-rep-points, or mean squared error representative points, aim to minimize the expected distance from random point drawn from the distribution to its closest representative point. These points, also known principal points, have been effectively used in various applications, including quantizer design and optimal stratified sampling. In practice, they are often generated by performing k-means clustering on large batch sample from the distribution, and then using the converged cluster centers as representative points. However, a primary drawback of mss-rep-pore preptpointste that they do not necessarily converge to the underlying distribution, which limits their effectiveness in some scenarios (MacQueen et al.
1967 Graf & Luschgy 2007).
The second class of coresets, energy-rep-points aim to minimize a measure of statistical potential. This class includes methods such as minimum-energy designs (Joseph et al. 2015) and minimum Riesz energy points (Borodachov et al. 2014). Although these point sets converge in distribution, their convergence rate is slow, and the construction of such point sets can become computationally expensive as the dinensionality of the data increases. Mak & Joseph| (2018) improved these methods and try to generate support points by minimizing the energy distance. However, energy distance, regarded as a kernel discrepancy for nonuniform distributions with the specific kernel choice of the negative Fuclidcan norm, may not be well-suited for high- dimensional manifold data. Such data is commonly encountered in contemporary machine learning tasks, where the properties of the negative Euclidean norm can pose challenges.
2

---

Published in Transactions on Machine Learning Research (02/2025)
A third class of task-agnostic coresets, the mmd-rep--poniee aims to minimize the maximum mean discrepancy (MMD) between the empirical distribution and the distribution of the representative points. The MMD is a kernel-based metric that quantifies discrepancies between distributions, and mmd-repp-ponts are constructed by selecting points that reduce this measure. However, the effectiveness of mmd-rep-points is highly sensitive to the choice of kernel, which can significantly influence the resulting approximation. Moreover, the method proposed by Dwivedi
Mackey (2021) requires min,m, p)) memory storage and does not support mini-batch operations, making it challenging to apply to large-scale datasets frequently encountered in contemporary machine learning settings.
In contrast, the Wasserstein distance, the core of optimal transport (OT, Villani (2009)); is recognized as the natural geometry for probability measures for its efficacy, and has gained widespread use across multiple machine learning domains, including generative modeling (Arjovsky et al., 2017) and domain adaptation Courtyy et al. 2017). Let (X(, be a metric space that Polish. For p E 10),), the p- p-Wasserstein 1/p distance is defined as () (inlzen(aw) Sexx d(a, (dx,, dy)) joint distributions with marginals μ and v. According to Theorem 1.1 in Kloeckner (0012) we know that if X' = Rd, and let v denote the minimizer of W (5,))))) in the space of all measures supported on at most B points in X, i.e.,
where [(,, v) is the set of all
VB= =
argmin
Wp(S,
B E vw
Σ :-1 }
where [wi) is any set of weights, X; E X and 6x indicates the Dirac measure al X, then it holds that W Vp) = (( 1/4).
Despite the various appcaling thcoretical properties, one major barrier for the wide applications of the Wasserstein distance the diffieulty in computation, especially when the size of the dataset is large. For two discrete distributions, OT solves a linear programming problem of nm variables, where and m are the number of Diracs that define the two distributions. Assuming n m, standard linear programming solvers for OT have complexity of On lon)logn) (Pele & Werman 2009), which quickly becomes formidable as n increases, except for some special cases Peeré et al. 2019).
To resolve this issue, many approximate solutions to the Wasserstein distance have been proposed, among which the Sinkhorn loss has gained massive popularity {Cuturi 2013). The Sinkhorn loss can be viewed as rngulpic-regullrizedd Wasserstein distance, which adds smooth penalty term to the original objective function of OT. The Sinkhorn loss is attractive its optimization problem can be efficiently solved, least in exact arithmelics, via Sinkhorn's algorithm (Sinkhorn 1964 [inkhorn & Knopp, 1967), which merely involves matrix-vector multiplications and some minor operations. Therefore, it is especially suited to modern computing hardware, such as GPUs. Recent thcoretical advancements indicate that Sinkhorn's algoritht achieves an e-approximation to the unregularized OT problem with a computational complexity of On2 lo (Dvurechensky et al. 2018). Meanwhile, the accelerated Sinkhorn algorithm improves this complexity to O(n779 173 Lin al. 2022)..
However, one critical pain point of the Sinkhorn loss, though typically ignored in theoretical studies, is that Sinkhorn's algorithm is numerically unstable (Poyré et al., 2019). We show in numerical experiments that even for very simple settings, Sinkhorn's algorithm can quickly lose precision. Various stabilized versions of Sinkhorn's algorithm, though showing better stability, still suffer from slow convergence in these cases. Moreover, since modern deep generative models mostly rely on the gradient-based learning framework, it crucial to use the Sinkhorn loss with differentiation support. One simple and natural method for differentiating the Sinkhorn loss is to unroll Sinkhorn's algorithm, adding every iteration the auto- differentiation computing graph (Genevay et al. 2018 Cuturi et al.) 2019). However, this approach is typically costly when the number of iterations are large. The slow convergence rate resulted from numerical instability of Sinkhorn's algorithm may further exacerbate this issue.
In this paper, we solve the coreset problem using the Sinkhorn loss. We rigorously analyze both the computation and differentiation of the Sinkhorn optimization problem, designing corresponding algorithms that are provably efficient and stable. Utilizing the analytic differentiation of Sinkhorn loss, we have designed based algorithm for the strategic selection of the coreset, which is task-agnostic and suitable for
3

---

Published in Transactions on Machine Learning Research (02/2025)
high-dimensional manifold data. As a result, it can be applied to wide range of machine learning tasks. We have demonstrated the superior performance of our method in extensive simulation studies and practical applications with image data. Our major contributions are:
(i) We derive an analytic expression for the derivative of the Sinkhorn loss, which can be efficiently computed during the back- propagatin phase in machine learning.
(ii) We have rigorously analyzed the advocated algorithms for Sinkhorn loss computation and differentiation (SLCD), and show that they have desirable efficiency and stability properties.
(iii We designed task-agnostic coreset selection method called Wasserstein coreset via Sinkhorn loss (WCSL) and highlighted its superior performance in low-budget scenarios, in extensive simulation studies and practical applications with imaging data.
This paper is organized as follows. In Section 2 we introduce the notation and provide the necessary background on the Sinkhorn loss. The construction of a Wasserstein coreset can be viewed as a two-layer optimization problem: the outer layer selects the coresct points, by minimizing the Sinkhorn loss. The inner layer computes the Sinkhorn loss, which requires optimization of the transport plan. In Section we present the outer layer, WCSL algorithm, and explain why computing the derivative of the Sinkhorn loss is essential. Section g addresses numerical challenges in the Sinkhorn algorithm and introduces the SLCD algorithm to efficiently compute and differentiate the Sinkhorn loss, resulting in the inner layer. Section 5 provides theoretical analysis on the convergence and stability of the SLCD algorithm, the consistency of the Wasserstein coreset, and the convergence properties of the WCSL algorithm. Finally, Section & presents numerical experiments that demonstrate the performance of WCSL in comparison to several alternative methods.
2 Background
Throughout this article we focus on discrete OT problems. We denote the (n - 1)ddmmension probability simplex by A {w E PT w-1,, 1}, and let μ Σi a;os; and v Σラ by6y, be two discrete probability measures supported on data points {r:] and {リ, b respectively, where a = (a1,-.., a,) e {T E R T1m = T M(<x)]_ (u, ( E Raxm be cost matrix with entries Mij- = d(z;, Mj) for any selected distance d(.,)) Without loss of generality, we assume that 2 m, as their roles can be exchanged. Then we can characterize the OT problem by the following optimization problem,
A and b (b1....,bbm, BT) E A.. Define I[(ab)
b}, and let
W(,p) = min (P,M), Pell(a.))
(2)
where (A, B).
For the Sinkhorn loss, we use the following notation. For x. E RR, x A means min v). For vector = "t....,,u_1)) and use diag(t) to denote the diagonal matrix formed by v. Let u= (ui,..., w)" be another vector, and use u @ o to denote the | x k matrix with entries (u, + my). For a matrix A= = (ar) ) = (A AA) with column vectors , let A= = (A and [x[ be the matrix with entries ed The symbol @ denotes the clementwise multiplication operator between matrices or vectors. || and || || stand for the Euclidean norm for vectors and Frobenius norm for matrices, respectively.
(U,...., 2T)TT let v -- (o ...., e " ê
An optimal solution to [2) denoted as P", is typically called an optimal transport plan, and can be viewed as a joint distribution whose marginals coincide with and v. The optimal value W W(u,v) = (P*,)) /) is then called the Wasserstein distance between and v if the cost matrix M satisfies some suitable conditions (Proposition 2.2 of Peyré et al. 2019).
Solving the optimization problem can be difficult even for moderate and m. One approach regularizing the optimization problem to add entropic penalty term to the objective function, leading
4

---

as an abbreviation for Published in Transactions on Machine Learning Research (02/2025)
to the entropie-regularrizd OT problem (Cuturi 2013):
(x = min Tell(a,b)
- (T,M)-X (TT)
(3)
where h(T) Ei=1 -i=1 [1 loggooi Tiş) is the entropy term. The new objective function is - -strongly convex on ll(a, 6), (3) has a unique global solution, denoted Tx i.e.
T argmin (T, h(T) TEll(a.b)
In this article, Tx is referred to as the Sinkhorn transport plan. To simplify the notation, we omit the subscript in 7x hereafter when no confusion is caused. The entropic-regularized Wasserstein distance, also known as the Sinkhorn distance or Sinkhorn loss in the literature (Cuturi 2013), is then defined as
Vx = (T",
(4)
Tt worth noting that in the literature, Sx and 5x are sometimes referred to as the sharp and regularized Sinkhorn loss, respectively. this article we focus on the sharp version Sx, and simply call it the Sinkhorn loss for brovity, as S^ achieves a faster rate at approximating the Wasscrstein distance than Sx. suggested by the following proposition.
Proposition 2.1 (Luise al. 2018). There cxist conslants C₁, C2 > O such that for any X 0, (S, ((, S C and [Sx(,,0) ((,,v)) s C2/2. The constants C₁ and C2 independent of X, and depend μ and v.
3 Wasserstein Coreset via Sinkhorn Loss (WCSL)
In this section, we will use the sinkhorn loss to generate the Wasserstein coresct. We call this method the Wasserstein coreset via Sinkhorn loss. The objective of WCSL is to find the empirical distribution of the Wasserstein coreset:
L/B argmin VB) mgEAg
(5)
As we have mentioned earlier, the standard solver of Wasserstein distance quickly becomes formidable as n increases. To address this issue, we propose to use the Sinkhorn loss to approximate the Wasserstein distance in the coreset problem. The objective of the empirical distribution of the Wasserstein coreset via Sinkhorn loss is revised as:
1/B argmin Sx vB).
(6)
"BEAg
The algorithm for WCSL (Algorithm is an iterative procedure that seeks to minimize the Sinkhorn loss between the full measure and the measure uniformly supported on the coreset.
Given a mcasure &, the basic idea the algorithm is randomly sample an initial set D(0) consisting of B i.idd sample points from ç, and then update the points we get with gradient descent based on the gradient of the Sinkhorn loss between the measure uniformly supported our selected points and the measure &, until convergence. Then output the optimized points D* as the coreset, given by
1 B
<
D argmin SA {x
&,
6x
With a slight abuse of notation, we denote Sx (x((x) Br {y By
B2
5x h EB, &x., 1 E B. E where { and B2 two sets of points.
In detail of WCSL, we use one mini-batch of i.id sample yll,t) from & update our coreset in each iteration. The preparation step within each iteration is to compute the matrix MDDI,.) y(.)), which ((D.))),,)..)))
5

---

Published in Transactions on Machine Learning Research (02/2025)
Algorithm Algorithm for Wasserstein corcset via Sinkhorn loss (WCSL)
iid. 1: Sample D = [x; 1 M=1 &. Set 1 = 0.
2: repeat
Initialize D(1.0) t- D())
for - 0, ...,, T do (1.1) 6: Compute M MD((l.1, y((tt) the distance matrix based on the selected distance.
5: Resample y1(.1) = {xm mm=1 ~ 5
7: Update D(I,++1) c- D(l.) + Du..MVi (D(l..)) (Or Adam update)
8: end for
Update D(++1) &- DDU.T)
10: Set. 1 4 1+11 1.
11: until (x(D()) ((D()) D() <8
return DII)
12:
holds the pairwise distances, typically the squared L2 norm, between the elements of D(I.) and y(h,)) The squared L2 norm is chosen to ensure compatibility with the definition of the 2- Wassersst distance (W2). This distance matrix M is then utilized to update D(G.) using the chain rule, as V Sx ((())) V (D(l.)) V computation of V Sx ( poscs a significant challenge, and this will be discussed in detail in Section where the Algorithm for the compulation and differentiation of the Sinkhorn loss is presented.
In general, we can get Da.. M by automatic differentiation. However, the
4 Computation and Differentiation of the Sinkhorn Loss
To use the Sinkhorn loss in deep neural networks or other machine learning tasks, it is crucial to obtain the derivative of Sx(u, with respect to the distance matrix. Differentiating the Sinkhorn loss typically involves two stages: the solution stage and the differentiation stage. In the solution stage, the Sinkhorn loss the transport plan is computed using some optimization algorithms, and in the differentiation stage, the derivative is computed using either analytic expression analyze both stages in detail, we first present the issues of Sinkhorn's algorithm.
automatic differentiation technique. Before
4.1 Issues of Sinkhorn's 's Algorithm
In the existing literature, one commonly used method for the computation of the Sinkhorn loss is Sinkhorn's algorithm (Sinkhorn, 1964 Sinkhorn & Knopp, 1967). Unlike the original linear programming problem ( the solution to the Sinkhorn problem has special structure, Cuturi (013 shows that the optimal solution T can be expressed as
T = diag(u u " MM diag(u")
(8)
for some vectors x² and v", where M, = (e-AMy). Sinkhorn's algorithm starts from an initial vector ((0) e R., and generates iterates (k) E RT and (1)) E RT as follows:
(++1) t a 2 [Mau
(1+1) + bo [M. (1+1)-1)
(9)
It can be proved that (*) - n² and (1)) + v", and then the Sinkhorn transport plan T can be recovered
by S
Sinkhorn's algorithm is very efficient, as it only involves netrix-vector multiplication and other inexpensive operations. However, one major issue of Sinkhorn's algorithm that the entries of M = (e-AMs) may easily underflow when À is large, making some elements of the vectors Mk)(k) and MT (++11) close to zero. As a result, some components of (k++)) and (((+1)) would overflow. Therefore, Sinkhorn's algorithm in its original form is unstable, and in practice the iterations which we call the Sinkhorn-log algorithm for simplicity. In addition, there are some other works that also attempt to improve the numerical stability of Sinkhorn's algorithm (Schmitzer 2019 Cuturi et al. 2022
9 are typically carried out in the logarithmic scale,

---

Published in Transactions on Machine Learning Research (02/2025)
50
15
Truth
Sinkhom
Sinkhom-log
Stabilized
Greenkhomn
L-BFGS
Figure 1: Visualization of Sinkhorn plans computed by different algorithms.
Remark: Due to numerical instability and slow comvergence, the transport plan obtained through Sinkhorn's algorithm and verious widely used algoritoms, deviates significantly from the ground truth. In contrast, the transport plan computed using L-BFGS demonstrates a near-perfect visual abigument with the true solution. Detailed experimental settings are provided in Appendir B.1
Algorithm Sinkhorn Loss Computation and Differentiation (SLCD)
Solution Stage:
Use L-BFGS to solve 3 argming f(3), where f(3) is defined in 12 and its gradient Vg in
13
Compute a" using 11 and compute 7 using 14
4:
Compute the Sinkhorn loss using
5:
Differentiation Stage:
Compute the analytical derivative using 15 from Theorem 4.1
Despite the advancements of Sinkhorn's algorithm, one critical issue observed in practice is that Sinkhorn- type algorithms may be slow to converge, especially for small regularization parameters. This would severely slow down the computation, and may even give misleading results when the user sets a moderate limit on the total number of itcrations.
In Figure 1 show a motivating example to highlight these issues. Consider two probability measures for the Sinkhorn problem, and the true 7" matrix is visualized in Figure ! along with the solutions given by various widely-used algorithms, including Sinkhorn's algorithm, Sinkhorn-log, the stabilized scaling algorithm (Stabilized, A Algorithmm 2 of [cchniteerr 2019), and the Greenkhorn algorithm (Altschuler et al. 2017 Lin et al. 2022). In Figure [ it is clear that the plans given by Sinkhorn's algorithm and Greenkhorn are farthest to the true value, and Greenkhorn generates NaN values reflected by the white stripes the plot. In contrast, the stable algorithms Sinkhorn-log and Stabilized greatly improve them. Sinkhorn's algorithm and Sinkhorn-log equivalent in exact arithmetics, their numerical differences highlight the need for numerically stable algorithms, However, Sinkhorn-log and Stabilized still have visible inconsistencies with the truth even after 1000 iterations.
4.2 The Advocated Alternative for Sinkhorn Loss Computation
To solve the numerical instability and slow convergence issue of Sinkhorn's algorithm, we advocate an alternative scheme to solve the optimal plan T*, as shown in the solution stage of Algorithm 2 Then we show both theoretically and empirically that this method enjoys great efficiency and stability.
Consider the dual problem of [3] which has the following form (Proposition 4.4 of Peyré al. 2019
max C(a, 8) max a 8Tb X a
Q Rd 8 mm
(10)
a

---

Published in Transactions on Machine Learning Research (02/2025)
Let a* = (ai....,,ç)) o and 3 (3i..... Pm)) be one optimal solution to (10) and then the Sinkhorn transport plan 7 can be recovered as 7 = ex ex(a* @ 8 M.. Remarkably, 110) is equivalent to an unconstrained convex optimization problem, so a simple gradient ascent method suffices to find its optimal solution. But in practice, quasi-Newton methods such as the limited memory BFGS mothod (L-BFGS, Liu
& Nocedal 1989) can significantly accelerate the convergence.
It is worth noting that we can reduce the number of variables to be optimized in (1O) based on the following two findings. First, as pointed out by Cuturi et al. (2019), the variables (a, B) have one redundant degree of freedom: if (a*,8*) is one solution to (10), then is (a" cln, 8 clm) for any c. Therefore, we globally set Bm 0 without loss of generality. Second, let a" (3) arg max, L(a, B) for a given 3 = ( _
35=,0²= and define ((3) - (3), Then we only need to minimize f(3) with (m - 1) free variables to get 3r, and a can be recovered as a = (3*).. As a result, we have the following lemma, with proof given A Appendix C.2
Lemma 1. The dual problem 10 is equivalent to ming f(3), where a² simple closed-form expressions:
a (8), f(3), and Vgf have the following
m 100s Mss)
a'(8), X log a; x
(11)
S(3) -a (3)) a- 3dh + +
(12)
Vaf = 1, b.
(13)
Then for specific B, we can recover the transport plan T(3) by:
T(3) ex(a (aala"(()) (3) o M).
(14)
With f(B) and Vaf. the L-BFGS algorithm be readily used to minimize f(3) and obtain 8r. Each gradient evaluation requires O(mn) exponentiation operations, which is comparable to Sinkhorn-log. Although exponentiation is more expensive than matrix-vector multiplication as in Sinkhorn's 'iss algorithm, this extra cost can be greatly remedied by modern hardware such as GPUs. On the other hand, we would show in Section that L-BFGS has a strong guarantee on numerical stability, which is critical for many scientific computing problems.
Using L-BFGS to solve (10) is known practice (Cuturi & Peyré 2018 Flamary al., 2021), but little is
known about its stability in solving the regularized OT problem. For the motivating example in Section
4.1 we demonstrate the advantage of L-BFGS by showing its transport plan in the rightmost plot of Figure We limit its maximum number of gradient evaluations to 1000, and hence comparable to other methods. Clearly, the I-BFGS solution is visually identical to the ground truth.
To study the difference between L-BFGS and Sinkhorn's 's algorithm more depth, we compute the objective function value of the dual problem [10 at each iteration for both Sinkhorn-log and L-BFGS. The results are visualized in Figure 2 with three different values of --1, X-1 0.1, 0.01, 0.001. Figure 2 gives clear clue to the issue of Sinkhorn-log: it has surprisingly slow convergence speed compared to L-BFGS when A-1 is small. Thcoretically, Sinkhorn algorithms will eventually converge with sufficient iterations, but in practice, a moderate limit on computational budget may prevent them from generating accurate results. To this end, L-BFGS appears to be a better candidate when one needs a small X for better approximation to the OT problem. In Appendix B.2 We have designed more experiments to study the stability and accuracy of the Sinkhorn loss computation across different algorithms.
4.3 The Analytic Differentiation of Sinkhorn Loss
For the differentiation stage, one commonly-usee method is unrolled Sinkhorn's algorithm, which is based on the fact that Sinkhorn's computation algorithm differentiable with respect a, b, and M. Therefore, one can use automatic differentiation software compute the corresponding derivatives in the differentiation stage. This method is used in Genevay et al. (018)) for learning generative models with the Sinkhorn loss,
8

---

Published in Transactions on Machine Learning Research (02/2025)
Objective Function Value
-
1-1- 0.01
-01000001
0.5
-05
0.0
-05
-0.5
-10
-1.0
Sinkhorn-log
-1.0
-15
1.5
L-BFGS Gradient Eval.
L-BFGS Iteration
-20
-.0
-25
2.5
-2.5
-3.0
-30
100
200
1000
Iteration Count
Figure 2: Comparing the convergence speed of Sinkhorn-log and L-BFGS.
Remark: Sinkhorn-log converges slowly compared I-BFGS, especially when the regularization parameter X-1 of the Sinkhorn loss is small. Since each L-BFGS iteration may ineolve than gradient evaluation, for L-BFGS plot the values against both the outer iteration and gradient evaluation counts.
but in practice it may be extremely slow if the solution stage takes large number of itcrations. To avoid the excessive cost of unrolled algorithms,, various implicit differentiation methods have been developed for the Sinkhorn loss (Feydy al. 2019: Campbell et al. 2020 Xie al.) 2020 Eisenberger al. 2022), but they still do not provide the most straightforward way to compute the gradient.
To this end, we advocate the analytical differentiation of the Sinkhorn loss, which uses the optimal dual variables (a) from the solution stage to directly compute the derivative, as presented in the differentiation stage in Algorithm| 2 The analytic form for v., (u, v) has been studied in Luise et al. (2018), and to our best knowledge, few result has been presented for Sx v). Our first main theorem, given in Theorem 4.1 fills this gap and derives the analytic form for Sx,,
Theorem 4.1. For a fixed X > O,
SMSS(H, A(8a ① 8, M) 9 T
(15)
r. Mr (MOT""I mo= - OM 8, D- je (a @ pr),, and D diag(b) - diag IT", In addition, D is positive definite, and hence invertàble, for any À > 0, a A", be Am, and M.
where T ex(a" 33 s, - (vr -T*s,,, 8. = GT 0)T
Assuming n 2 m, the main computational cost is in forming D, which requires O(m?n) operations for matrix-matrix multiplication, and in computing & which costs O(m²) operations for solving a positive definite linear system. Theorem 4. provides simple and efficient way to compute the gradient of the Sinkhorn loss with respect to the cost matrix M.
5 Stability Analysis and Convergence
With the computation and differentiation algorithms advocated for the Sinkhorn loss in Algorithm are able to find the Wasserstain coreset using Algorithm [ efffiieentl. In this section, focus on the theoretical properties of these algorithms, and show that they enjoy provable efficiency and stability.
5.1 Stability Analysis of Algorithm 2
As first step, consider the objective function 12 and show that f(ß) has well-behaved minimizer, which does not underflow or overflow.
Thcorem 5.1. Let f* denote the ménimum value of f(3) and 5 an optimal solution, and let a* a² (5*). Then f -00, 8 unique, ||a*1| oc, and || co. In porticular, let I arg max; Tim: Qmax max; di, and log(n/bm)- Then Lo, S TTm. of & Ua, and Lg, S B3 < U3, < a, for all 1,....,
9

---

Published in Transactions on Machine Learning Research (02/2025)
and j 1, 1..... m, where
Ua, Min + x- log(a; A bmn);
Us, M15- p0- [oog(a A bb)) c),
La= =x- .log(aag//)) mm) max
Lg = 1-1 .looog(b///)) max (Ua,- -Mij), L.. -x- ..Mooo((a//m) max (U 8,- -Mi))-
)
U3:- [ma -MiM A-1 [log(Amax A b1) + c],
Theorem 5. shows that the optimal dual variables 7" are bounded, and more importantly, the bounds are roughly at the scale of the cost matrix entries My and the log-weights log(a;) and log(b;). Therefore, at least the larget of the optimization problem is well-behaved. Moreover, one useful application of Theorem 5.1 is to impose a box constraint on the variables, adding further stability to the optimization algorithm.
After verifying the properties of the optimal solution, more interesting and critical problem is to seek stable algorithm to approach the optimal solution. Indeed, in Theorem 5.2 we prove that the solution stage of Algorithm| 2 is one such method.
Theorem 5.2. Let { be a sequence of iterates generated by the L-BFGS algorithm starting from a fized initial value 310), and define (ih) = (B(K),, O)T, T(k) exfa" (3(8) @ g(k) M), and f(k) = f(3(k)). Then
there exist constants 0 S r < 1 and C1, C2 > O, such that for each k 0, let e(k) = - f*)::
(a) f(k) - f S e(k)
(Linear convergence for the objective value)
(6) [10(k) - 12 s CGE(1)
(Linear convergence for the iterates) (Exponential decay of the gradient)
(c) T(1) 1,m -a, IV
1, - -b112 - (2e(1).
(d) T < <minfa, bj + (C }. is m - 1.
(T(() does not overflow)
(e) max; T >a/m m, max, T 1/11)(b))))) -26 /n, 1 15 S m m.
(I(() does not underflow)
The erplicit erpressions for the constants C1, C2.r are given in Appendir
A
Theorem 5.2 reveals some important information. First, both the objective function value and the iterates have linear convergence speed, so the solution stage using L-BFGS takes ((log(1/E)) iterations to obtain an s-optimal solution. Second, the marginal error for μ4, measured by ||TT(* 1m - all, is exactly zero due to the partial optimization a. The other marginal error ||7(k)T1,, 1,, bll, which is equal to the gradient norm, is also bounded at any ileration, and decays exponentially fast to zero. This validates the numerical stability of the L-BFGS algorithm. Third, the estimated transport plan at any iteration k. 7<(), is also bounded and stable. This result can be compared with the formulation in S : it is not hard to find that u", v", and M when computed individually, can all be unstable due to the exponentiation operations, especially when X-1 is small. In contrast. T and T(k),, thanks to the results Theorem 5.2 do not suffer from this issue.
We emphasize that Theorem 5.2 provides novel results that are not direct consequences of the L-BFGS convergence theorem given in Liu & Nocedal (1989). First, classical theorems only guarantee the convergence of objective function values and iterates as in (a) and (b), whereas we provide richer information such as the marginal errors and transport plans specific to OT problems. More importantly, our results are all nonasymptotic with computable constants. To achieve this, we carefully analyze the eigenvalue structure of the dual Hessian matrix, which is of interest by itself.
Likewise, we show that the derivative of V Sx v) as in Theorem 4.1 can also be computed in a numerically stable way. Let S be the k-step approximation to VMS == SM v/ using the L-BFGS algorithm, i.., replacing every T in S by T(k) Then we show that the crror on gradient also decays exponentially fast.
Theorem 5.3. Using the symboos defined in Theorems 4-1 and 5.
let o 1/a (D), where a (D) is
the smallest eigenvalue of D. Assume that for some ko,
min 1, (6a D FF) -11 117
(ka) <
4X
10

---

Published in Transactions on Machine Learning Research (02/2025)
and then for every k 2 ko; - S Csve(k) Cs Csv p where the explicil expression for Cg is given in Appendix A ko always evists e(k) decays to zero exponentially fast, which is ensured by Theorem 5.2 a)
5.2 Consistency of Wasserstein Coreset
We examine the distributional convergence of the empirical mcasure induced by the Wasscrstein coresct to the target distribution
Theorem 5.4. Let X~ ~ ,e, XB ~ uB, where pB is the empirical distribution of Wasserstein coresets, which is defined in 50. Then as B 00, XB 4 X.
To summarize, the theorem shows that the Wasserstein coreset is indeed representative of the target distribution, when the number of points B grows large. Consoquently, this result implics the consistency of the Wasserstein coresel, critical property for any coreset algorithm.
5 (a) If X R is continuous, then ((B)) 4 g(X). (b) Ifg g is continuous and bounded, then EEgg( X B)] + [x)]
Corollary 5.5. Let X ,, Xn μn, with un defined in
Part (a) Corollary 5.5 follows from the contimuous mapping theorem and Theorem 5.4 and part (b) holds due to the Portmanteau theorem. The corollary establishes the consistency of the Wasserstein coreset for integration, thus validating its use in range of applications. Specifically, part (a) indicates that the Wasserstein coreset can be effectively employed for uncertainty propagation in stochastic simulations. Part (b) further demonstrates that any continuous and bounded function can be accurately estimated using the Wasserstein coreset, the sample average of g(xB) converges to the expectation of g(X).
5.3 WCSL Algorithm Convergence Analysis
Without assuming convexity, as convexity may not hold in practical scenarios, our aim is to demonstrate the convergence of WCSL to a stationary point and establish a robust convergence criterion. First, we demonstrate that (E, converges to stationary point in Theorem 5.6 We only provide the convergence guarantee when using the Adam update, as Adam a better performed optimization algorithm practice. We then provide a convergence criterion for the WCSL algorithm Proposition 5.7
Theorem 5.6. Using the Adam update, after N iterations for some N N such that N
ft T- where
B1 is the first moment decay rate in Adam, the Sinkhorn loss Sx D(M)) converges to stationary point. Specifically, there crists a constant C' > such that
< log N VN
E DA
Theorem 5.6 establishes the thcoretical convergence property of WCSL after a finite number of iterations. Now we need a practical method for monitoring the convergence of the WCSL algorithm and ensuring that the algorithm is progressing towards a stationary point of the Sinkhorn loss.
As mentioned in the WCSL algorithm (Algorithm D the stopping criterion is (D(()), (((())((-)))) < 8. We now provide justifications for this stopping criterion to be proper convergence criterion.
Proposition 5.7. Let & be the target distribution, and (120) be a sequence of representative samples updated by Algorithm 4 If Sx (DI-11) 8, then [Sx(E, D()) SX Dil-1)) <
The convergence of the WCSL algorithm can be characterized by (S, Sx D(1-1) D < 6, for chosen 8 0. This criterion implies that the difference in the objective function, specifically the Sinkhorn loss between successive iterations, is sufficiently small, indicating that the algorithm is approaching a stationary point. Given Proposition 5.7 we know that Sx (D()),, D(l-1) (D())) serves as an upper bound for [Sx D()) SA, D((-))) Thus, we can ascertain that the algorithm is progressing towards a stationary point of (X(E, if (D("), (,(D)), which can be computed efficiently, small. This provides a practical method for monitoring the convergence of the WCSL algorithm and ensuring that the iterations yield diminishing improvements, thus verifying the convergence of the algorithm.
11

---

Published in Transactions on Machine Learning Research (02/202))
Sandom
Semara
Normal
Normal Thimn
VIC
-2
-2
Spiril Rancom
Spiral Kimeans
Spiral
Sahall Semal Thinring
Figure 3: Comparative visualization of samples selected from 2-dimensional datasets.
6
Numerical Experiments
6.1
Simulations
Several simulations are presented to demonstrate the effectiveness and efficiency of the WCSI. method. We evaluate the sample selection quality of the WCSL method visually and in terms of distance metrics. We also performed a runtime analysis for the WCSL method. The results provide insights into the practical applicability of the WCSL method in various data settings.
To evaluate the performance of the WCSL method, we compared it against several alternative sampling techniques: random sampling, the centroids obtained via k-means clustering, the SCCP method (Mak & Joseph 2018), and kernel thinning (Dwivedi & Mackey 2021). Random sampling, which selects samples uniformly at random from the dataset, serves straightforward benchmark. The k-means centroids approach involves partitioning the dataset into clusters and then using the resulting cluster centroids representative points. The SCCP method compresses continuous probability distribution into a finite set of support points that minimize the energy distance. This approach guarantees that the support points converge in distribution to the original distribution and often yiclds improved integration error rates for a broad class of functions compared to standard Monte Carlo methods. Finally, the kernel thinning method selects representative points by minimizing the maximum mean discrepancy (MMD), thereby providing another kernel-based strategy for dataset compression.
To visually show the sample selection quality difference between random samples and coreset techiniques, we present a comparative visualization of subsamples selected from synthetic 2-dimensional datascts. The experiment settings include datasets of n = 1000 points each, reduced to coresets of B = 100 points. The three configurations tested are a normal distribution, a 2-dimensional spiral, and a circular distribution. Figure 3 illustrates the results of coresct techniques. All the coresct techniques demonstrate superior performance in capturing the original data distributions compared to random sampling. Random sampling is not consistent and often fails to accurately represent. the underlying distríbution, highlighting the necessity
12

---

Published in Transactions on Machine Learning Research (02/2025)
Waszeratein Dretance by Coreest Size Distnbbtton
MMD by Coreset Size Dietribution
2.25
2.00
2.00
1.75
1.90
Nomalzze
1.50
Nornalze
1.00
1.00
0.75
0.90
0.50
150
200
200
400
Wasserstein
Careset
by Size Distributior
1.0
Nomelladd Wisseesstenn Chetence
MNO
0.8
0.0
0.8
0.5
Coreant Size
Coreeet Size
Figure 4: Comparison of different sample selection techniques across various datasets. Lower is better.
of using coreset techniques. Among the methods tested, the WCSL method performed well, effectively preserving the essential data characteristics the reduced sample sets. This visual assessment underscores the potential of the WCSL method in maintaining the integrity of the original data structure in various data reduction scenarios.
However, visual comparisons alone insufficient to rigorously determine which coresct method performs belter. Therefore, compute both the Wasserstein distance and MMD between the full dalasel and the coreset obtained using various techniques. We use 100-dimensional t-distribution and 100-dimentional normal distribution with n = 10000 points as our datasets, and reduce to corcsets of various sizes, ranging from B = 25 to B = 400. In our analysis, we seek to ensure meaningful comparisons across different sample selection techniques. Consequently, we normalize the computed Wasserstein distance and MMD values by the distance metrics between the randomly sampled subsets and the full dataset. This normalization procedure aligns the distance scales, allowing direct comparisons to be made between the various coreset construction methods.
The results in F demonstrate that the coresets generated by WCSL consistently achieve low normalized MMD and Wasserstein distance values. This indicates that WCSL more effectively preserves the original data distribution and structure. In contrast, random sampling leads to persistently high distance metrics, underscoring its limited ability to capture the underlying distribution. While the k-mcans centroids approach sometimes outperforme random sampling, it displays instability and can yield results even worse than random sampling in certain instances. The SCCP and kernel thinning methods perform moderately well, but their effectivencss varies across different metrics and scenarios. Overall, these findings suggest that WCSL provides a more robust and reliable representation of the original data, ensuring that its essential statistical properties are well-maintained.
Before evaluating the computational efficiency of the WCSL method empirically, we first present its time complexity analysis. The WCSL algorithm involves two layers of optimization: the inner SLCD problem and the outer coreset construction. To analyze the SLCD problem, we decompose it into solution stage
13

---

Published in Transactions on Machine Learning Research (02/2025)
Wasscrstoin Corosct Runtimo Analysis (p=10, B-250)
Runlime per lteration [5]
30
20
15
10
5
200,000
400,000 p0n
600,000
500 000
1,000,000
Number of Itorations
00
1,000,000
Full Datasel Sample Size n
Figure 5: Runtime analysis of the WCSL method with respect to full dataset size
and a differentiation stage. In each SLCD step, the solution stage involves computing the distance matrix, which requires (mm2p) time, where m represents the batch size and denotes the data dimension. Fach iteration of the L-BEGS optimization procedure has a complexity of O(m²), and the algorithm typically requires ((log(//e) iterations converge, where denotes the convergence tolerance. In the differentiation stage, computing the derivative incurs an additional computational cost of 0(m²). As result, the overall time complexity of cach SLCD step is O(m² (ol(pppoooo////// (0/((1c) +m)) Since < n, where is the full dataset size, treat negligible. Consequently, because each iteration processes the entire dataset in mini-batches. the time complexity of each WCSL iteration simplifies to (mp)).
Multiple iterations are typically needed for WCSL to converge, and the number of these iterations depends on the chosen stopping criterion. Based on our empirical observations, for a fixed slopping criterion, we assume that O(n) iterations are required. Under this assumption, the overall time complexity (n²pp) which remains tractable for large-scale data scenarios.
To further demonstrate the computational efficiency of the method, measure the execution times for range of dataset sizes while fixing the data dimension p 10 and coreset size B = 250. We vary n from 100,000 to 1,000,000, and record the per-iteration runtime. Figure 5 shows that the runtime per iteration increases proportionally with n, reflecting the anticipated linear scaling in computational demand. Moreover, the number of required iterations also grows roughly linearly with n, underscoring the suitability of the WCSL method for large-scale data applications. Together, these results confirm that WCSL provides an efficient and practical solution for large-scale problems.
We also provide the rumning time of the WCSL with respect to coreset size B and data dimension p in Figure G The running time is measured with the coresct size varied from B = 100 to B = 1000 given sample size It 10000 and data dimension
25. The rumning time is also measured with the data dimension varying from 10 to p = 100 given the sample size n 10000 and the coreset size = 250. The running time per iteration increases linearly with respect to the coresct size and stays almost flat with respect to the data dimension p. while the number of iterations required for convergence decreases as the coreset size B increases and increase as the data dimension p increases. These findings highlight the trade-offs bet
14

---

Published in Transactions on Machine Learning Research (02/2025)
Wasserstein Corcsct Runtime Analysis (n=10000, p=25)
Runtime per leeraion [8)
20
10
250
500
750
1000
1250
1500
1750
2000
Numbar of Iterations
130
120
110
Coresel Size B
Wassorstein Coresct Runtimo Analysis (n=10000, B=250)
2.4
2.2
Rumlime per lleratoon
20
1.8
1.6
10
200
400
600
800
1000
120
100
Number of lteraaions
80
60
40
Data Dimensionallly p
Figure 6: Runtime analysis of WCSL with respect to coreset size B and data dimension
coreset size and computational efficiency, providing valuable insights for optimizing the WCSL method for large-scale data applications.
15

---

Published in Transactions on Machine Learning Research (02/2025)
Accuracy vs Coreset Size with 95% Confidence Interval
0.9
0.8
0.7
Accuracy
0.6
0.5
Method
Random
0.3
SCCP
WCSL
50
100
150
200
250
300
350
400
Coreset Size B
Figure Classification accuracy across various coreset sizes, showcasing the performance of prediction models trained subsets selected via different methods.
6.2 Active Learning
One of the potential applications of the Wasscrstein coreset is in the field of active learning. In active learning, involves the data points X E X' and their corresponding labels y = (1,...,C), C), with C being the number of label categories. A labeling function f : X + y determines the accurate label for any given sample in active learning scenarios. Our primary goal is minimize the expected risk given data distribution with a parameterized loss function ((), This objective is formally expressed min. [E[((X, Y; w;)).. Tdeally, we have access to labeled dataset of N labeled pairs {(X,, Yi) to reduce the expected risk.
However, in practice, we have an unlabeled dataset of N observations {X.] and budget to get 9(X)), the labels of X,, for B times, where B < N. Building on the work of Sennr & Savarese (2018), we recognize that the expected risk is bounded by the summation of three factors: the generalization error from training with fully labeled dataset instead of the true population measure, the empirical risk from training on a selected subset of B points, and the discrepancy in loss between training on the entire dataset versus the selected subset, the subset loss. The mathematical representation of this constraint is:
EI w) & B( 1 N
B Σ ((x,, ((x)): j=l
N Σ ((x,,,
B
generalization
empirical risk
4x, (()):w)
(( ((X))w)
B
subset loss
The generalization error remains unaffected by the chosen subset, and when B is small, the empirical risk is typically negligibly small. This suggests that active learning strategies should prioritize minimizing the loss associated with the chosen subset.
16

---

Published in Transactions on Machine Learning Research (02/2025)
MNIST and FashionMNIST Coresets
7
2
0
O
1
6
5
2
5
9
4
3
4
8
4
6
4
7
O
9
9
4
8
7
2
0
6
/
0
4
8
7
4
9
9
5
4
4
9
4
6
O
8
/
3
0
4
6
6
2
O
/
6
3
8
O
2
7
3
/
7
7
7
9
3
8
0
5
1
3
5
1
/
/
1
9
3
7
5
9
/
6
2
0
7
7
O
6
5
3
8
7
6
4
2
8
3
6
Figure 8: Comparative visualization of coreset samples on MNIST and FashionMNIST datasets.
Notice that if E F, then the subset loss would be bounded by the coreset loss:
B
Σ (xx, w)
-EExx, - j=1
1
Q(xj) w) S d)(E,vB)
This motivates us to use coreset technique in enhancing the elficiency of choosing subsets in active learning, We consider family of 500-dimensional I-distributions, categorized into 25 classes. Fach class is defined by one-unit shift to the right in a randomly selected dimension of the t-distribution. The full dataset comprises 20,000 points. To establish a baseline, we train a logistic regression model on the entire dataset and use this model to relabel the dataset, treating these labels as the ground truth for subsequent classification tasks. Subsequently, we generate coresets using random selection, SCCP, and WCSL method. Logistic regression models are then trained on these coresets of sizes from 25 to 400, with their classification accuracies evaluated thereafter.
Given the complexity of the datasct with 25 classes, achieving high classification accuracy with small sample sizes is challenging. Especially for coreset sizes of 25 or 50, il is likely that some classes will not be represented in the coresets unless an effective coreset selection technique is used. This can lead to high classification error rates. Despite this challenge, Figure [ shows that models trained on Wasscrstein corcsct-selected samples demonstrated superior performance compared to those trained on coresets generated by random sampling or SCCP, particularly with smaller coreset sizes. This improvement classification accuracy underscores the ability of the Wasserstein coreset approach capture the underlying distribution of the full dataset effectively. Consequently, models trained using Wasserstein-solectee coresets exhibit enhanced classification performance, validating the strategic advantage of this method in active learning contexts.
6.3 MNIST & FashionMNIST
To further demonstrate the practical applicability of the WCSL method, applied it to image datasets, the MNIST and FashionMNIST datasets. Representative samples generated by the WCSL method are depicted in Figure
g demonstrating the ability of the method to retain characteristic data features in the reduced sample set.
We also evaluate the performance of the WCSL method in terms of distance metrics, comparing with other subsampling techniques. The kernel thinning method proves to be very slow on image datasets; thus, we do not include its results. The results, presented in Figure indicate that the WCSL method consistently
17

---

Published in Transactions on Machine Learning Research (02/2025)
Wesseretein Detance by Coreeat Size MNIST
MMD by Coreset Size -NIST
ACU
Noralrzdd Viaceersee
1.50
0.0
1.25
0.7
1.00
0.75
0.0
300
360
200
400
by
1.0
5.0
Norallze Cletance
0.8
2.5
0.8
2.0
1.5
100
Coreant Size
Coneeet Size
Figure 9: Comparative visualization of Wasserstein distance and MMD by coreset techniques and size. Lower is better.
exhibits lower Wasserstein distances and reduced MMD values compared to random subsampling. This affirms the method's superior capability in approximating the original data distribution, ensuring that the roduced samples are highly represontative of the full dataset.
7 Discussion
In conclusion, this paper has introduced a novel and robust Wasserstein coreset construction method that utilizes the entropic-regularized framework of the Sinkhorn loss to effectively approximate large datasets while preserving their intrinsic geometric structure. Through rigorous analysis and exlensive experimentation, we have demonstrated the stability and efficiency of our computation and differentiation algorithms in settings that have traditionally posed numerical challenges. Our findings reveal that Wasscrstein corcsets yield superior performance in model accuracy, particularly under budget constraints where label scarcity is issue in active learning scenarios. This work not only ad advances the methodology of data distillation, but also offers tangible benefits for active learning and other machine learning tasks that require fast, accurate, and cost-cffective data approximation techniques. The success the WCSL method in varions practical scenarios, including imaging data, confirms its potential as valuable tool big data analytics, affirming role in enhancing the scalability and accessibility of machine lcarning applications.
In light of these advances, it is pertinent to reflect two specific aspects of the coreset methodology within the broader context of data representation and active learning strategies. First, while the coreset is designed to mirror the underlying data distribution, it does not inherently ensure that the samples are drawn from the original datasel. This distinction becomes critical in applications where authenticity is paramount, such as dataset splitting for training and validation in machine learning models, financial audits where transactional integrity is necessary, or medical studics that require exact patient data without generalizations. Ensuring the representativeness of coresets drawn from the original dataset might involve additional constraints in the selection am post-hoc verification procedures confirm the origin of the data points. Secondly, the
18

---

Published in Transactions on Machine Learning Research (02/2025)
active learning framework benefits from the Wasserstein coreset's efficiency, particularly in the initial stages of model training. However, as learning progresses, there might be need for more explorative approach to sample selection. This could involve adaptively adjusting the coreset as new labels are acquired or integrating uncertainty measurements to guide the selection process. Exploration strategies such as entropy maximization or variance reduction could be particularly beneficial in iterative learning settings, ensuring that each new query contributes maximal information the learning model. In both cases, the challenge lies in striking balance between representativeness and practical constraints. Future work may explore hybrid strategies that leverage the Wasserstein coreset's strengths while addressing its limitations in specific application scenarios. The goal is to refine the technique to become a versatile tool that is able to handle a wide array of tasks in the burgeoning field of machine learning.
Acknowledgements
Xiao Wang's research was supported in part by the NSF Grant supported in part by National Natural Science Foundation of China (12101389), Shanghai Pujiang Program (21PJC056), MOE Project of Key Research Institute of Humanities and Social Sciences (22JJD110001), and Shanghai Research Center for Data Science and Decision Technology.
Yixuan Qiu's work was
References
Jason Altschuler, Jonathan Niles Weed, and Philippe Rigollet. Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration. In Advances in neural information processing systems 30, 2017.
Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasscrstein GAN. arxiv.org, 2017. ISSN 2331-8422.
Sergiy V Borodachov, Douglas P Hardin, and Edward B Saff. Low complexity methods for discretizing manifolds via riesz energy minimization. Foundations of Computational Mathematics, 14:1178-1208, 2014.
Richard H Byrd, Jorge Nocedal, and Ya-Xiang Yuan. Global convergence cass of quasi-Newton methods on convex problems. STAM Journal on Numerical Analysis, 2475) 1987.
Dylan Campbell, Liu Liu, and Stephen Gould. Solving the blind perspective---poitt problem end-to-end with robust differentiable geometric optimization. In European Conference on Computer Vision, pp. 244-261. Springer, 2020.
Trevor Campbell and Tamara Broderick. Bayesian coreset construction via greedy iterative geodesic ascent. arXiv.org, 2018. ISSN 2331-8422.
Sebastian Claici, Aude Genevay, and Justin Solomon. Wasserstein measure coresets. or arXiv (Cornell University), 2018. ISSN 2331-8422.
Michael B Cohen, Cameron Musco, and Christopher Musco. Input sparsity time low-rank approximation via ridge leverage score sampling. arXiv.ory, 2016. ISSN 2331-8422.
Cody Coleman, Christopher Ych, Stephen Mussmann, Baharan Mirzasolciman, Peter Bailis, Percy Liang, Jure Leskovec, and Matei Zaharia. Selection via praxy: Efficient data selection for deep learning, arXiv. org, 2019. ISSN 2331-8422.
Nicolas Courty, Remi Flamary, Devis Tuia, and Alain Rakotomamonjy. Optimal transport for domain adaptation. IEEE transactions on pattern analysis and machine intelligence, 3059):1853-1865, 2017. ISSN 0162-8828.
Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. volume 26, 2013.
Marco Cuturi and Gabriel Peyré. Semi-dual regularized optimal transport. arXiv preprint arXiv: 2018.
19

---

Published in Transactions on Machine Learning Research (02/2025)
Marco Cuturi, Olivier Teboul, and Jean-Philippe Vert. Differentiable ranking and sorting using optimal transport. volume 32, 2019.
Marco Cuturi, Laetitia
Yingtao Tian, Charlotte Bunne, Geoff Davis, and Olivier Teboul. Optimal transport tools (OTT): A JAX toolbox for all things Wasserstein. arXiv preprint arXiv:2201.12324, 2022.
Alexandre Défossez, Léon Bottou, Francis Bach, and Nicolas Usumier. A simple convergence proof of adam and adagrad. arXiv preprinl Xiv:20003 02395, 2020.
Petros Drineas, Malik Magdon-Ismail, Michael W Mahoney, and David P Woodruff. Fast approximation of matrix coherence and statistical leverage. The Journal of Machine Learning Research, 13(1):3475-3500 2012.
Pavel Dvurechensky, Alexander Gasnikov, and Alexey Kroshnin. Computational optimal transport: Complexity by accelerated gradient descent is better than by Sinkhorn's algorithm. In International conference on machine learning, pp. 1367-1376, 2018.
Raaz Dwivedi and Lester Mackey, Kernel thimning. arXiv preprint arXiv:2105.058: 05544, 2021.
Marvin Eisenberger, Aysim Toker, Laura Leal-Taixé, Florian Bernard, and Daniel Cremers. Scalable sinkhorn backpropagation. OpenReview preprint idzuR7707L555,, 2022.
Yizheng Fan. Schur complements and its applications to symmetric nonnegative and Z-matrices. Linear algebra and its applications, 353()-3):289 307, 2002.
Kai-Tai Fang and Yuan Wang. tumber-th methods in statistics, volume 51. CRC Press, 1993.
Dan Feldman. Coresets: An updated survey. Wiley interdisciplinary reviews. Data mining and knowledge discovery, 10(1):e1335 n/a, 2020. ISSN 1942-4787.
Jean Foydy, Thibault Séjourné, François-Xavee Xavier Vialard, Shun-ichi Amari, Alain Trouvé, and Gabriel Peyré. Interpolating between optimal transport and MMD using Sinkhorn divergences. In The 22nd International Conference Artificial Intelligence and Statistics, pp. 2681 221-2690, 2019.
Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar 7. Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, et al. Pot: Python optimal transport. Journal of Machine Learning Research, 22(78):1-8, 2021.
Aude Genevay, Gabriel Peyré, and Marco Cuturi. Learning generative models with sinkhorn divergences. In International Conference on Artificial Intelligence and Statistics, pp. 1608-1617. PMLR, 2018.
Siegfried Graf and Harald Luschgy. Foundations of quantization for probability distributions. Springer, 2007.
Nicholas Higham. Accuracy and slubilily of numerical algorithmss SIAM, 2002.
V Roshan Joseph, Tirthankar Dasgupta, Rui Tuo, and CF Jeff Wu. Sequential exploration of complex surfaces using minimum energy designs. Technometrics, 57(1):64-74, 2015.
Benoît Kloeckner. Approximation by finitely supported measures. ESAIM. Control, oplimisation and calculus of variations, 18(2):343-38) -59. 2012. ISSN 1292-8119.
Tianyi Lin, Nhat Ho, and Michael Jordan. On the efficiency of entropic regularized algorithms for optimal transport. Journal of Machine Learning Research, 23(137):1 42, 2022.
Dong C Liu and Jorge Nocodal. On the limited memory bigs method for large scale optimization. Mathematical programming, 45(1):503 528, 1989.
Mario Lucic, Matthew Faulkner, Andreas Krause, and Dan Feldman. Training gaussian mixture models scale via coresets. orv. 2018. ISSN 2331-8422.
20

---

Published in Transactions on Machine Learning Research (02/2025)
Giulia Luise, Alessandro Rudi, Massimiliano Pontil, and Carlo Ciliberto. Differential properties of sinkhorn approximation for learning with wasserstein distance. Advances in Neural Information Processing Systems, 31, 2018.
James MacQueen et al. Some methods for classification and analysis of multivariate observations. In Proceedings of the fifth Berkeley symposium mathematical stotistics and probability, volume 1, pp. 281-297. Oakland, CA, USA, 1967.
Simon Mak and Roshan Joseph. Support points. The Annals of statistics, 46(6A))2562-255) 2018. ISSN 0090-5364.
Albert W Marshall, Ingram Olkin, and Barry C Arnold. Inequalities: Theory of Majorization and Its Applications, Second Edilion. Springer Science & Business Media, 2011.
Jorge Nocedal and Stephen J Wright. Numerical optimization, 2nd edition. Springer, 2006.
Jorge Nocedal, Annick Sartenaer, and Ciyou Zhu. On the behavior of the gradient norm in the steepest descent method. Computational Optimization and Applications, 22(1):5 35, 2002.
Ofir Pele and Michael Werman. Fast and robust earth mover's distances. 2009 IEEE 12th international conference on computer vision, pp. 460-467. IEEE, 2009.
Gabriel Peyré, Marco Cuturi, et al. Computational optimal transport: With applications to data science. Foundations and Trends@ in Machine Learning, 11(5-6):355-077, 2019.
R. Tyrrell Rockafellar. Couvea Analysis. Princeton University Press, 1970.
Seref Sagiroglu and Duygu Sinanc. Big data: A review. In 2013 International Conference on Collaboration Technologies and Systems (CTS), pp. 42 47. IEEE, 2013. ISBN 1467364037.
Filippo Santambrogio. Optimal transport for applied mathematicians. Birkäuser, NY, 55(58-63)) 63).94 2015.
Bernhard Schmitzer. Stabilized sparse scaling algorithms for entropy regularized transport problems. SIAM Journal on Scientific Computing, (1(3):A1443-A144 2019.
Ozan Sener and Silvio Savarese. Active learning for convolutional neural networks: A core-set approach. arXiv.org, 2018. ISSN 2331-8422.
Richard Sinkhorn. A relationship between arbitrary positive matrices and doubly stochastic matrices. The annals of mathematical statistics, 35()):877-887 879, 1964.
Richard Sinkhorn and Paul Knopp. Concerning nonnegative matrices and doubly stochastic matrices. Pacific Jowrnal of Mathemalics, 21(2):343-3 348, 1967.
Gui-Xian Tian and Ting--Zu Huang. Inequalities for the minimum eigenvalue of M-matrices. The Electronic Journal of Linear Algebra, 20:291-302, 2010.
Cédric Villani. Optimal transport: old and new, volume 338. Springer, 2009.
HaiYing Wang, Rong Zhu, and Ping Ma. Optimal subsampling for large sample logistic regression. Journal of the American Statistical Association, (182) 844, 2018
Yujia Xie, Hanjun Dai, Minshuo Chen, Bo Dai, Tuo Zhao, Hongyuan Zha, Wei Wei, and Tomas Pfister. Differentiable top-k with optimal transport. In Advances in Neural Information Processing Systems 33, 2020.
21

---

Published in Transactions on Machine Learning Research (02/2025)
Appendix
A Explicit expressions for constants
We first define few user constants for the L-BFGS algorithm. Let my be the maximum number of correction vectors used to construct the BFGS matrix BR, and Ci E (0,1/2), e2 E (c1, 1) are two constants related to the Wolfe condition: we assume the L-BFGS algorithm uses some line search algorithm to select the step sizes ag that satisfy:
f(xx + axds) S f(xx) + clarg dx, g( + +add)) d 2 c2 de.
where f(.) and g(-) stand for the objective function and gradient function, respectively, Iy is the k-th iterate, 9k = g(zx), dx = =B 9k is the search direction, and By is the BFGS matrix that approximates the Hessian matrix. mu, C1, and 02 are selectod by the user.
For Theorem 5.2 let ((0) be the initial value, and let p Mra and n4 max, Mij: i == 1, ..... n. Then we define the following constants:
( max
U. = bml
15/5m 1// + -
11)+
aj log ai - --
1=1
m- 1
A = x- loga; 0; - 2-1 e A(MintUa
Σ
-AMes
1, .....,n
j=1
M1 X. + 2
min A(A.-Mam) 1si5n
M2 > 1 Σ
n
A(A.- Mim)
2(n - m + 1)
M3 = M₂ log M1 1,
=4 m m-= + m2 M2 - mg [log M1 - log(1 + M2))
C = 1/M C2 = 2M, M2
r= c (1 M / MM2e - MM)+
For Thcorem 5.3
Cs = V | + C.)]
C. 20 + 3 o w.) +3
C₂ || || + 2C,, diag(a ")²" || ||a o s,)l
B Additional experiment details
B.1 Settings of the motivating example
Consider a small problem of n = 90 and m = 60. Let x; = 5(i - 1)/(n - 1), i = L,...,, n be equally- spaced points on [0, 5], and similarly define y/ = 5(j - 1)/(m - 1), j
j= = L,..., m. The cost matrix is set to Myj (x; - ,2 and the weights a and are specified as follows. Let fi be the density function of an exponential distribution with mean 1, and fa be the density function of a mixture of two normal 0.25).
distributions, 0.2 (1, 0.04) + 0.8 0.250. And then we set (I,), b, faly;), a; = à Σk= de. and by = b, / b be.
For example showed in figure I We fix the regularization parameter - to be 0.001. This value is selected such that the resulting Sinkhorn plan T is visually close to the OT plan p The maximum number of iterations is 10000 for Greenkhorn and 1000 for other algorithms. Other
In Figure 10 we show the Wasserstein and Sinkhorn transport plans under different values of A It can be seen that. when A-1 S 0.001, Tx is visually indistinguishable from P*.
22

---

Published in Transactions on Machine Learning Research (02/2025)
40
79
50
Wasserstein
Sinkhorn 0.1
Sinkhorn 0.01
Sinkhorn 0.001
Sinkhorn 0.0001
Figure 10: Visualization of Sinkhorn plans with different X values.
We compute the true Tx using the E-scaling algorithm (Algorithm 3 of |Schmitzer (009)). This algorithm is typically accurate, but it requires solving sequence of Sinkhorn problems with increasing X's, where the solution corresponding to the previous À is used as a warm start for the next one. Therefore, its computational cost is typically large, and it does not compare fairly with other methods. Due to this reason, in this article we mainly use the E-scaling algorithm to compute high-precision reference values, and do not include it for method comparison.
B.2 Sinkhorn Loss Computation stability and accuracy
To further compare the numerical stability and accuracy of different algorithms for computing the Sinkhorn loss, we consider the following experiment. First, we simulate data XI...., X, ~ fi(x) and Yi...., Yn f2(y) from some specific distributions (x) and f2(y), x, y E RR, and construct the cost matrix as M = (My)) where Mj = || 12 The weights are fixed as a = n ',, and compared the motivating example in Section 4.1 let T (Tiz) be the computed Sinkhorn transport plan and 7 be the true value (computed using the e-scaling algorithm), and then we compute two types of errors: the error the transport plan,
mm' For each of the five methods
Errptan (T)
Σ
(Tej TT2)², s,j
and the error on the Sinkhorn loss value,
Eriioss (T) = (T, (T" MY = -T- -T",
For each configuration of the experiment, we simulate the data 100 times, and visualize the distribution of the errors using boxplots.
In our experiment, we fix n = 150, m = 200, and consider varying dimensions p= 1, 10, 50. The Sinkhorn regularization parameters compared are - 0.01, 0.001, and for cach X-1 we set specific maximum number of iterations for all algorithms. Two data generation models considered:
(a) Both fz( and fz(y) are multivariate normal distributions NN0, 1p):
(b) Both fi) and fy) have independent components. Each marginal distribution of fi is an exponential distribution with mean 1, and each marginal distribution of f2 is mixture of two normal distributions, 0.2. (1, 0.04) + 0.8 (..00.5).
The final results demonstrated in Figure 11 where all the errors shown in the log-scale.
23

---

Published in Transactions on Machine Learning Research (02/2025)
Errors on the Sinkhorn transport plan
Dim: p
"pf--
Distribution: normal
-2
Log scale Emor
n/Maxite
0.01/300
0.01/1000
Distribution:: exp mxx
0.001/600
0.001/2000
5
Sinkhomn
hint.tom-log
Statilized
Greenkhom
L-BFGS
Sinkhern
Sinihern-log
Stablized
Greenkhom
L-BFG5
Sinkhorn
Sinthamn-og
Stabilized
Greenkhom
L-BFGS
Algorithm
Sinkhorn transport plan
Errors on the Sinkhorn loss value
Dim: -1
Dim: 10
Dim: 50
2
0
Distribution: normal
Log Log-scale Error Enor
n/Maxiter
-8
0.01/300
0.01/1000
2
Distribution: exp mbx
0.001/600
0.001/2000
-2
Sinktern
Sinktorn-log
Stabiized
Greenkham
L-BFGS
Sinkter
Sinthern-log
Stablized
Greenkhom
L-BFGS
Sinkhom
Sinihorm-iog
Stablized
Greenkhom
L-BFGS
Algorithm
(b) Sinkhorn loss value
Figure 11: Comparing different algorithms on the errors of Sinkhorn transport plan and loss value. The missing boxplots indicate that the corresponding results are NaNs.
In Figure I many buxplots for the Stabilized and Greenkhorn algorithms are missing, since they produce all NaN values in the 100 simulations due to numerical overflows. For Sinkhorn, even if it generates no NaNs values explicitly for - 0.001, it does not give any meaningful results either. This implies that numerical stability a critical issue in computing the Sinkhorn loss.
For Sinkhorn-log, it gives reasonably small crrors when the regularization large (X-1 0.01) with sufficient number of iterations, but its accuracy quickly deteriorates when X-1 decreases to 0.001. Moreover, we can find that Sinkhorn-log is sensitive the limits number of iterations. For example, when p and
24

---

Published in Transactions on Machine Learning Research (02/2025)
X- 0.01, the loss value error can be small as 10-6 given 1000 maximum mumber of iterations, but if we restrict the limit to 200, the error can be as large as 10-2 even 10", depending on the data distribution.
In contrast, the difference on maximum number of iterations has minor effect on the L-BFGS algorithm, indicating that it converges fast, and additional iterations are not needed. These findings demonstrate the advantage of the advocated T-BFGS method in both numerical stability and accuracy.
B.3 Running time of Solution Stage and Differentiation Stage
In this section we compare the running time of three algorithms on differentiating the Sinkhorn loss.
(a) The proposed algorithm "Analytic": L-BFGS in the solution stage, and analytic differentiation in the differentiation stage.
(b) "Implicit" implemented in the OTT-JAX library (Cuturi et al. 2022 Sinkhorn-log in the solution stage, and implicit differentiation in the differentiation stage.
(c) "Unroll" implemented in the OTT-JAX library: Sinkhorn-log in the solution stage, and unrolled automatic differentiation in the differentiation stage.
We the second data generation model in Section B.2 simulate data points, and use the three algorithms above to compute the Sinkhorn loss and its derivative with respect to the cost matrix. For each configuration, randomly generale the data 100 times, and compute their mean solving time and mean total time. The stopping rule implemented in the OTT-JAX library is ||7,m all + |TT 1, - bll < Sott, and one of the terms would be exactly zero in the last iteration. The stopping rule for L-BFGS is ||T I, bllso < Etffas To account for such difference, we set Stbfg= = 10-6, and let 6ott / m] Stbfgs- In fact, this setting favors the competing method, its stopping criterion is strictly weaker than the proposed onc. To test whether the algorithms actually converge under the given criteria, we also report the number of converging cases within the 100 repetitions.
Results for different data dimensions and regularization parameters are given Table and Table 2 the former 1000 maximun number of iterations, and the latter uses 10000.
where
C Proofs of theorems
C.1 Technical Lemmas
In this section we state few technical lemmas that are used to prove our main theorems. Lemma 2 to Lemma f below are standard conclusions in vector caleulus, Lemma 5 and Lemma G are derived from the cigenvalue theory, and Lemma and Lemma g are related to the Sinkhorn problem.
(2. 0). For a matrix (a:y) E g=xm the vectorization operator, vec(A), creates vector by stacking the column vectors of A together, i.e.,
We first introduce the following notations. For any x IR, let [x]+ = max
vec( A) = (an,.... Aml, 012, Un2,
01m.
anm
For two matrices A (aiz) Raxm and B E Rexg, the Kronecker product of A and B is defined
B B
Qim B
A = =
:
On
apm B
For differentiable vector-valued function defined
m the partial derivative of with respect to is
(x)
afiz) az.
0f(x)
02
T
8 0/
0 (w)
8x
25

---

Published in Transactions on Machine Learning Research (02/2025)
Table Running time three algorithms for differentiating the Sinkhorn loss, with maximum 1000 iterations.
Mean solving
Std. of solving time
Mean total time (ms)
Std. of total time
Converged
time (ms)
Analytic
19.07
2.11
19.33
2.13
100
64
Implicit
27.77
0.69
39.66
3.41
A-1 0.1
Unroll
27.53
0.13
101.30
0.36
0
Analytic
23.21
3.45
23.48
100
m = r 64
=
Implicit
27.33
0.05
50.52
5.19
0
A-1 0.01
Unroll
27.39
0.12
100.66
0.39
0
Analytic
23.79
2.96
24.19
2.96
100
m n 128
= 16
Implicit
29.81
0.61
58.21
4.91
0
X- 0.1
Unroll
29.75
0.06
107.18
0.19
Analytic
32.24
4.99
32.68
5.00
100
128
= X- 0.01
Implicit
29.73
0.05
85.56
7.48
Unroll
29.65
0.08
106.66
0.46
0
Analytic
28.60
3.82
29.40
3.85
100
256 p
Implicit
33.69
0.32
82.53
4.21
0
A
0.1
Unroll
33.60
0.23
123.25
0.07
0
Analytic
43.94
6.08
44.86
6.12
100
m = n 256
32 A 0.01
Implicit
33.75
0.20
137.46
11.66
0
Unroll
33.78
0.10
122.99
0.27
Analytic
41.35
5.64
43.05
5.53
100
m 7
= 61
Implicit
46.06
0.05
130.96
11.97
0
_ 0.1
Unroll
46.02
0.04
167.64
0.22
Analytic
67.76
9.65
69.42
9.65
100
512
p =
Implicit
46.08
230.34
11.72
X 0.01
Unroll
46.09
0.05
167.68
0.20
0
We use J to denote the n x n identity matrix, and (.ax(.) and Tmia stand for the largest and smallest eigenvalues of some symmetric matrix, respectively.
Lemma 2. Given two matrices A E pmxx and B R
vec(AB) (I, @ A)vec(B) = (BT @ Im)vec(A).
Lemma 3. Let f 1 + RD and R + R" be two vector- valued differentiable functions of x E R=.
Then
d (1) (T)))(<) g(x) 0f(z)
[(x)T 0g(z)
T
Ox.
Lemma 4. Let f Rm + R and 1m + Rr be two vector-valued differentiable functions of e R
Then
i veel (()((()))()))) (g(x) 1)) 0f(x)
d.T + (I, Φ f(x)) Oglx)
T
Lemma 5. Let A and B be two x n positive definite matrices, and let a1 2 ...
2 an > 0 and B1 2 ...
B, be the ordered eigenvalues of A and B, respectively. Then for any x E RR,
T 1A1/22 11
26

---

Published in Transactions on Machine Learning Research (02/2025)
Table 2: Running time of three algorithms for differentiating the Sinkhorn loss, with maximum 10000 iterations.
Mean solving
Std. of solving
Mean total time (ms)
Std. of total time
Converged
time (ms)
time
Analytic
19.01
2.10
19.27
2.10
100
64
Implicit
267.95
24.57
277.86
24.58
A-1 0.1
Unroll
269.49
24.76
986.46
90.67
12
Analytic
23.28
3.46
23.54
3.48
100
m = r 64
=
Implicit
273.42
0.34
284.84
4.74
0
A-1 0.01
Unroll
273.54
0.49
1011.70
0.75
Analytic
23.76
2.92
24.16
2.90
100
m n 128
= 16
Implicit
297.36
0.34
320.29
4.32
0
X- 0.1
Unroll
297.93
0.63
1074.47
1.33
Analytic
32.21
5.00
32.68
5.01
100
=
Implicit
297.79
0.30
324.95
X- 0.01
Unroll
297.90
0.32
1073.24
0.84
0
Analytic
28.72
3.84
29.57
3.83
100
256 p
Implicit
342.23
1.53
386.80
9.22
0
0.1 A
Unroll
343.18
0.70
1256.10
2.94
0
Analytic
44.45
6.14
45.32
6.18
100
m = n 256
32 A 0.01
Implicit
343.57
2.01
411.12
18.92
0
Unroll
341.00
2.02
1253.67
3.33
Analytic
41.50
5.60
43.27
5.53
100
m 7
= 64
Implicit
469.04
0.77
552.64
15.31
0
_ 0.1
Unroll
465.03
1.82
1705.51
6.36
Analytic
67.97
9.69
69.68
9.67
100
512
p =
Implicit
464.59
624.28
26.48
x-1
0.01
Unroll
464.16
0.55
1707.04
3.31
0
Proof. Let Uixn x A1/2, and then = UUT xFAr Ax < an 1x]3, and UBUT A/22//2/ By Theorem A.4 in page 788 of Marshall et al. 2011 we immediately get
UBUT tr(UBUT) S Biu S 12112
Lemma 6. Let and be two symmelric matrices of the same size. Then
Omin (A) + umin (B) S Omial (A B) S m (A + B) S T (A) Omnx (B).
Proof. Using the well-known identity Fmax (A) max ia --1 x Axx we have
Omax Omax(A+ + B) = max x =T(A+ + B) max x" T Ax + max ||]|==
TBr = Omax(A) + T (B).
Lar 1-1
||||=
Applying the inequality above -(A B), we get the result on the opposite direction.
27

---

Published in Transactions on Machine Learning Research (02/2025)
Lemma 7. Let (33 be defined as in 12 and let μs M a and u, max, Maj, i= 1, ................................ ff sc
then we have max, B3 s Ua and miny $3 2 The where
Uc bm
max
Mi
0j log a;
-- +
(16)
1SySm-1
+
Le =
min 15ssm-1
bj
Σ a + - Σ
ag loga, a: - x +c
(17)
1=1
Proof. By definition,
f(3) = -- Σ aj log
e 100 M,,)
- --11
a log a; 6 + +
1=1
If f(3) S c, then
n
co = c+ Σ aj log 04 - -x-1 2 - Σ aj log
MB, ( Ma)
gT
By definition we have 8n 0, and let J arg maxiaajsmmmmm Bz. Then
m-1 1
8 a log ( - -3 = Σ (,3 - -M1)) Σ
Byb
j=1
m m-1
BJ Mi - -B3
Σ
bj
3y ps 3m
max P4:
15/5m-1
j=1
which verifies 16 by noting that max, Bg [3,]+.
Next, let K = arg mini minis)sm-111 By. We can assume that 3K < 0, since otherwise the trivial bound min, Bg = Bm = 0 is already met. Considor the sets S+ = {j: (j:3, B, > 0} and S_ = {j 33 < 0}. Then clearly,
8 = Bxbx= +
Σ B,bJ + Σ B,b₁ S 3x + 13,1+ Σ
by 0 < 8xbx + [851+ (1 - bm).
H jES+
J JES
i JESt
Also note that log
プØ e > max; J.y for any In R, so
X
0g log
Σ
MB8 Mis)
2 21 Σ
log
Σ
AlB uc)
x
log Σ eB3
Aua
t=1
1=1
1=1
1=1 l
---log
Σ ABs
ai44 2 --1 max (AB₁) - Σ
auu
1=1
n Σ
n
max B, 3
0 = [55]+ Σ
a;
1=1
As a result,
n
m
n
3o 2 111
a; log
e A18y
-33²b 2 185]+ - Σ
aæu; - 3xbx - -135]+ - bm.)
1=1
1mm113) [3,]+
Σ
a.44 - 8xbx 2
Σ
a,84 8x
1=1
1=1
28

---

Published in Transactions on Machine Learning Research (02/2025)
and then
n
Bx
-b
Σ
a
co
min 15/Sm-1
bg
Σ
a + Co
which verifies
17
Lemma 8. Let be an n m matria with strictly positive entries, and suppose that n 2 m. Define p T1m; In, and
(dag(())
T
11
T
diag(e)
D diag(D) TT -TTdiag(0 diag(u) 'T.
Then
max (1a 5
max Vjj 15/5m-1
m +
mmiu (D) 2 Omin H 2 2(n m + +)
min Tims 1si5
m-1
6 min (D) 2 min
Σ Dij 15jsm-1
min
Σ
14 Timm
13i5m-1
==1
1=1
Proof. Consider the matrix H sJ, where an (n + m 1) (n m - 1) matrix filled ones, and
positive scalar. Let
R
Σ
h == 1, 1...., + m 1. Skyl:
jz
Suppose S ISs T and then for k- = I,...,n, n, it easy to find that
m-1
R4 = (n-1)s+ +
Σ
(Tk 8) = (n 1)s - uh Tkm (m - 1)$ = (n m)s + th Tkm
j=1
and for k = n + ,,..., m
Rs = (m- 2)s +
Σ
(Ti 8) = 2)s + Vk-n ns = (m - n - 2)8 + Ok
1=1
Then it is easy to see that
R4 km (n ms)s k 1, ......7 Rs (n + 2 m)s,
Skk
R4
Vk-n
k = n + 1,
m -
Let
min Tim - (n - m)s = (n + 2 -m)8,
15i5n
and then minisis, Tmm/(2n) 2m + 2), and
- m7 2
SEL Rs 2 L =
min Tim. > 0
2(n +
13i5n
for all k. By the Gershgorin circle theorem, every eigenvalue of S must be greater than L. Since H = S + SJ and J is nonnegative definite, we also have amin (H) > 0, implying that. H is positive definite.
For the sccond formula, it casy find that the matrix the Schur complement of the block diag(u) of the H matrix. So by Theorem 3.1 of Fan (2002), we have Tmin (D) 2 J (H) and F (D) S Jmax (diag()) maxisjsm-11 Vj-
29

---

Published in Transactions on Machine Learning Research (02/2025)
Finally, let maxig/Sm-1 My, and then D can be expressed D clm-1-B, B, where B = T diag(u)) diag(c1m- m-1 2) is a matrix that have nonnegative entries. In addition, we have proved that D is positive definite, 50 D is a nonsingular M-matrix by the definition in Tian & Huang (2010). Then Theorem 3.2 of Tian & Huang (2010) shows that
m-1
min (D) 2
min
Σ Dis- j=1
Isism
Let 8 = D1m-1, and then clearly mini minisiim-1
m 2j=1 Dij min; 8, Note that
8 = D1m-1 T diag((( 17' 1 1 - T d'dag(()) ( - mT
- diag(()) Tm) T 1, + 7 diag(())))
"T -T "Tm ( -TTm
where Tm stands for the m-th column of T. Therefore,
n Σ 14 Tii Tim- 15755-
min ISSm-
5, min
C.2 Proof of Lemma 1
Let = ex + 3 - M), and then it easy to find that VaL(a, 3)= = a Tim and VzL(a, = T 1n Since a" (B) arg maxa [(a, B), we find that a; (3); is the solution to the equation a T1m 0. By definition, we have
m
Moc+8,-Maj.
da:
AB,- (18--M3 Ms), ,,...,n,
0:
j=1
j=1
so the solution is a; = x- log - 101log EE/ M(BB Mas) 1, so
Since T1m = a, we immediately get TT,,
(()) = (²) a + 8 x 1 a (B) -
and we get the expression for f(3) L(aa" (B), 3). Finally,
(3) (a.((a,))aaar ViL(a, 33 a=a*(
Vg((a'(3)
03T
Since (3) = arg max, L(a, 3) implies that Valla, = 0, we have Vz((a*(3),3)
V Vgl(a, B) aa)*(()) and hence
V8f(3)
C.3 Proof of Theorem 4.1
By definition we have
Sx(1,)) = (T*, = "ee(TT")Tvec(A
so Lemma 3 gives
(S
vec(M)) dvec(T")
Ovec( + vec(TT)) dvec(M)
(18)
Ovec(M)
Ovec MYT
Obviously, (rec(M))Ovec((M)) is the (nm) (nm) identity matrix [(mm). the second term of 18 is
essentially vec(T")" and the remaining task to derive
Let R = a* 0 - M = a IT + 1T - M, and then T = [R Using the chain rule of derivatives, we have
Ovec(T") )vec(T") Ōvec(R)
(19)
@vec(M)T dvec(RyT @vec(M)T*
30

---

Published in Transactions on Machine Learning Research (02/2025)
It is easy to find that @vec(T")/Ovec((R)) is an (nm) (nm) diagonal matrix with diagonal elements vec(XT"),
voc(M)) Övec(T") Övec(R)T
= Avec(M o TT)².
(20)
Furthermore,
Ōvec(R)
@vec(a" 1m + 1,, - M)
@vec(M))
Ōvec(M)²
03* Óvec(M)) + (Im & 1n)
Oa
(1mm @ In)
[(mm)
(21)
(V)
where the second identity is an application of Lemma Combine
19
20 and
21 and then we get
vec(M) (vec(T") @vec(M)T
Xvec(M o TT) @vec(R)
Ōvec( (M)
da*
@ In)
Ovec MyT
03*
>vvecMM @ (T @ 1n)
Övec(M)T
Xveec (T))
which is the first term of 18
Using the identities in Lemma
we have
vec(M o T & In) = [1mm @ vec( o T*) T = [1 L.)ec(M
= μr
vec(M 2 @ TT' 2 1n) = [Imm @ 1,,) " vec o TT]== = [Im @ vec
= Me [ T
Since we have set 3r = 0, 18 simplifies to
0S*(H,)) v) Övec(M)T
da
a
(rr Övec(MyT + aTT
Övec(M)T vec(M TT)) veo(TT)T
(22)
Let w² (a*T, (TTTT and then the main challenge is to calculate Ow" /ovec(M)T.
First, note that the optimality condition for (a" , = arg maxa,3 L(a, is
V B,L(a,B))
O, VaL(a, () ((
0.
(23)
Section C.2 has shown that Val(a, 3) = a - Tim and V VaLía, 3) = TT
1n. Moreover,
V2,Lla, B3 Adiag(T1,))
Vzlla, B) = - Adiag(TT
V (Vo((a, = -XĀ.
Define the function
Flw, V.Lla. ViL(a,
T1 TT,,
(24)
where w = (aT BTYT, and then w* satisfies the equation F(w", = o, indicating that wt is implicitly a function of M, written as w = w(M). By the implicit function theorem,
Ow(M) -F- FM FM =
aF M)
OF w. M)
Dvec(M)T
TwT
Ovec MTT
"=w"
31

---

Published in Transactions on Machine Learning Research (02/2025)
Note that
(
(
BT B)=
Fw FT
[v8 (V.C(a,)))
vzl(a,
Then by the inversion formula for block matrices, we have
(
A B
B) =-1-1 ( A A BÂ
B D
BTA
A
-A -A-BBÂ-
)
=
- -Â-1 1TTAAA
where Á = D - ATA-LB. B. For 9 (F)=,)) the vector 5 (8 ST) expression:
- -AF has the following
5, = - #4 + Â Acc
8u A Mr A BÂ BT A pr A- Pe.
= -A Bs.
After some simplification, we obtain
Å = diag( [ 1n) - T diag (iaag((mm) (T'I)
8, Å Phe Á 311f"(((Tmm))) O (t)
Su = (TIm) 3 Hr (T1m) о (TBy).
Gm H
Next, partition Fas as Fas
where GM E T and HM E (m ())mm)))(mm) By definition,
8G. Tvecl MyT
GM
G =
Ti
OG
j=1
j=1
@ve
S0
aG,
0,
i
0M
ATRI. i=
This indicates that Gm X (diag(Ti), diag((m))) where Ti...., Tm the column vectors of T. Similarly, for Hj Eim1 Tij.
0H avec(
OH,
0.
j XTRI,
HM-
0MH
j=1 1'
0Hm-1 AF)
implying that
Tf
HM =
TT
m- 0
As result,
da- Htth Tve + 1/t &vec(M))
03-
(ut,ptt) Ow
(1it M.. ) FM
Övec(M))
X GM
x- HM ) = -sTT, sT)) (
X G
XF
Pr
x-1 Has
((sn & TT) ,..., (ou ④ Tm 17) + (5o, T.. :..., 80,m TT - 0
[vec (diag(s.)) + Tdiag( [T)))
(25)
32

---

Published in Transactions on Machine Learning Research (02/2025)
Finally, substitute 25 into 22 and we have
(S = |vec (u + T* didag(()))) - ( + vec( @vec (M)T
Transforming back the matrix form, and we obtain
@Sx(p.)) = Ms+ - 0M
Replacing T with T and noting that T my b .T
get the stated result. The positive
definiteness of the Å matrix is direct consequence of Lemma
C.4 Proof of Theorem
In the proof of Thcorem have alrcady shown that
(diag(T1 (1)m 7T
T diag(TTT,,) )
V2 a C(a B) -HH = =
where T ex(a 3 MM. Plugging a" (3) into L(a, B), and then (3((a(()) is the Schur complement of the top left block of V2 v. f(a, B), given by
VzLla (3), - - [diag 11)) 1,,) - TT
Since f(3) = -((a (3), 3), by Lemma we find that (25(5) positive definite, so f(3) strictly convex on 3, and hence 3 is unique.
The optimality conditions for (a*, ,88 are T11, and 7 b, where T exlaa @ fr M). Since T exp{ β5 M1j)] 2 0 and a; = -ズ= Tij. bj Σi= Tīj, we have Tī minfa, b,} for all and
j. implying that
a 5 = Ms + x min (log(a:). log(b,))
Since By = 0 by design, we have of S Ua, +00, = 1,...,n, n. where Ua Umm - min(moo((ee (a;), log(bm)) This indicates that at upper bounded.
Mim
Next, let I = arg max; Tim Since T 1,, b implies that bm Ei=1 Tú nIimm we have
Tim exp[M(af M1m)) 2 bmn/n,
and hence af 2 Mim + 1 log(bm/ /n). Again, since ot + 33 s Usi for all i and j, it holds that
8 a7 < U1j Mi - log(bmm/ (f/m//mnn) /n)
Mrj + x min (onflog(ur), log(b;)) - Mrm x log(bm//))
= Míj Mr + 1 min (oin(log(vuul (nar /bm), (og((nbs/bm)))
= UR, < +00, j = I....,m. m.
(26)
On the other hand, TT b implios that b, = EE The = - En Maa MM]-Mu) for any j. so
logb, A33 log Σ Ma; Mb,)
$ SA87 + log
Σ AU.:--Mi
It is well-known that
max maxf21,...,.2 xn S log
Σ 6
< max max(f1....,.
2.. log
33

---

Published in Transactions on Machine Learning Research (02/2025)
for any 21,..., In R,
n A(Un; Ales)
83 X log bj x- log
2 - log bj- - max (Ua. - Mij) x logn
Lg == x- (og(b,)(b))) max (Vas M2j) > -00, j=
(27)
Then 26 and 27
together show that |871 < xc.
Similarly, mm a, so
logar a; S Aa; + log Σ AeUny -My)
of Luu; x- log(as/m) max (Ug, Maj) > Oc, i= 1=1,...;1
The trivial bounds La,, and Ua, are obtained by removing the unknown index 7.
The results above verify that lail < oo and [3j| < co, and hence ||a*|| < oo and [|8*|| < co. Finally, plugging in f to the objoctive function, and we immodiately get f > -o0.
C.5 Proof of Theorem 5.2
Claims (a) and (b) are direct consoquences of the convergence property of the L-BFGS algorithm (Theorem 7.1, Liu Nocedal (1989)) and we only need to verily its three assumptions. The new results here are explicit expressions for the constants C1, C2, and
First, is twice continuously differentiable, so Assumption (1) of Liu & Nocedal 1989) is verified. Second, [ is closed convex function, and we define the level set of f as Le {B E Rm=
::ma-T::1 ((3) < c). Theorem
5.1 has shown that f > -00, and when f*, obviously L. = { is non-empty and bounded. Then Corollary 8.7.1 of RRockafellar (1970) shows that Lg bounded for every c. In particular, for a fixed initial value ((0), define L {3 f(3) S { and then L is a bounded, closed, and convex set, which verifies Assumption (2) of Tiu Nocedal 1989). Third, let H(3) = (3)(),, and then the proof of Theorem
5.1 we have already shown that
H [iag (T
1,,) 1,)) 7TH diag (1,m (
where T cx/a (3) M]. Lemma verifies that
m (H(3)) 2 X. 2(n m + 1)
min Tim 18
Tmax mmax(H(B)) S X
max 1/j, 15i5m1
where v In, On the L sct, Lemma shows that max, Bg $ Ue and miny 35 2 Les with c = f(3(0).
Therefore,
0; = -- lollo - A-'log
MB3 - M,y) Σ 1=1
> --1 log a; - >-'log e AMem +
Σ
e A(Uc-M,,)
j=1
m-1
=A, = - log a; - Ua --
e A
+
e
-XMi5
1=1
Tim = Mla--Mull 2 M
On the other hand, Ti must satisfy Ti > 0 and have Tig ag Tim Therefore,
EØ= Tıj = a; for any i and j, so for j = 1,..., m 1, we
Tim) = Σ 0ラ = Σ Tí S Σ (a:
Ti S Σ é M(d Min),
== 11,...,,m
T=
i=1
i=1
34

---

Published in Transactions on Machine Learning Research (02/2025)
This implies that there exist constants M1, M2 >> 0 such that
Milka S x² (B)x S M2|1a112
for all r E Rm-1 and B E L, with
Mi = X.
11 - m +2
min e A(As Mem) 15:5n
2(n m + 1)
M2 = A 1 CA(Ai- - Mem)
Am
This verifies Assumption 7.1(3) of Liu & Nocedal 1989)
Next, we derive the explicit constants in the theorem. Following the notations in equation (7.3) of Liu Nocodal 1989), the BFGS matrix Bş for the L-BFGS algorithm has the following expression
B₁ B(m, B(++) = = B B
y19i
8f B(ss
vf5st 8t
where m = m + 1, mo), me is a user-defined constant explained in Section A and {wr} and {si} are two soquences of vectors. We also choose B(0) to be the idcntity matrix.
Fix m, and let cos O2 of B st/([ss (1BW)ll)) of (()ss// ||>s| Then it be verified that
Pk: vi st /1/1111113 Th 11aa11/// 81, and qu
|B00 1r
tr(B(+1)) = tr(B()
+
tr tr(B())
Ts: cos? 0x
+ Tk-
(28)
87 B()ss
87's1
dct(B((+1) = det dcc(B())) (B()))(((ee) pp//
(29)
Define vb( B) tr(B)- det [B,, and it is known that (BB) > 0 for any positive definite matrix B. Equation (6.50) of Nocedal & Wright (2006) shows that
< (B(+1) tr(B())
Qk
Tk dct log Pk + log ql
cos 8g
=B)( + (Th- log ph -1) (1- cos² 0g
+ log cos2 0g
+ log cos? Os-
Under the assumptions verified above, equations (7.8) and (7.9) of Liu & Nocedal 1989 show that M1 S of ///s/|/ /|/l2 < M2 and ||ux|2/bj 91 & M2 for every I, so Mi I Pk I M2 and Th S M2. Since h(x)
x + log(x)= < for all x > 0, we have
(Bx) + M3 log cos² 0
(30)
where Ma M2 log Mi Now we show that (Bx) can be upper bounded. First, Lemma 6 that for 1 = 0, ..... m 1,
implies
max (((+)) maax (B)) +0 ||ull///7T st S mBx + M2,
so mmax (Bx) < 1 + M2 This also implies that da - 1 + my. Next,
28
shows that tr (r(BB) S
m 1 + mo and (29 gives
log deget(B's (let(B((+))) 1)) = log (c)( + log - log qk 2 log det( + log M1 log(1 + my)
implying that log det (B)) 2 mi [log Mi - log(1 + m2) As a result, we gel
(Bx)) S M1 m mo mo [log M1 - log(1 + mgM2))
(31)
Combining 30 and 31 we have cos² 0g e Mc)+
35

---

Published in Transactions on Machine Learning Research (02/2025)
Finally, using the argument in Byrd et al. (1887) have f(k+1) - ff* -f*), where
r= -=1-c (1-c)))) cos 0₂ < 11 Ma0ll-MM)MMMMe (MMs++MM)
and C1, 02 two constants for the Wolfe condition explained in Section 2/M:. For (c), we follow the analysis in Nocedal et al. 20002). Let g(3 = V (f), and then g(5*) = 0. By Taylor's
The constant C1 is simply
theorem we have
(()) - f 2 (3 3*) 3*);
(32)
where H1 (E) for some & the line segment connecting 3 and 3 Also,
g - - 9(3)) = H2(B - 3²). H2 [ H + t(B - -
(33)
Combining 32
and 33
we get
||ا(3)|1 = - H 2 - B-)
(f(B) - f*)
(3 (3 3r)
for any 3 L. Let x 23 Br), and then
||r(3)11 2x 11, 1/2 11 11 1/2
(J( --)-
|
It is easy to find that
maax (1T m
min (H1)) S. MT max (12) m1 max (112)1 S M2.
By Lemma we have x (f(33 for all B E L, where C2 = 2M M2. Since f(k) I f(o) due to claim (a), we find that
1/2 < M 1 for any x E Rm-1 and hence [[)(3)] <
g(k) L for all k O. Therefore,
5 - 8r) s (r10 -8),
and claim (c) is proved.
Claims (d) and (c) can be verified follows. For any B L, recall that T= c8) MM] and then T1m 0 and
[T 1n bll00 I TT - bl [af(B)]l
indicating that
n
Tτ bi <1Va5(3311. j= 1.....,
(34)
Then 0 Ti S プド Ti ag and 0 < T. S Ty
ff(]]] V
The gradient [|V]f( can be
bounded using claim (c), (d) is also proved. On the other hand, 34 shows that
0
T4 S m max Toj,
j=l
and similarly we have
by - -V S Σ To S n max Tij-
Replacing ||V ((3) by its upper bound, and claim (e) is verified.
36

---

Published in Transactions on Machine Learning Research (02/2025)
C.6 Proof of Theorem 5,3
For matrix Anxm (at)), define ||A||>> maxig luyll and the notation A 2 0 means aij 2 0 for all i and j. First, it is easy to show that ||A @ B|i - || [l] since
I B
(aisbiz)²
||AI=
Next, we show that if Baxm 0, Cpxn and 0, where is vector, then [IC(AQ ||||s>
Proof. Let us be the i-th element of (A @ B)v, and then
m 02
m
aiy Biy
|ud|
s ||||| 20
Σ
bygg.
j=1
j=1
Consequently,
p
2
n
2
y
m
2
|CCA @ Byul
Chilli
Ckil uil
S || az
Σ
Cki
Σ
big
j=1
bo ||CBll
=
Similarly, if 12220 and
0, then
A(
Ujj
S ||
Σ
Uz
=
| |
Moreover, for matrices Buxm and Cmxp 0, let Cj be the j-th column of C, and then
p
p
I B)C
I(A B)C,I' &Al
|B, ||A |BC||
Let + f and B = 3r + for some perturbation vectors f and g, and define T cxla --MM] Then it is easy to find that
T-T==exle + 5* M + (f - = T
Let E, = cxlf @g 1, so - T = Er o T Since 1e" - 1| < 2| for |æ| < 1, we have
(Er) = (()()))))) - < 2X|,f + gýI
long Alfi + 93| < 1. This can be achieved by assuming 2X( fl + [|g]lo)) < 1, since in this case Alf. + 9;1 < fll + (1glloo) < 1/2. Then clearly | <.
Consider 8, = 8n + 5n and 8, 8, + 8, for some perturbation vectors 8, and 8,, and let.
VMS T + X(S 8 GT
T X(s, @ s, M) o T 165, O
VMS+ ET 3 V 1(o @ (T (T
37

---

Published in Transactions on Machine Learning Research (02/2025)
Then we have
S SMSSFF
1 Er O V [6 06.) о (T + Ey @ I
doll Er o l
< < -rlls
$
< 5 X(1 + 3 1
< V + 2 7 [1 [1
+ ((0olo0)
Therefore, we just need show proper bounds for || and |16ullo.
Consider fr, (M o 1m and Aic (M о 1) and then
Pr = (M @ (T* + Er O T ) = pr + == pr + (Er O M M 1
ce (M o (T" + Er o T*)) " Me + te = Me (Er O
It can be easily verified that
M /bll ) ( < <s aa ② d = O
(a o 6,]1| = diag(a (((ag(a)((Er)) @ M o 1m s | diag(a (M @ 1mll
< < o pr)l:
||&.1| = [| @ M & T*) 1m | || o Im < <elluel-
Define by Me TT @ pr) and b. f (a (r), and have
by = Me 6. (T + ET T T( @ (v++ 6,))
= He 8e T (a- o (pr 8,)) (Er o o (M, +,))
= Mh T (a p - T 0,) (Er 2 Ta @ (px+ 5,))
As result,
16, bull = 160 - T 1 +,))
< 6elll + 06)1 0 + +e| 26,1
(1, +,)))
< <el|uall (c 0 +
< <E||mell +E|lT o p>)l + @ + e (a @.,)
<El|mell + Ta o ,
On the other hand, let byg be the (i,j) element of the matrix TT diag(a )Ý, and tij be the (i,j) clement of 7-T diTaaa(("")) "T* Then
Eg
ThiTy aa Σ + (ET) (Tr)
k=1
n
[6ij ligl
> [(TT (ET)aj + (FT) k=1
< (2e + s2)
T < 3Etig.
k=1
This implies that
| - D| = =117*T diag(a -/ -TT TTTdaag(a a 1|l|
-Iag'diag((
Σ
1.g12 s 3ell D|
i.
38

---

Published in Transactions on Machine Learning Research (02/2025)
Then by Theorem 7.2 of Higham (2002). if 3EO D|l 1, we have
|6,1
80
|ell + (3 o @x+)l + +3IDD))
S
|
3ea D||
||801|
where o -| |D- ol = 11/0mia() Assume that & << min 1/(6o) D((e),, then with slight simplification,
we have | loc I 11o.1 S Cve, where
C. 33 @ + + +3)
On the other hand,
16011100 S ||8,11 < ||a
② 6,
+ |a 3 (Ts, s=)
Σ =||u,|||
o T ]8 s,)l
X 1||m||
④ (T" o ((ET : 7" )8,)l
S |||m~|
(T* |diag( a (ET @ 8.
S E||æ~||
(T* =||diag( sell
|||/~||
(T* +e||dia(( yT + 6,)1l @
S æ||μ~|| + (1++))
3 ( ||ddaga "T""
+ diag
S +|| +1 (1
+ =|la= @ (T*s,))
S al1 + 2C, diag
| > ||a @ (T")s)))).
Combining the results together, get
V SM (n +
where C.. ||,||| + /diaga a ) ||a o s.)l
Finally, Theorem (.2(b) shows that |19||0 |F ||0(k) B8 S (ie(k) In addition, by
11
we have
af (k) = a a3((0)) x- loga; -- log m
-co(t) Mes)
Then using Lemma (??), we obtain
|a - -a"l --
log Σ
(60 - llog Σ c MB Myy)
_X-"1><(800 - -3,) = -1glloc
1=1
implying that |x = |1a0 - a0 S 1lo1lo S V(e((t). As a result, = 2X(1]]1 + (1gllx)) S
4Xv (1())
C.7 Proof of Theorem 5..
Lemma 9. Let (X, d) be Polish space, and p (1). Denote the space of Borel probability measures
X as P(X'). The Wasserstein space of order p is defined as Pp(X)= {u P(X); fx d(zo, x)2 u(dx) < oo), where E X is arbitrary. Let (μu E ()(N))) be a sequence of probability measure, and E Py(X). If Wp(,, + 0, then pkk 4 p.
This lemma a direct result of Theorem 6.9 in Villani (2099).
Lemma 10. Let AN = {uNlun = EN wixx, VX, E X w, 2 0. EN = 1}, meaning the set of all
probability measures on X' supported on at most N points. Then
min Wo(u, + ANEAN
39

---

Published in Transactions on Machine Learning Research (02/2025)
This lemma is proved in Theorem 1.1 of Kloeckner 2012
Now we use the above two lemmas to prove Theorem 5.4
Let Py= = argmin Wp(n, PN), PB = argmin (p(
WNEAN
PBEAg
By Lemma 10 we have (x, + O. Therefore, Ve > 0, NN E N, st. Wy(,, < c/22. Note that with B large enough, we can always construct μB E AB to approximate PN well. For example, we can construct the An E An as empirical distribution of samples generaled from X ji,, Le for each N, 3Bx N, s.t.VB N, BN. W, PB) </22.
Since P3 is the minimizer of Wp Wp, PB) μD AB, Ve, EN, BN EN, B BN, s.t.:
Wp(,, I Wp(u, PB)
(35)
S Wpl,, (Tx) + W (B)
<.
Therefore, (1) + 0.
By lemma 9
we have uB 4 p.
C.8 Proof of Theorem 5.6
The proof is a direct result of the Theorem of Défossez et al. (020)_ order to apply the theorem, There are three conditions in Section 2.3 of Défossez et al. (2020). Slightly revised to more suitable our situation, the three conditions are:
(a) The function S)) bounded below, i.e.33C, R, s.t.
vy VD Sx(V, C..
(b) The gradient of Sx()),)) D) is uniformly almost surely bounded, i.e.,RR Ve (Ve is used here simplify the final bounds), s.t.
DDS D) SR- Ve.
(c) The gradient of Sx(2, D) is L-Lipschitz continuous with respect to the 1o-norm, ie., 3L > 0, s.t.
Sx01, D2) SA,, s L[[D; - D2ll2-
Since VDSA(D,D)- = VMSX(), and notice that VpM is just derivative of point-wise 1o-msrms, it is sufficient to verify that SX), D) satisfies these properties.
According to the analytic form for VMSx provided in Theorem reviewing the required properties by definitions, it is trivial to verify the three conditions.
Hence, we have, 30, s.t.:
SA D(T)) s CR + RL + 22)
log N VN
Let C, = CRP RL L2), and we proved the theorem.
C.9 Proof of Proposition 5.7
The proof is a direct consequence of the triangle inequality of the Sinkhorn loss (Luise al. 2018 1f
S2 (D())
D(l-1)) < 8, then
(S, 5.0)DD0-111 Si(()) D(-11) (DDD)) < 8.
40


Reviews:
[{"review_id": "QyOLORTPJc", "reviewer": "Reviewer_2HMD", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\nThe paper is well-written overall with clear explainations in issues faced under Sinkhorn stability. Further the convergence analysis and stability proofs strengthen the arguments for the paper's main proposal against solving instability issues of traditional Sinkhorn algorithm settings. As part of experiments, the authors showcase better coreset sample quality selection than other sampling approaches (e.g. Fig 4) as well as achieving better efficiency.\n\nFurther the authors showcases new theoretical results showcasing stability analysis for Sinkhorn loss along with numerical analysis to support the same.\n\n**Weakness**\n\n- Not clear about the essence of Proposition 2.1. Why is this proposition important?\n- Eq 4 arguments are not clear. How is Eq 3 can be directly refined to Eq 4? Since Sinkhorn loss has 3 arguments in previous definition.\n- I might have missed the definition of BFGS from the paper. Could authors elaborate on it more in their writeup? In terms about the description of the algorithm.\n- Can you elaborate a bit more on the Sinkhorn log algorithm. Is it the entire u is transformed to log scale or just $M_e$.\n- I am not sure in terms of baselines, if you could also have tried out submodular based sample selection, as they also provide better guarantees", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes an important problem towards coreset selection particularly Wasserstein coresets using Sinkhorn loss. To adress stability issues in Sinkhorn algorithms, authors propose stable algorithms and analyse via experimentation. The paper extensively studies the stabiliy analysis of L-BFGS from 1989 paper which is utilised in several recent works to solve for alternative to forward pass. Further through their experiments, the authors showcase how their approach performs better than other existing methods in terms of sample selection quality and efficiency.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- Fig 5. Runtime analysis if it is possible to overlay Coreset Runtime analysis w.r.t Full data training for a better comparison.\n- Improve Figure descriptions: For example Fig 2: pls include the definition of $\\lambda$ and high level observations from the plot in the figure description, show that at a glance it becomes clear what the figure tries to convey.  \n- In the same page , description of Fig 2. its written $\\eta$ instead of $\\lambda$.", "broader_impact_concerns": "There are no broader impact concerns in this work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "4OMUtEBaRw", "reviewer": "Reviewer_Gxub", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\nNovel analytical derivation of Sinkhorn loss derivatives\n\nWeaknesses:\n- Limited comparison with existing coreset algorithms beyond SCCP and random sampling\n- Insufficient justification of the task-agnostic nature of the coresets.\n- Lack of clarity of presentation.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper is on selecting a representative subsample, also called a coreset, that summarizes a given dataset. The authors use Wasserstien metric to formalize the notion of a representative subset. Since computing Wasserstein distance between two measures is expensive, this paper works with an approximation of Wasserstein distance, called the Sinkhorn loss. In Algorithm 1, the authors give a gradient based algorithm to construct what they call the Wasserstein coreset via Sinkhorn loss (WCSL). The authors then give analytic expressions for the gradient of the Sinkhorn loss function, state theorems on stability and convergence properties of L-BFGS/ADAM  algorithm for minimizing this loss function. The paper contains various experiments on synthetic and real data to prove the validity of the algorithm's output.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I have several concerns and I cannot provide a reocmmendation for acceptance without them being addressed.\n\n(1) How does the coreset algorithm compare with other coreset algorithms? The authors early on make a distinction of task specific vs. task agnostic methods, please define clearly what is a task agnostic coreset. What about baselines like centroids of k-means, Kernel Thinning (https://arxiv.org/abs/2105.05842),...\n\n(2) The authors claim that their coresets are task agnostic. Why ?\n\n(3) It is not clear to me how line 7 in Algorithm 1 type checks. It essentially reads:\n set_new <- set_old + gradient_term . \nPlease explain what this means. What is '+' for a set, if the gradient_term is a set.\n\n(4) The authors talk of forward and backward passes of Sinkhorn loss function. Please define what this means to someone who has seen forward and backward passes only in the context of neural networks.\n\n(5) How does gradient descent compare with L-BFGS for the regularized OT problem. What is known in terms of stability and convergence of gradient descent?\n\n(6) In Theorem 5.2, what is epsilon? It is not defined before it appears in the guarantees. \n\n(7) In theorem 5.4, what is B fixed or increasing? In the explanation that follows \"*when the number of points n grows large\", what is 'n'? *\n\n(8) Algorithm 1 has only gradients of Sinkhorn loss. What is the implication of the results on L-BFGS/Adam and duality on Algorithm 1?\n\nThe current exposition has a lot of gaps and several things are not clear. I recommend a thorough revision.", "broader_impact_concerns": "I don't have any concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "9clixZVUAQ", "reviewer": "Reviewer_y7hu", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strong Aspects:\n\n1. The problem of coreset selection is important, and Sinkhorn's algorithm,\n   while powerful, has the issue of numerical instability. The paper tries to\n   address this significant instability issue, which is crucial for making the\n   method reliable in practical applications.\n2. The derivation of an analytic expression for the Sinkhorn loss derivative\n   enables its smooth application in practical machine learning frameworks.\n3. The proposed WCSL is very general, making it applicable across different\n   machine learning tasks without the need for task-specific modifications. \n4. The numerical experiments involving both simulation studies and real-world\n   applications are comprehensive. \n\n### Weaker Elements:\n\n1. The presentation of the paper is dense in mathematical notations, which is\n   acceptable for a theoretical investigation. However, for some key concepts,\n   such as Wasserstein distance and Sinkhorn Loss, the notations are not always\n   clear. The paper could benefit from more consistent and clear presentation of\n   notations.\n   - For example, in the introduction, what exactly are $\\mathcal{X}$ and $X$?\n     While it may be easy to infer their meanings, an explicit statement would\n     be helpful.\n   - The term \"p-Wasserstein distance\" is introduced on page 2, but the meaning\n     of $p$ is not clearly explained. Similarly, the notation $W_p(\\cdot,\n     \\cdot)$ is not well-defined at this point.\n   - What is the variable η referring to in \"we globally define\n     $\\eta=\\lambda^{-1}$ for simplicity\"?\n   - As noted below Equation (1), the Wasserstein distance is presented using\n     varied notations, which may make it easier for the authors, but it is not\n     reader-friendly. This could confuse readers who are less familiar with the\n     topic. More intuitive explanations would make the paper more accessible.\n  \n2. There is limited discussion on alternative methods. The proposed method,\n   though general, is closely related to existing approaches in the literature,\n   such as leverage sampling (1) and optimal subsampling (2). A comparison with,\n   or mention of these methods, would broaden the scope and add value to the\n   paper.\n\n   (1) Drineas, P., et al. *Fast approximation of matrix coherence and\n   statistical leverage.* JMLR, 2012.\n\n   (2) Wang, H., Zhu, R., & Ma, P. *Optimal subsampling for large sample\n   logistic regression.* JASA, 2018.\n\n3. There is no time complexity analysis, which hinders understanding of the\n   computational efficiency in more detail.\n\n4. While the paper focuses on improving computational efficiency, it lacks\n   exploration of scalability to extremely large datasets. Expanding on the\n   scalability in larger datasets would improve the paper.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper studies generating Wasserstein coresets using the Sinkhorn loss, which\nis a technique inspired by optimal transport theory. The authors try to address\nthe numerical instability by applying stochastic gradient decent. The paper\nrigorously analyzes the stability of the forward and backward computations from\na theoretical perspective, and provide extensive experiments, including\nreal-world applications with image data, to demonstrate the performance of the\nproposed methods.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. Improve the presentations, especially on the consistence and definitions of\n   notations. See the detailed comments in the \"Weaker Elements\".\n2. Discuss and/or compare with more benchmarks.\n3. Provide some time complexity analysis. This would help readers to understand\n   the possibility of scaling the method to larger datasets.\n4. If possible, consider a case of larger data, either in simulation or real\n   data.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper is a solid contribution, and definitely TMLR-worthy after the revisions made during the review process.

There are three minor fixes that need to be made to the manuscript prior to publication.
1. Figure 9 does not have kernel thinning results, but kernel thinning is in the figure legends. In your response to Reviewer Gxub you state that "For large image datasets MNIST and FashionMNIST, the kernel thinning method proved very slow; thus, we did not include the results for kernel thinning on image datasets." Please fix the legends in Figure 9 and add this explanation to the paper.
2. Léon Bottou's name did not typeset correctly in the reference for the Wasserstein GAN paper.
3. Benoît Kloeckner's name did not typeset correctly in the reference for "Approximation by finitely supported measures"
Decision: Accept (poster)