Efficient, Stable, and Analytic Differentiation of the Sinkhorn Loss
Authors: Yixuan Qiu
Haoyun Yin
Xiao Wang
Venue: Submitted to ICLR 2023 (2023)
Keywords: Optimal transport
Wasserstein distance
Sinkhorn loss
differentiable optimization
deep generative model

Abstract: Optimal transport and the Wasserstein distance have become indispensable building blocks of modern deep generative models, but their computional costs greatly prohibit their applications in statistical machine learning models. Recently, the Sinkhorn loss, as an approximation to the Wasserstein distance, has gained massive popularity, and much work has been done for its theoretical properties. To embed the Sinkhorn loss into gradient-based learning frameworks, efficient algorithms for both the forward and backward passes of the Sinkhorn loss are required. In this article, we first demonstrate issues of the widely-used Sinkhorn's algorithm, and show that the L-BFGS algorithm is a potentially better candidate for the forward pass. Then we derive an analytic form of the derivative of the Sinkhorn loss with respect to the input cost matrix, which results in a very efficient backward algorithm. We rigorously analyze the convergence and stability properties of the advocated algorithms, and use various numerical experiments to validate the superior performance of the proposed methods.

EFFICIENT, STABLE, AND ANALYTIC DIFFERENTIA- TION OF THE SINKHORN Loss
Anonymout authors Paper under double-blind review
ABSTRACT
Optimal transport and the Wasserstein distance have become indispensable build- ing blocks of modern deep generative models, but their computational costs greatly prohibit their applications in statistical machine learning models. Recently, the Sinkhorn loss, as an approximation to the Wasserstein distance, has gained massive popularity, and much work has been done for its theoretical properties. To embed the Sinkhorn loss into gradient- gradient-based learning frameworks, efficient algorithms for both the forward and backward passes of the Sinkhorn loss are re- quired. In this article, we first demonstrate issues of the widely-used Sinkhorn's s algorithm, and show that the L-BFGS algorithm is a potentially better candidate for the forward pass. Then Sinkhorn loss with respect to the input cost matrix, which results in an efficient backward algorithm. We rigorously analyze the convergence and stability proper- ties of the advocated algorithms, and use various numerical experiments to vali- date the performance of the proposed methods.
derive an analytic form of the derivative of the
INTRODUCTION
Optimal transport (OT, Villani 2009) is a powerful tool to characterize the transformation of proba- bility distributions, and has become an indispensable building block of generative modeling. At the core of OT is the Wasserstein distance, which measures the difference between two distributions. For example, the Wasserstein generative adversarial network (WGAN, Arjovsky et al., (017)) uses the 1-Wasserstein distance as the loss function minimize the difference between the data distribution and the model distribution, and a huge number of related works emerge afterwards.
Despite the various appealing theoretical properties, one major barrier for the wide applications of OT is the difficulty in computing the Wasserstein distance. For two discrete distributions, OT solves a linear programming problem of num 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 log n.) (Pele & Werman 2009). which quickly becomes formidable gets large, except for some special cases (Peyré et al. 2019).
To resolve this issue, many approximate solutions OT have been proposed, among which the Sinkhorn loss has gained massive popularity (Cuturi, 2013). The Sinkhorn loss can be viewed as an entropic-regularized Wasserstein distance, which adds a smooth penalty term the original ob- jective function of OT. The Sinkhorn loss is attractive as its optimization problem can be efficiently solved, least in exact arithmetics, via Sinkhorn's 's algorithm (Sinkhorn) 1964 Sinkhorn & Knopp, 1967]. which merely involves matrix-vector multiplications and some minor operations. There- fore, especially suited to modern computing hardware such the graphics processing units (GPUs). Recent theoretical results show that Sinkhorn's 's algorithm has a computational complexity of O to output an e-approximation to the unregularized OT (Dvurechensky al. 2018).
Many existing works on the Sinkhorn loss focus on its theoretical properties, for example Mena this article, we are mostly concerned with
& Niles-Weed (2019) and [enevay et al. (2019). the computational aspect. Since modern deep generative models mostly rely on the gradient-based learning framework, crucial to use the Sinkhorn loss with differentiation support. One simple and natural method to enable Sinkhorn loss in back- propagation is to unroll Sinkhomn' algorithm, adding every iteration to the auto-differentiation computing graph (Genevay et al. 2018)) Cuturri et al. 2019). However, this approach is typically costly when the number of iterations are large.

---

Instead, in this article we have derived analytic expression for the derivative of Sinkhorn loss based on quantities computed from the forward pass, which greatly simplifies the back prooparopggttion of the Sinkhorn loss.
More importantly, one critical pain point of the Sinkhorn loss, though typically ignored in theoretical studies, is that Sinkhorn's algorithm is numerically unstable (Peyré et al. (019). We show in numer- ical experiments that even for very simple settings, Sinkhorn's 's algorithm can quickly lose precision. Various stabilized versions of Sinkhorn's 's algorithm, though showing better stability, still suffer from slow convergence in these cases. In this article, we have rigorously analyzed the solution to the Sinkhorn optimization problem, and have designed both forward and backward algorithms that are provably efficient and stable. The main contribution of this article is as follows:
We have derived an analytic expression for the derivative of the Sinkhorn loss, which can be efficiently computed in back-propagation.
We have rigorously analyzed the advocated forward and backward algorithms for the Sinkhorn loss, and show that they have desirable efficiency and stability properties.
We have implemented the Sinkhorn loss as an auto-differentiable function in the PyTorch and JAX frameworks, using the analytic derivative obtained in this article.
The code to reproduce the results in this article is available at https: htmps://100 //1drv ms /u/s! ArsORq8a2 424 WWooF W
2 THE (SHARP) SINKHORN LOSS AS APPROXIMATE OT
Throughout this article we focus on discrete OT problems. Denote by A, = {w E R : w w = }} the n-dimensional probability simplex, and let μ = Σi= 0,62, and v = Σド= by be two discrete probability measures supported on data points {-:)==1 and {u, [y respectively, where Am. and 82 is the Dirac at position Define [1 b) = (T R .Tmm a. 1",, 1, b). and let M E Raxm be cost matrix with entries M,j. i = ,....,n,7 n, I,...,
(a1,... an)T E Ans b = (b1, 6m)
,....,m.. m. Without loss of generality we assume that n 2 m, as their roles can be exchanged. Then OT can be characterized by the following optimization problem,
W(M,a,b) = min (P, Pell(a
(1)
where (A, B) B). An optimal solution to D.. denoted 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 (M, a, b) = [P*,M] is then called the Wasserstein distance between a and v if the cost matrix M satisfies some suitable conditions (Proposition 2.2 of Peyré et al. 2019).
As is introduced in Section solving the optimization problem can be difficult even for moderate and m. One approach to regularizing the optimization problem is to add an entropic penalty term the objective function, leading to the antropic-regularizee OT problem (Cuturi) 2013):
3x(M, = min TETI(a,))
Sx(T) == min TETI(a,))
(T, X-T)(T)
(2)
where ((T)= ET= (1 - log Tıy) is the entropy term. The new objective function Sx(T) X- -troonly convex on l1(a, b). so 2 has unique global solution, denoted as Tx. In this article, Tt is referred to as the Sinkhorn transport plan. The peropic-regularizer Wasserstein distance, also known as the Sinkhorn distance or Sinkhorn loss in the literature (Cuturi 2013), then defined as
Sx(M a, b) (T),M).
To simplify the notation, we omit the subscript À in T5 hereafter when no confusion is caused. It is worth noting that in the literature, Ss and S5 are sometimes referred to the sharp and regularized Sinkhorn loss, respectively. The following proposition from Luise et al. (2018) suggests that Sx achieves a faster rate at approximating the Wasserstein distance than Ss. Due to this reason, in this article we foeus on the sharp version, and simply call S2 the Sinkhorn loss for brevity.
Proposition (Luise et al. 2018). There exist constants C1, C2 > 0 such that for any X > O, |SxM,, a, b) - W(M a, b)l < Cee and (S(M,a. a, b) W(M, a, b) & C2/X X. The constants C₁ and C2 are independent X, and depend on μ and
2

---

3 DIFFERENTIATION OF THE SINKHORN LOSS
To use the Sinkhorn loss in deep neural networks or other machine learning tasks, it is also cru- cial obtain the derivative of Sx (M, a, b) with respect its input parameters. Differentiating the Sinkhorn loss typically involves two stages, the forward and backward passes. In the forward pass, the Sinkhorn loss or the transport plan computed using some optimization algorithm, and in the backward pass the derivative is computed, using either an analytic expression or the automatic differentiation technique. In this section we analyze both passes in details.
Throughout this article we use the following notations, For y, E R. x A means mi, v). For vector v = (u,,..., ex let wr = (vi ..... Ux "T ū (UI,..., U 1)-1)², and denote by diag(v) the diagonal matrix formed by v. Let u = (ui...., 4) be another vector, and denote by u the 1 k matrix with entries (u, + vÿ). For matrix A (ary) (A...... Ax) with column vectors AI...., Ag. let À (AI..... As-1), and ex[A] be the matrix with entries eaar, The symbol @ denotes the elementwise multiplication operator between matrices vectors. ||- and || || stand for the Euclidean norm for vectors and Frobenius norm for matrices, respectively. Finally, we globally define 1) A for simplicity.
3.1 ISSUES SINKHORN'S S ALGORITIIM
In the existing literature, one commonly-used method for the forward pass of the Sinkhorn loss Sinkhorn's 's algorithm (Sinkhorn) 1964, [Sinkhorn & Knopp 1967]. Unlike the original linear programming problem ( the solution to the Sinkhorn problem has special structure. Cuturi 2013) shows that the optimal solution T can be expressed as
T diag(u" (aaagg"""Meecii ). diag(v")
for some vectors u" and v², where Me = (e (e)) Sinkhorn's algorithm starts from an initial vector (0)) IR,, and generates iterates (kk) 17 and ((k)) ITT follows:
(k+1) + a 2 [M (b+1) + b [M ((tt)+)))11
(4)
It can be proved that (k)) + u² and (k) + v², and then the Sinkhorn transport plan T" can be recovered by (3).
Sinkhorn's algorithm is very efficient, as it only involves matrix-vector multiplication and other inexpensive operations. However, one major issue of Sinkhorn's 's algorithm that the entries of M, = (e-AMy may easily underflow when À is large, making some elements of the vectors M, ((k) and M (++11 close to zero. As a result, some components of (1+11 and ((++1) would overflow. Therefore, Sinkhorn's algorithm in its original form is unstable, and in practice the iterations (4p are typically carried out in the logarithmic scale, which we call the Sinkhorn-log algorithm for simplicity. Besides, there some other works also attempting improve the numerical stability of Sinkhorn's 's algorithm (Schmitzer] 2019 Cuturi et al.) 2022).
Despite the advancements of Sinkhorn's '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 iterations. Below we show a motivating example to highlight this issue. Consider triplet (M, a, b) for the Sinkhorn problem, and fix the regularization parameter 77 to be 0.001, where the detailed setting is provided in Appendix B.1 The true T matrix visualized in Figure along with the solutions given by various widely- algorithms, including Sinkhorn's 's algorithm, Sinkhorn-log, the stabilized scaling algorithm (Stabilized, Algorithm 2 of Schmitzer 2019), and the Greenkhorn algorithm (Altschuler et al., 2017, Lin et al., 2002). The maximum number of iterations is 10000 for Greenkhorn and 1000 for other algorithms.
In Figure I 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 in the plot. In contrast, the stable algorithms Sinkhorn-log and Stabilized greatly improve them. Sinkhorn's 's algo- rithm and Sinkhorn-log the need for numerically stable algorithms. However, Sinkhorn-log and Stabilized still have visible inconsistencies with the truth even after 1000 iterations.
equivalent in exact arithmetics, so their numerical differences highlight
3

---

Truth
Sinkhorn
Sinkhom-log
Stabilized
Greenkhom
L-BFSS
Figure 1: Visualization of Sinkhorn plans computed by different algorithms.
3.2 THE ADVOCATED ALTERNATIVE FOR FORWARD PASS
To this end, we advocate an alternative scheme to solve the optimal plan T*, and we show both theoretically and empirically that this method enjoys great efficiency and stability. Consider the dual problem of ( which has the following form (Proposition 4.4 of Peyré et al. 2019):
max C(a, 3) = max a²a + 3 0
Σ 2m
-aM--a a E R" "RR,BER" (5)
Let a' = (ai..... aC,) and 3 = (Bf..... );,) be one optimal solution to and then the Sinkhorn transport plan T can be recovered as T ex|a"088* Remarkably. is equivalent to an unconstrained convex optimization problem, so a simple gradient ascent method suffices to find its optimal solution. But practice, quasi-Newton methods such as the limited memory BFGS method (L-BFGS, Liu & Nocedal, 1989) can significantly accelerate the convergence. Using L- BFGS to solve (5) is known practice (Cuturi & Peyré 2018 Flamary al.) (02I), but little is known about its stability in solving the regularized OT problem. We first briefly describe the algorithm below, and in Section 4 we rigorously prove that L-BFGS converges fast and generates
stable iterates.
is worth noting we can reduce the number of variables to be optimized (5) based on the follow-
ing two findings. First, pointed out by Cuturi et al. (2019), the variables (a, B) have one redundant degree of freedom: (a" (8*) one solution to (5), then so is (a* +cl,, 8" cllmm for any c. There- fore, we globally set 8m = 0 without loss of generality. Second, let a (B) = arg max. L(a, [()a((()))))))))) ())),B) B). Then only need to minimize f(3) with (m - 1) free variables to get 3 and a can be recovered as a = (89).
for given 3 (8T ) T (3T, (BT,0)T, and define f(3)
the appendix show that a² (3), f(3), and V Vzf have simple closed-form expressions:
= -
a ( = 17 17 - log
e M0a- - Ma))
(6)
Vzf T( b.
T(0) ca (8) M].
With f(B) and Vzf. the L-BFGS algorithm can be readily used to minimize f(ß) and obtain 3*. Each gradient evaluation requires ((mn) exponentiation operations, which is comparable to Sinkhorn-log. Although exponentiation more expensive than matrix-vector multiplication in Sinkhorn's is algorithm, this extra cost can be greatly remedied by modern hardware such GPUs. On the other hand, we would show in Section g that L-BFGS has strong guarantee on numerical stability, which is critical for many scientific computing problems.
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 eval- uations to 1000, and hence comparable to other methods. Clearly, the L-BFGS solution is visually identical to the ground truth. To study the difference between L-BFGS and Sinkhorn's 's algorithm in more depth, we compute the objective function value the dual problem 6 at each iteration for both Sinkhorn-log and L-BFGS. The results visualized in Figure 2 with three different 17 values, 7 = 0.1, 0.01, 0.001. Since each L-BFGS iteration may involve more than one gradient evaluation, for L-BFGS plot the values against both the outer iteration and gradient evaluation counts.
For the motivating example Section 3.1
4

---

--1
x-1=
0.000
6.1-
0.0
Function
-05
-03-
Sinkhorn-log
-15
-15
L-FGGS Gradient Eval.
L-BFGS iteration
-1.3
-25
-20
38
Iteration Count
Figure 2: Comparing the convergence speed of Sinkhorn-log and L-BFGS.
Figure 2laree2givve gives clear clue to the issue of Sinkhorn-log: it has a surprisingly slow convergence speed compared to L-BFGS when v small. Theoretically, Sinkhorn algorithms with sufficient iterations, but in practice, moderate limit on computational budget may prevent them from generating accurate results. To this end, L-BFGS appears to be a better candidate when needs small 7 for better approximation the OT problem. In Appendix [B22 we have designed more experiments to study the forward pass stability and accuracy of different algorithms.
eventually converge
3.3 THE ANALYTIC BACKWARD PASS
For the backward pass, one commonly-used method is unrolled Sinkhorn's algorithm, which is based the fact that Sinkhorn's forward pass algorithm differentiable with respect a, b, and M. Therefore, one can use automatic differentiation software to compute the corresponding deriva- tives in the backward pass. This method is used in Genevay al. (2018) for learning generative models with the Sinkhorn loss, but in practice it may be extremely slow the forward pass takes large number of iterations. To avoid the excessive ferentiation methods have been developed for the Sinkhorn loss (Feydy et al. (2199) Campbell et al., 2020 Xie et al. 2020 Eisenberger al. 2022). but they still not provide the most straightforward way to compute the gradient.
of unrolled algorithms, various implicit dif-
To this end, we advocate the use of analytic derivatives of the Sinkhorn loss, which solves the optimization problem (2) in the forward pass, and use the optimal dual variables (a", 8*) for the backward pass. The analytic form for Sx a, b) has been studied in Luise et al, (2018), and to our best knowledge, few result has been presented for V SXM,, a, b). Our first main theorem, given Theorem fills this gap and derives the analytic form for V SX
Theorem 1. For a fixed >
V SAr,A(M, T A(sa &8 M)OT.
(5,, Mr (M ④ 1 Te = (M @ T T 3, = D- - -T + (4) and D = diag(b) - - diag(a "y A., b E Am, and M.
where T ex/a" + 3r MI, 8n a- ② (1r 705.).
In addition, is positive definite, and hence invertible, for any
Assuming 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 positive definite linear system.
CONVERGENCE AND STABILITY ANALYSIS
In Sections 3.2 and 3.3 we have introduced the advocated forward and backward algorithms for the Sinkhorn loss, respectively. In this section, focus the theoretical properties of these algo- rithms, and show that they enjoy provable efficiency and stability. As a first step, we consider the target of the optimization problem (6). and show that ((3) has well-behaved minimizer, which does not underflow or overflow.
Theorem 2. Denote by f" the minimum value of f(3) and 8 an optimal solution, and let a² o (3*). Then f > -c, 8 is unique, ||a"|| 2o, and ||,"1|| co. In particular, let I
5

---

arg m Tim 0max max; dg and log(n/bm). Then Lo. Las og S Uas and
Lg, S 8j s U3, I UU, for all ,,.,,, n and j = 1, ....,, m, where
Ua, Mm + +1 log(a; A bm),
Up, M15- M +1/ [oog(a A bb) + c]
Lo. 17 log(a: (a(a:/m) /m/m- J (Up,- Mij),
LØ, n log(bb/n) max (Ua, -Mis);
Up max - 3 + 1) (log( A bj) + +cl, La. = n loo((o,//)) max (U MMi)).
Theorem 2 shows that the optimal dual variables y" are bounded, and more importantly, the bounds are roughly at the scale of the cost matrix entries Me and the log-weights log(a;) and log(b;). Therefore, least the target of the optimization problem is well-behaved. Moreover, one useful application of Theorem is to impose a box constraint on the variables, adding further stability to the optimization algorithm.
After verifying the properties of the optimal solution, a more interesting and critical problem is to seek a stable algorithm to approach the optimal solution. Indeed, in Theorem 3 we prove that the L-BFGS algorithm for minimizing (6) is one such method.
Theorem 3. Let {300] be a sequence of iterates generated by the L-BFGS algorithm starting from a fixed initial value 300), and define g(*) = (3(K),,0) ,O)T, T(k) = ex(a" ex(a)(3(() @ 3(0) M). and f(h) = =((3t)). Then there exist constants 0 I r < and C1, C2 > O such that for each > 0:
(a) (16) -r <+ - = =(k).
(Linear convergence for the objective value)
(b) [13(0)) - 13- s C2e(k)
(Linear convergence for the iterates)
(c) T(a) 1m=a, a. |1 11= b112 < C2e(1).
(Exponential decay of the gradient)
(d) Tilh) < min bi + (2s( J. 1 smsm-- m-
(T(k) does not overflow)
(e) max; T >a./m, max; T (k) (b, (b,--/2eee (2e(ki) nn 15 ssssm- (T(K) does not underflow)
The explicit expressions for the constants C1, O2,r r are given in Appendix A
Theorem reveals some important information. First, both the objective function value and the iter- ates have linear convergence speed, so the forward pass using L-BFGS takes O(lo [[log(1/)) iterations to obtain e-optimal solution. Second, the marginal error for μ., measured by (|] 1.) all. is exactly zero due to the partial optimization equal to the gradient norm, is also bounded at any iteration, and decays exponentially fast to zero This validates the numerical stability of the L-BFGS algorithm. Third, the estimated transport plan any iteration k (k)), is also bounded and stable. This result can be compared with the formulation in By: it is not hard to find that u", v", and Mc, when computed individually, can all be unstable due to the exponentiation operations, especially when 1/ is small. In contrast, T and (()), thanks to the results in Theorem 3 do not suffer from this issue.
The other marginal error [ b11, which is
We emphasize that Theorem 3 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 impor- tantly, 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 (MS.,M a, b) as in Theorem I can also be computed in a numerically stable way. Let S be the k-step approximation V SM = V SM (M, a, b) using the L-BFGS algorithm, i.e., replacing every T in V by (k)) Then we show that the error on gradient also decays exponentially fast.
Theorem 4. Using the symbols defined in Theorems I and the smallest eigenvalue of D. Assume that for some Kip
let a 1/ (D), where a (D)
min 1, (6a D
(Ka) <C
4X

---

and then for every k 2 ko. V V S Cyve(k) C -ft ,k/2, where the ex-
plicit expression for Cs given in Appendix A Ki always exists g (8)) decays to zero exponentially fast as ensured by Theorem 3ta).
APPLICATION: SINKHORN GENERATIVE MODELING
The Sinkhorn loss useful in unsupervised learning tasks that attempt match two distributions pe and p", where p" stands for the data distribution, and is the model distribution. If pa and p can be represented or approximated by two discrete measures 6e and v, respectively, then one would wish to minimize the Sinkhorn loss (Mo 0g, bo) between 140 and v, where the cost matrix M and weights a, b may depend on learnable parameters In gradient-based learning frameworks, the key step of seeking the optimal parameter 0 that minimizes (MM,, as, be) is compute the gradient VaSx(Ma, ag, by). which further reduces to evaluating Va.bSx a, b) and V SMMMMM a, b).
Luise et al. (0018) assumes that M is fixed, and studies the gradients Va,5Sx(M b). However, in many generative models is more important to derive SA a, b), as the weights a and typically fixed, whereas the cost matrix M is computed from the output of some parameterized layers. Consider a data set ,,...., Xn, that follows the distribution p'. X; E Rº, and our target to fit deep neural network go R" R such that (o approximately follows p", where Z ~ (00 I,). To this end, we first generate random variates Zz.... Zm ~ N(0, I,), and let Y - (2,). The two weight vectors are simply taken to be
$1, b 1n mm and the cost matrix between X, and Y, is given by M3 E Raxm where (Ma)ss |X, g0 (2)z,, 13 Then learn the network parameter 0 by minimizing the Sinkhorn loss {(0) = SM n 1 m 1m).
We show such applications in Section
6 RELATED WORK
In T we list some related work on the differentiation of the Sinkhorn loss, and provide a brief summary of the contribution and limitation of each work. The target of our work is closest to that of Genevay et al. (2018), i.e., to differentiate the Sinkhorn loss with respect to the input cost matrix. However, they use unrolled Sinkhorn's 's algorithm, so there is no analytic form for V Sx. Our work mostly motivated by Luise et al. (20T8), but they consider the derivative with respect to the weights instead of the cost matrix. Similarly, Cuturi et al. (2019) and Cuturi et al. (2000) consider gradients on weights and data points, but not the general cost matrix. Campbell et al. (2020); Xie et al. 22200;; Eisenberger et al.) (2022) all consider the derivative of the transport plan T using implicit differentiation. Although this is more general problem than computing V as SA, it loses the special structure of the Sinkhorn loss. As result, the compact matrix form of the derivative presented in Theorem is unique. Furthermore, the storage requirement for our result is (nm)), whereas some existing works need to store much larger matrices. Finally, very few work has rigorously analyzed the stability of the algorithm, in both the forward and backward passes.
7 NUMERICAL EXPERIMENTS
7.1 RUNNING TIME OF FORWARD AND BACKWARD PASSES
In our subsequent experiments, we compare three algorithms for differentiating the Sinkhorn loss: the proposed method ("Analytic") using L-BFGS for the forward pass and the analytic derivative for backward pass, the implicit differentiation method ("Implicit"), and the unrolled Sinkhorn algorithm ("Unroll"". Both Implicit and Unroll use Sinkhorn-log for the forward pass, and are implemented in the OTT-JAX library (Cuturi et al. 2022).
We simulate cost matrices of different dimensions, and compare the forward time and total time of different algorithms with both 7) 0.1 and y = 0.01. The detailed experiment setting and results are given in Appendix B.3 Table and Table B The main conclusion that under the same accuracy level, the forward pass of Analytic is typically faster than those of the other two algorithms, and the backward time is significantly faster. Overall. the proposed analytic differentiation method demonstrates visible advantages computational efficiency.
7

---

5 and hidden Table A brief summary of existing literature on differentiation of Sinkhorn loss.
Derivative w.r.t the cost matrix
Minimum storage requirement
Analytic solution
Compact matrix form
Stability analysis
Reference
Target
Genevay et al.
2018
Sx
X
X
Luise et al.
2018
Va
/
X
Cuturi et al.
2019
Va
x
Cuturi ct al. 2020
Vo p
/
X
Camphell et al.
2020)
V MMT T
/
X
V
X
Xie et al 2020)
V MT*
/
/
X
/
X
Eisenherger al.
2022
V T
/
X
This article
SA
/
7.2 GENERATIVE MODELS ON TOY DATA SETS
In this section we apply the Sinkhorn loss to generative modeling, and test the accuracy and effi ciency of the proposed algorithm. Following the methods in Section we consider a toy data set with n 1000 and p (shown in Figure that go approximately pushes forward Z1. Zm (0, I,) to the data distribution. In our ex- periments, go is a fully-connected ReLU neural network with input dimension dimensions 64-128-64. The number of latent data points is m 1000.
3 and we attempt to learn neural network ge such
In the first setting, we intentionally keep both the observed data and the latent points Zi.... Zm fixed, so that the optimization of is deterministic process without randomness, and the optimiza- tion variable obtained from each forward pass used as the warm start for the next training iteration. At every ten iterations, we compute the 2-Wasserstein distance between the observed data and the pushforward points {ge( (. and we train ge for 2000 iterations using the Adam optimizer with a learning rate of 0.001. This setting, though not common in generative modeling, helps to monitor the computation of gradients without the impact of random sampling. Moreover, the Wasserstein distance has achievable lower bound of zero if go is sufficiently expressive, accuracy and efficiency of gradient computation by plotting the metric against running time.
study the
The comparison results for the three algorithms shown in the second and third plots of Figure 3 from which we have the following findings. First, when y = 0.5, the 2-Wasserstein distance will increase after certain number of iterations, indicating that the Sinkhorn loss is not aceu- rate approximation to the Wasserstein distance when 17 is large. Second, it is clear that the Unroll method has the slowest computation, as it needs to unroll potentially large computational graph in automatic differentiation. Third, the proposed analytic differentiation shows visible advantages on computational efficiency, thanks to the closed-fome forsed-fomm expression for the derivative.
Finally, to examine the performance of differentiation methods in genuine generative modeling setting, we use a regular training scheme that randomly samples the observed data and latent points at each iteration. Due to the first finding above, we choose the smaller y value for the Sinkhorn loss, and use a mini-batch of size n = m = 256 for training. We run 5000 iterations in total, and the metric curves are shown in the last plot of Figure 3 It can be found that the performance of the three algorithms are similar to the fixed-2 2 case: all three methods properly reduce the Wasser- stein distance over time, but the proposed algorithm uses less time to accomplish the computation. Additional experiments on simulated data sets are given in Appendix B.4
7.3 DEEP GENERATIVE MODELS
Finally, use the Sinkhorn loss to train larger and deeper generative models the MNIST (Le- Cun et al. 1998) and Fashion-MNIST (Xiao et al., 2017) data sets. In this experiment, we do not pursue training the best possible generative model; instead, primarily validate our claims on the
8

---

Thes dall ports
Food Zi 2m 0.5
Tood 2, gm0.
Samples 2, gm
impict
Impkt Unal
inst et
lime Gecunzal
Expoed Time Geunnzs)
Eippoll Tini (500.1dd))
Figure 3: From left to right: the true observed data; the measured Wasserstein distance over time in the fixed Z setting with 77 0.5; the similar plot for 7 = 0.1; the random Z setting with nj 0.1.
VIII D303
NNBSTT Data
9
G
2
2
2
0
6
6
7
Expted Time (eecont))
Elppse me Coccendd
(a)
(c)
Dvction NMEST Dabs
Pection MNIST DATS
The Dac ouan Diaggnng
Expeed Thme (eccone)
Dispeed Time concender
(d)
Figure 4: (a)() Sinkhorn divergence between randomly sampled test data and model-generated images over iterations. (b)(e) Similar (a)(d) but using 2- 2-Wasserstein distance. (c)(f) Randomly generated images from the models trained using the analytic differentiation.
efficiency of the proposed analytic differentiation method. The architecture of neural networks and details of the training process are given in Appendix [B.5] During the training process, we com- pute the Sinkhorn divergence (Feydy al.) 2019) and 2-Wasserstein distance between randomly sampled test data and model-generated images over iterations, and Figure + show the results of Analytic, Implicit, and Unroll algorithms. The conclusion is similar: the proposed analytic differ- entiation method achieves the same learning performance as existing methods, but is able to reduce the overall computing time. This once again demonstrates the value of the proposed method.
CONCLUSION
In this article we study the differentiation of the Sinkhorn loss with respect to its cost matrix, and have derived an analytic form of the derivative, which makes the backward pass of the differentiation easy to implement. Moreover, we study the numerical stability of the forward pass, and rigorously prove that L.BFGS is a stable and efficient algorithm that complements the widely-used Sinkhorn's algorithm and its stabilized versions. In particular, L-BFGS typically converges faster for Sinkhorn problems with a weak regularization. It is worth noting that the proposed analytic differentiation method can be combined with different forward algorithms, and reasonable scheme is to use L- BEGS for weakly regularized OT problems, and to choose Sinkhorn's algorithm otherwise. The differentiable Sinkhorn loss has many potential applications in generative modeling and permutation learning (Mena et al., 2018), and we anticipate that the technique developed in this article would boost future research those directions.
9

---

REFERENCES
Jason Altschuler, Jonathan Niles- Weed, and Philippe Rigollet. Near-linear time approximation algo- rithms for optimal transport via Sinkhorn iteration. In Advances neural information processing systems 30, 2017.
Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In International Conference on Machine Learning. 2017.
Richard Byrd, Jorge Nocedal, and Ya-Xiang Yuan. Global convergence of cass of quasi-Nous Newas- methods on convex problems. SIAM Journal on Numerical Analysis, 24(5):1111-119,, 1987.
Dylan Campbell, Liu Liu, and Stephen Gould. Solving the blind perspective-n-pinn problem end- to-end with robust differentiable geometric optimization. In European Conference on Computer Vision, pp. 244-261. Springer, 2020.
Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in neural information processing systems 26, 2013.
Marco Cuturi and Gabriel Peyré. Semi-dual regularized optimal transport. arXiv preprint
181 2018.
Marco Cuturi, Olivier Teboul, and Jean-Philippe Vert. Differentiable ranking and sorting using optimal transport. In Advances in Neural Information Processing Systems 32, 2019.
Marco Cuturi, Olivier Teboul, Jonathan Niles-Wer Weed, and Jean-Philippe Vert. Supervised quantile normalization for low rank matrix factorization. In International Conference on Machine Learn- ing. 2020.
Marco Cuturi, Lactitia Teboul. Optimal transport tools (OTT): JAX toolbox for all things Wasserstein. arXiv preprint arXiv:2201.12324, 2022.
Yingtao Tian, Charlotte Bunne, Geoff Davis, and Olivier
Pavel Dvurechensky, Alexander Gasnikov, and Alexey Kroshnin. Computational optimal transport: Complexity by accelerated gradient descent is better than by Sinkhorn's 's algorithm. Interna- tional conference machine learning, pp. 1367-1376, 2018.
Marvin Eisenberger, Aysim Toker, Laura Leal-Taixé, Florian Bernard, and Daniel Cremers. Scalable Sinkhorn backpropagation. OpenReview preprint id:uR7707LL55h, 2022.
Yizheng Fan. Schur complements and its applications to symmetric nonnegative and Z-matrices.
Linear algebra and its applications, 353(1-3):289-33) 0009993307200000 2002.
Jean Feydy, Thibault Séjourné, François-Xane Xavier Vialard, Shun-ichi Amari, Alain Trouvé, and
Gabriel Peyré. Interpolating between optimal transport and MMD using Sinkhorn divergences. In The 22nd International Conference on Artificial Intelligence and Statistics, pp. 2681-2690, 2019.
Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z Alaya, Aurélie Boisbunon, Stanis- las Chambon, Lactitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, et al. POT: Python optimal transport. Journal of Machine Learning Research, 22(78): 88, 2021.
Aude Genevay, Gabriel Peyré, and Marco Cuturi. Learning generative models with Sinkhorn diver- gences. In International Conference on Artificial Intelligence and Statistics, 2018.
Aude Genevay, Lénaic Chizat, Francis Bach, Marco Cuturi, and Gabriel Peyré. Sample complexity of Sinkhorn divergences. In International Conference on Artificial Intelligence and Statistics, 2019.
Nicholas J Higham. Accuracy and stability of numerical algorithms. SIAM, 2002.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE,
1998.
Tianyi Lin, Nhat Ho, and Michael Jordan. On the efficiency of entropic regularized algorithms for optimal transport, Journal of Machine Learning Research, 23(137): 2022.
10

---

Dong C Liu and Jorge Nocedal. On the limited memory BFGS method for large scale optimization. Mathematical programming, 45(1):503-555 1989.
Ziwei Liu, Ping Luo, Xiaogang Wang, and Xiaoou Tang. Deep learning face attributes in the wild.
In Proceedings of the IEEE international conference on computer vision, pp. 3730-3738, 2015.
Giulia Luise, Alessandro Rudi, Massimiliano Pontil, and Carlo Ciliberto. Differential properties of Sinkhorn approximation for learning with Wasserstein distance. In Advances in Neural Informa- tion Processing Systems 31, 2018.
Albert W Marshall, Ingram Olkin, and Barry C Arnold. Inequalities: Theory of Majorization and Its Applications, Second Edition. Springer Science & Business Media, 2011.
Gonzalo Mena and Jonathan Niles-Weed, Statistical bounds for entropic optimal transport: sample complexity and the central limit theorem. In Advances in Neural Information Processing Systems 32, 2019.
Gonzalo Mena, David Belanger, Scott Linderman, and Jasper Snoek. Learning latent permutations with Gumbel-Sinkhorn networks. arXiv preprint 18X2v:880..00660 2018.
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. In 2009 IEEE 12th inter- national conference on computer vision, pp. 460-467. IEEE, 2009.
Gabriel Peyré, Marco Cuturi, et Computational optimal transport: With applications to data science. Foundations and Trends@ in Machine Learning, 11(5-6)355-607, 2019.
Tyrrell Rockafellar. Convex Analysis. Princeton University Press, 1970.
Bernhard Schmitzer. Stabilized sparse scaling algorithms for entropy regularized transport problems. SIAM Journal on Scientific Computing, 410):A1443-A11481,2000 2019.
Richard Sinkhorn. A relationship between arbitrary positive matrices and doubly stochastic matri- ces. The annals of mathematical statistics, 35(2):876-87, 1964.
Richard Sinkhorn and Paul Knopp. Concerning nonnegative matrices and doubly stochastic matri- ces. Pacific Journal of Mathematics, 21(2):343-345 1967.
Gui-Xian Tian and Ting-Zhu Huang. Inequalities for the minimum eigenvalue of M-matrices. The Electronic Journal of Linear Algebra, 20:291-302, 2010.
I Toolstikhin, O Bousquet, S Gelly, and B Schölkopf. Wasserstein auto-encoders. In 6th International Conference on Learning Representations (ICLR 2018), 2018
Cédric Villani. Optimal transport: old and new, volume 338. Springer, 2009.
Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-MNIST: a novel image dataset for bench- marking machine learning algorithms. arXiv preprint arXiv:
2017.
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.
11

---

A EXPLICIT EXPRESSIONS FOR CONSTANTS
We first define a few user constants for the L-BFGS algorithm. Let my be the maximum number of correction vectors used to construct the BFGS matrix Bg, and с (0, 10,1/2), 2 E (c1, 1) two constants related to the Wolfe condition: we assume the L-BFGS algorithm uses some line search algorithm to select the step sizes og that satisfy:
[( + ards) S (1xs) + crang de,
g( + =x 2 c29 di.
where S()) and g(-) stand for the objective function and gradient function, respectively, is the k-th iterate, 9s g(xx), dN = = Ih is the search direction, and B, is the BFGS matrix that approximates the Hessian matrix. mig C1, and 02 are selected by the user.
let (300) be the initial value, and μ Ma and u, max, Mej. 1,...." n. Then we define the following constants:
For Theorem 3
Uc =b
[
(1575m-11 11))
max
Σ
Oe log -1 1 1(310)
m-1
A, = 1) log a; - Ue - 100g c (-M(Mit
+
Σ
e
-AMes
i i= 1......
2=1
M1 = X m7 + 2
min 1sisn
((m4-Mim)
M2 =
1
Σ
n
e A(A, Mim)
2(n m + 1)
i=1
M3 M₂ log M1 1,
M4 = m - mo - mg log M1 log(1 + mg.)2)
C = 2/M1,
C2= = 2 M2
r= 1 e1 - M2 /2e (M)++MM)
For Theorem 4
Cs 4) - [r + C.)]
-
C. 20 (1luell /I + 04)1 + +3|p|||]]
C, |//,||| |l + 2C,
||a の (T*s
B ADDITIONAL EXPERIMENT DETAILS
B.1 SETTINGS OF THE MOTIVATING EXAMPLE
Consider a small problem of n 90 and m 60. Let x; 11,
be m. The
equally-spaced points on j0, 5). and similarly define y, 5(j - 10)11mmmm
cost matrix is set to M,j = (x; y²)²² and the weights a and are specified as follows. Let fi be the density function of an exponential distribution with mean and /2 be the density function of mixture of two normal distributions, 0 (..8))008--80 N...-N((,, 0.25). And then we set ă. = fi (xx), b, = f2 a; 0 >h- dg., and bj b, Σ m2 Ā- be-
We fix the regularization parameter 7 Sinkhorn plan 7 is visually close to the OT plan p*. In Figure 5 we show the Wasserstein and Sinkhorn transport plans under different values of n. can be seen that when 7 S 0.001, TX visually indistinguishable from P
0.001. This value is selected such that the resulting
We compute the true T using the e-scaling algorithm (Algorithm 3 Schhiizerr 2019). This algorithm is typically accurate, but it requires solving sequence of Sinkhorn problems with in- creasing A'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.
12

---

Wasserstein
Sinkhorn 0.1
Sinkhorn 0.01
Sinkhorn 0.001
Sinkhorn 0.0001
Figure 5: Visualization of Sinkhorn plans with different v) values.
B.2 FORWARD PASS 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 X1..... X, fi (x) and Y1...., Y. ~ f2(3) from some specific distributions fi (I) and fz(3), z, y E R², and construct cost matrix as = (Mij), where M,j = || Yşll². The weights are fixed as a = n -1,, and b m 1m- For each of the five methods compared in the motivating example Section 3.1 let T' = (Ty) 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 on the transport plan,
Errptan (T)
(Ty TT2)2
and the error on the Sinkhorn loss value,
ErTloss (T)) ((T, M) (T* MY = I
For each configuration the experiment, we simulate the data 100 times, and visualize the distribu- tion 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 y 0.01, 0.001, and for each set a specific maximum number of iterations for all algorithms. Two data generation models are considered:
(a) Both (x) and f2(y) are multivariate normal distributions N(0, Tp);
(b) Both Ji (xx and f2(v) 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 N (1, 0.04) + 0.8 (3, 0.25).
The final results are demonstrated in Figure 6 where all the errors are shown in the log-scale.
In Figure 6 many boxplots 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 generates no NaNs values explicitly for v) = 0.001, does not give any meaningful results either. This implies that numerical stability is critical issue in computing the Sinkhorn loss.
For Sinkhorn-logg gives reasonably small errors when the regularization is large (n 0.01) with sufficient number of iterations, but its accuracy quickly deteriorates when ń) decreases to 0.001. Moreover, we can find that Sinkhorn-log sensitive to the limits on number of iterations. For example, when maximum number iterations, but we restrict the limit to 200, the error can be as large as 10-2 even 10". depending the data distribution.
1 and
0.01, the loss value error can be small as 10 given 1000
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 not needed. These find- ings demonstrate the advantage of the advocated L-BFGS method in both numerical stability and accuracy.
13

---

Errors on the Sinkhorn trans plan
Dim
Dim:
m: 50
Disailaddo donnn nemaal
Log Scale Emorr
nyMaxiter
-6-
0.01/800
0.01/1000
Disalladion oig mxx
0.001/600
-2-
0.001/2000
Swithorn
Sinktom- Ang
Stabilind
Greenthem
L-BFGS
Sininorn
Siatham--oe
Statalied
Greenshom
L-BPGs
Smihhem
Sindtm- Rog
Sabatod
Greenkham
L-BFGS
Algorithm
(a) Sinkhorn transport plan
Errors the Sinkhorn loss value
Dimep-1 m1
Dim 10
Dim: p 50
Dicoibution namal
Log -cw---aaaeeene Emor
-6-
ryMaxiter
-8
0.01/300
0.01/1000
Distribttons mk
0.0016600
0.001/2000
-2
Sinieom
Sinihamt-og
Stabilired
Griendeen
L-BFGS
Sinkion
Sintharn-ig
statiliied
Greektem
L-OFGs
Sistten
Stowham-log
Statatired
Greenbbam
L-BFGS
Algorithm
(b) Sinkhorn loss value
Figure 6: Comparing different algorithms on the errors of Sinkhorn transport plan and loss value.
The missing boxplots indicate that the corresponding results are NaNs.
B.3 RUNNING TIME OF FORWARD AND BACKWARD PASSES
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 forward pass, and analytic differentiation in the backward pass.
(b) "Implicit" implemented in the OTT-AXX library (Cuturi et al. 2022): Sinkhomn-og in the for- ward pass, and implicit differentiation in the backward pass.
(c) "Unroll" implemented in the OTT-JAX library: Sinkhorn-log in the forward pass, and unrolled automatic differentiation in the backward pass.
We use the second data generation model in Section to simulate data points, and use the three algorithms above to compute the Sinkhorn loss and its derivative with respect the cost matrix. For each configuration, we randomly generate the data 100 times, and compute their mean forward time and mean total time. The stopping rule implemented in the OTT-JAX library is [|T,,, al +
14

---

|T",, bll Eott. and one of the terms would be exactly zero in the last iteration. The stopping rule for L-BFGS 1TT boloc < Etbsgs- To account for such a difference, we set Etb =
and let Coet / max n. flbggs. In fact, this setting favors the competing method, as its stopping criterion strictly weaker than the proposed one. 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 in Table 2 and Table 3 where the former uses 1000 maximum number of iterations, and the latter uses 10000.
Table 2: Running time of three algorithms for differentiating the Sinkhorn loss, with maximum 1000 iterations.
Mean forward time (ms)
Std. of forward time
Mean total time (ms)
Std. of total time
Converged
Analytic
19.07
2.11
19.33
2.13
100
m = n 64
p 8
Implicit
27.77
0.69
39.66
3,41
0
7) = 0.1
Unroll
27.53
0.13
101.30
0.36
0
Analytic
23.21
3.45
23.48
3,46
100
64
p = 8
Implicit
27.33
0.05
50.52
5.19
0
7 0.01
Unroll
27.39
0.12
100.66
0.39
0
Analytic
23.79
2.96
24.19
2.96
100
128
p = 16
Implicit
29.81
0.61
58.21
4.91
0
? = 0.1
Unroll
29.75
0.06
107.18
0.19
0
Analytic
32.24
4.99
32.68
5.00
100
m = n 128
p 16
Implicit
29.73
0.05
85.56
7.48
0
7 0.01
Unroll
29.65
0.08
106.66
0.46
0
Analytic
28.60
3.82
29.40
3.85
100
256
p = 32
Implicit
33.69
0.32
82.53
4.21
0
1) 0.1
Unroll
33.60
0.23
123.25
0.07
0
Analytic
43.94
6.08
44.86
6.12
100
256
p = 32
Implicit
33.75
0.20
137.46
11.66
0
?) = 0.01
Unroll
33.78
0.10
122.99
0.27
0
Analytic
41.35
5.64
43.05
5.53
100
m = n = 512
p 64
Implicit
46.06
0.05
130.96
11.97
0
7 = 0.1
Unroll
46.02
0.04
167.64
0.22
0
Analytic
67.76
9.65
69.42
9.65
100
512
p = 64
Implicit
46.08
0.05
230.34
11.72
0
7 0.01
Unroll
46.09
0.05
167.68
0.20
0
B.4 MORE EXPERIMENTS ON STMULATED DATA
We experiment on more simulated data to compare the computational efficiency of Analytic and Implicit. We do not include Unroll since it is too time-consuming. The results given in Figure as analogues of Figure 3 Three more simulated data sets are studied, and we also include the Sinkhorn loss metric evaluate the model performance over time.
15

---

Table 3: Running time of three algorithms for differentiating the Sinkhorn loss, with maximum 10000 iterations.
Mean forward time (ms)
Std. of forward time
Mcan total time (ms)
Std. of total time
Converged
Analytic
19.01
2.10
19.27
2.10
100
64
p = 8
Implicit
267.95
24.57
277.86
24.58
12
? = 0.1
Unroll
269.49
24.76
986.46
90.67
12
Analytic
23.28
3,46
23.54
3,48
100
64
p = 8
Implicit
275.42
0.34
284.84
4.74
0
7) 0.01
Unroll
275.54
0.49
1011.70
0.75
0
Analytic
23.76
2.92
24.16
2.90
100
m n1
128
p 16
Implicit
297.36
0.34
320.29
4.32
0
7) = 0.1
Unroll
297.93
0.63
1074.47
1.33
0
Analytic
32.21
5.00
32.68
5.01
100
128
p 16
Implicit
297.79
0.30
324.95
9.99
7 0.01
Unroll
297.90
0.32
1073.24
0.84
0
Analytic
28.72
3.84
29.57
3.83
100
256
p = 32
Implicit
342.23
1.53
386.80
9.22
0
1) 0.1
Unroll
343.18
0.70
1256.10
2.94
0
Analytic
44.45
6.14
45.32
6.18
100
m n1
256
p 32
Implicit
343.57
2.01
411.12
18.92
0
7 = 0.01
Unroll
341.00
2.02
1253.67
3.33
0
Analytic
41.50
5.60
43.27
5.53
100
512
p = 64
Implicit
469.04
0.77
552.64
15.31
0
7 = 0.1
Unroll
465.03
1.82
1705.51
6.36
Analytic
67.97
9.69
69.68
9.67
100
512
p = 64
Implicit
464.59
1.04
624.28
26.48
0
7 0.01
Unroll
464.16
0.55
1707.04
3.31
0
B.5 ARCHITECTURES FOR DEEP GENERATIVE MODELS
For the deep generative models on MNIST and Fashion- data in Section 7.4, the architectures of the generators are given in Table where FC2 stands for fully-connected layer with output dimensions, and Convo Convc.k...piin.pov kernel size stride 8, input padding pin, and output padding pout.
the transposed convolutional layer with output channels,
Both models are trained using the Adam optimizer with learning rate 0.0001, on mini-batches of size 600. The Sinkhorn regularization parameter is set to = 0.1, and the training process consists the squared L distance to construct the cost matrix, and
of two stages. In the first stage, in the second stage we switch to the L' distance. The intuition that the squared L distance has smoother derivatives, thereby making the training more stable in early steps; the other hand, the L' distance makes the generated images sharper. The two stages are run for 20 and 30 epochs, respectively. The Wasserstein distances values in Figure are computed in the first stage.
16

---

31
Hond
Redd7
Matred
Motrod
Mithod
Mochoo
Diowe This
Diovel becones)
Elussed Time becomes"
2, ne
Hone Z. p-0.
Plos ge
RextZ gutt
Netted
Veenod
Wetted mnalcc
Methos
Disuse Thme taecanda
Dipoed Tre baconds)
Sesird Time denardul
Elossed Tixe becones
Andd 1-6.
Fink ps
ReddZ
Time Rexiz guddt
Bland co
Hund g-1.)
Matrodd
Midival
Mittred
Methoo And
Daziee Time becandi
Deseed Tim becones
Depied Time Oetands
Elicaed Trxe Seccomes
iolldt
thoiic
Time
Sampked mo-lll
Sampleelz que
Samuled meila
Sampled
Metrod
Metrod
Methos'
Dasond Timn (sscandal
Ulapeed Tiee Secoeds)
B.0zoo 157a Gracaniil
Elodaed Yeat (0conna
Figure 7: More experiments on simulated data. The plots are analogues of Figure 3 with three additional data sets and the extra Sinkhorn divergence metric to evaluate model performance.
B.6 APPLYING THE SINKHORN LOSS TO WASSERSTEIN AUTO-ENCODER
Finally, we illustrate the application of the Sinkhorn loss in an Wasserstein auto-encoder (WAE, Tolstikhin et al.] 2018) generative model, based on the CelebA data set (Liu et al.) 2015). The training data consist of 20000 human face images, each of size 64 x 64. The WAE model learns an encoder Q TRP IR" and decoder G R" RR solution to the following optimization
17

---

Table 4: Architectures of the neural networks for MNIST and Fashion-MNIST data.
MNIST
Fashion-MNIST
Layer
Output shape
Layer
Output shape
Input
16
Input 2
32 6272
FC3136 + ReLU
3136
FC6272 + ReL.U
Reshape
64 X 7 X
Reshape
128 x X 7
C,, ReLU
32 14 14
Conv64,4,2.100 + ReLU
64 x 14 14
Conv + ReLU
32 x 29 x 29
ConVe4,422.1,1 + ReLU
64 x 29 x 29
1,,21,4
x 28 28
ConV1,4,1.2.0
1 x 28 28
problem:
min E Exx |X- G(Q(()))++ & D(pz, Q.G
where Px. ()()), and pz are the distribution of data points X, the distribution of Q(X), and pre-specified latent distribution, respectively, and is a regularization parameter to balance the two terms. The first term in (7) is the reconstruction error, and the second term quantifies the divergence of the distribution of Q(X) to the latent distribution pz- We simply take pz to be multivariate standard normal distribution, and use the Sinkhorn divergence (Genevay et al. 2018) to define D(-):
D(u,v) = Sx(A,)) / Sx(,
1 Sx(,
where μ and are two discrete distributions, and with slight abuse of notation, Sx(u, v) is the Sinkhorn loss studied in this article. In actual implementation, and v are Diracs of data points from pz and P)(X), respectively. To generate new images, sample latent data points Zz,.... Zn from the latent distribution pz. and then pass them to the generator obtain images Y; G(Zz), i 1.,....,..
Since the focus of this article is the computation and differentiation of the Sinkhorn loss, do not attempt to build and train the model with full complexity. Instead, to compare computing time, only 10 epochs for illustration purpose. The architectures of the encoder and the decoder given in Table 5 and we set the latent dimension to = 64, the WAE regularization parameter to E and the Sinkhorn loss parameter to 1) 0.1. We use the Adam optimizer with a learning rate of 0.001 and mini-batch size of 500. We have found that the Implicit method generates NaNs after 72 mini-batch iterations, so we only show the results for the proposed Analytic algorithm and the existing Unroll method. Figure X () shows the training process of the WAE model based on the Sinkhorn divergence for 10 epochs, and Figure(8 shows the randomly generated images using the Analytic method after 50 epochs. It is clear that the proposed Analytic method is more efficient than the Unroll method in training.
Table 5: Architectures of the neural networks for CelebA data.
Encoder
Decoder
Layer
Output shape
Layer
Output shape
Input
3 x 64 x 64
Input
64
ConV32,422100 + ReLU
32 32 32
Convese.4.Luu
ReLU 256 x
ConVe4.22110 + ReLU
64 x 16 x 16
Convi28,42210 + ReLU 128 x 8 8
Conv 120,2,2.22 ReLU
128 x 8
Conve4.42,1.0 + ReLU 64
16 x 16
20. ReLU
256 x 4 4
Conv32.421100 + ReLU
32 x 32 x 32
ConV61,4..000
64 X X
Conv3,1.2.10
3 x 64 x 64
18

---

Training using Divergence
Analytic
Unro
10
Loss Fundtion Value
(secords)
(a)
(b)
Figure 8: (a) Training process of WAE for the CelebA data for 10 epochs, based the Sinkhorn
divergence. (b) Randomly generated images using the Analytic method after 50 epochs of training.
B.7 COMPUTING ENVIRONMENT
All of the experiments in this article were conducted on a personal computer with Intel 19- 1190K CPU and an Nvidia RTX 3090 GPU.
For the motivating example Section 3.1 and the experiments in Section B.2 results for Sinkhorn's 's
algorithm, Sinkhorn-log, Stabilized, and Greenkhorn are computed using the Python Optimal Trans-
port (POT) library (Flamary et al., (202T). For the experiments in Sections 7.3 B.3 B.4 and the Implicit and Unroll algorithms are implemented in the OTT-JAX library (Cuturi et al. 2022).
C PROOFS OF THEOREMS
C.I TECHNICAL LEMMAS
In this section we state few technical lemmas that are used to prove our main theorems. Lemma to Lemma below are standard conclusions in vector calculus, Lemma and Lemma (5) are derived from the eigenvalue theory, and Lemma 7 and Lemma (8) are related to the Sinkhorn problem.
We first introduce the following notations. For any æ E IR, let [x]+ = max(#, 0}. For a matrix A (Aig) R2 the vectorization operator, vec(A), creates vector by stacking the column vectors of A together, i.e.,
vec(A) = (a11,.... 4.1, 112 an2, 01m. anm )
For two matrices A = (aiz) Roxm and B E Rexq, the Kronecker product of A and B is defined as
a11 B
01m
B
A 0 B
Qx B
a B
For differentiable vector-valued function f IR" + is defined
Th the partial derivative of f with respect to
(fi(w)
0f((z)
0f(x) T
al (z)
2 fm (x)
We use L, to denote the n x n identity matrix, and max (.) and minn..) stand for the largest and smallest eigenvalues of some symmetric matrix, respectively.
Lemma 1. Given matrices A E mxx and B B
vec( (BB (T, 6 A) A)vec( B) (B" 60 ve)eec(A).
19

---

Lemma 2. Let f m
be two vector-valued differentiable functions of R² and g I I
I E R.. Then
dr / (x) g =g(x)² 0f(z)
①
f(x))
T
Or
dr
+ R$ and g Rmm Rp" be two vector-valued differentiable functions of Lemma 3. Let f R
E Rm. Then
0f(z)
+ (1, @ f(x)) 0g(x) dr
vec(f(z)(()))) (g(a) 11) dx
dr
Lemma 4. Let A and B be two n. X positive definite matrices, and let 2 ... 2 On 0 and 31 2 ... B, > be the ordered eigenvalues of A and B, respectively. Then for any x E IR"
T A BA S 0131111112
Proof. Let Uixn xTAA'/ and then u = UUT x a|||=| and UBUT xT 1 ///// By Theorem A.4 in page 788 of Marshall et al. (2011). we immediately get
UBUT trUBBU") S Biu 0//111211
Lemma 5. Let A and B two symmetric matrices of the same size. Then
mmin (A) - m (B) S Omin (Ain(AA + B) S max (A + B) S maax (A) + Faxx (B).
Proof. Using well-known identity Tmax (max(A) = max|| x' Ax, we have
0max (A + B) = max x ( B) max x" T Axr + max ||]|=1
x BB = Tmax + Tmax (B).
11
||||=|
Applying the inequality above to -(A+ + B), we get the result on the opposite direction.
Lemma 6. Let (I1,..., zn) and y (y1...., Vn) T be vectors. Define LSE(z) log Di= e"t, and then for any x, y E R
- | o| < min (x: vi) LSE(x) - LSE(y) & max (x:: - y1) S ||0 ylooo.
Proof. It is easy to find that vgLSE(() softmax(x) = (s ST) the softmax function, where 8; E#: Σn k= c (0, By the mean value theorem, have
LSE(x) LSE(v) softmax(uz + (1 - w)v) y)
for some < Let softmax(wz (wx + (1 w)y). and then
LSE(x) - LSE(y) = Σ 21(1; v1) < max (x: vi) Σ
max (x; vi) s 1| vllocs
iml
and similarly, LSE(x) - LSE(Y) 2 min;(a; - y:) 2 -|, ylloc.
Lemma 7. Let f(3) be defined as in 16). and let μ M a and u. max, Mij. = 1,... n. If f(3) < c, then have max; Bi < U. and min; Bj Le, where
U. = b-d
max 1
+
a; log +
(8)
1=1
Le
min 1575m-1
Σ
a32 + 17
a log a; 77 +
1=1
20

---

Proof. By definition,
12
m
n
f(3) = 11
Σ
0i log
Σ 1=1
c M00 1 Mes)
Σ
a log a4 - BT + n
1=1
If f(3) S c, then
Σ a4 log 04 1 17 Σ aj log
e M(3 MMy)
- -.
= + 1
jml
By definition we have B. 0, and let J = arg maxigjsm mm 11 By. Then
m
60 2 log (M(d)-Mal)) -3 Σ
a: MMs) Σ 3zb3
i=1
m-1
2 BJ M,Mis - 33 2 bj = 0mm33 - PJt 2 bm
max
Mj:
1575
j=l
which verifies ( by noting that max, 3g = [$
Next, let K arg mini --jmmm By.. We can assume that Bx < 0, since otherwise the trivial bound
miny Bj = p. = 0 is already met. Consider the sets S+ { : Bg > 0} and S_ = {j : B3 < 0}.
Then clearly,
b = 3KbK +
Σ
8yb5 +
Σ
Bybs S 3xbx + [Bul+
Σ b, + S BxbK + +B,]+ bm)
j
j jES-
j jES_
JES+
Also note that log(
25= E²5) 2 maxj Tj for any FI,..., Tn E R, so
m
Og log
10a- Mis)
2 1) i=l
aj log
Σ j=1
c
M(3--u)) uz)
17 02
log
Σ
iml
j=l
m
n
n
= 17 log
a 2 max (A3;)
Σ
a,u;
j=1
1=1
max 3g
02 [3,]+
a
1=1
As a result,
co 2
ag log
e MB, Mis)
- -8 [351+ - Σ
a.u; - G,ubKBkbK - 115]+
bm))
1=1
1=1
1=1
n
n
bmm(B11] [8,]+
Σ
0 - Bxbx 2
Σ
0g44 Bxbx,
and then
2 2-
2-
n aru; + Co
Bx 2 -bx
Σ
auu; + Co
min
by
Σ
15/sm-1
1=1
which verifies 99.
Lemma 8. Let T be an n x m matrix with strictly positive entries, and suppose that n 2 m. Define
T1 m. V = T and
H (iag(())
T diag(E)
= diag(z) T diaggu))
1
21

---

Then
maxx
Fma 5 max Vjj 15/5m-1
m +2
min (D) 2 0 min (H) 2
min T Isi5n
2(n
+)
m-1
Mmin (D)
min
Σ
Dig
min
M4
Ti TijTim- im.
15
15i5m-1
==1
Proof. Consider the matrix S = H - sJ, where J is an (n + m - 1) x (n + m - 1) matrix filled
of ones, and is positive scalar. Let
R
Σ
Skyl:
h - m 1.
j
Suppose & silnisisnsssssi and then for h: ,,..,,,h n, is easy to find that
m-1
R₂ = (n- - 1)8+ +
(Tkj 8) (n 1)8 + pk Tkm (m - 1)s = (n m)) + pa T Tamm
j=l
and for k = n1 + 1,..., + m - 1,
n
Rg (m 2)s
(Ti,k-n 5) s)= (m h 2)s v ns (m 2)8 + vk-
Then it is easy to see that
[ik R2 T (n m)8,
k 1, ,...,,
Suk Rs
Vn-n- Rs = (n + 2 m)),
k = n1 + 1,...,77 ,
Let
min Tim (n m)s = (n + 2 - -m)s, Isi5n
and then s minigisn Tim (2n 2m + 2), and
12 1 m 2
Ska R 2 L = 2(n
min Tim >> 0 1si5n
+ 1)
for all k. By the Gershgorin circle theorem, every eigenvalue S must be greater than L. Since
H = S sJ and J is nonnegative definite, we also have o (H) 2 L > 0, implying that H is positive definite.
For the second formula, it is easy to find that the D matrix is the Schur complement of the block diag(u) of the H matrix, So by Theorem 3.1 of Fan (2002). we have amin ()) 2 Gmin (H) and mmax ()) mmax (diag())
maxigjsm-1 Vj--
Finally, let e = max1 sinxissmm1120 V/. and then D can be expressed as clm-1 B, where T diag(μ) (T) diag(e1 (iag(c1m-11 /0) is a matrix that have nonnegative entries. In addition, we have proved that / is positive definite, so D is nonsingular M-matrix by the definition in Tian & Huang 2010). Then Theorem 3.2 Tia & Huang (2010) shows that
Mmin (D) 2
min Σ Dij- 13i5m-1
j=1
Let 5 D1 and then clearly min
ΣØ=1 Dij min; 8, Note that
= 8 = m1 - di)gg(( 111 8 - nT didig(()) - (TI (T1m - -Tm) (u) "T dTiag
TT -iag())) Tm 1 T diag
I diag(a)) -1T
where T,, stands for the m-th column T. Therefore,
min 1515m-1
8, min
Σ Mb Tim. 15/5m-
22

---

C.2 PROOF OF (
Let T = ex 2 0 - MI. and then it is easy to find that vaL(a, 3) - a 71, and V 3 [ B)
3 --T 1,, Since o" (3) = arg maxa, L(a, B), we find that as = a a3();; is the solution to the equation 71m O. By definition, we have
0 Max+8; -M,) e
MBs M M.,),
i= i=
so the solution is a; 1 log 0; 100g - A(0g-MMy)) Since 71mm a, we immediately get 1
so C(a = 8 b T1, = a² (8) 3T 7
and we get the expression for f(3) (3), 3). Finally,
VzL(a"(3), = [8a* [3)
Val(a, B) ())a= + ViL(a, 33 (a)a(())
0,T
Since a²(3) arg m L(a, (3) implies that Val(a, (a
0, we have
(3), = VzL(a, 3) (s) and hence
v Vaf(3) - ((aaa"()) -(())) -b. b.
C.3 PROOF OF THEOREM
By definition we have
Sx (TT,M) = vec(T*) vec(M).. ve)(AM.
so Lemma gives
a b) = V )vec(T")
@vec(M)T + vec(?")) TT Övec(M) dvec(M)T*
(10)
Ovec(M)²
Obviously, (M)
is the (nm) (nm) identity matrix ((m)). so the second term of 10 essentially vec(T*) and the remaining task is to derive /
Let R = a" @ 8 M a IT m + 1n - M, and then T = ex [R]. Using the chain rule of derivatives, we have
@vec(T") Tv)) Övec(
(11)
@vec(M)T @vec(RyT @vec(M)T
is easy find that ( is an (nm) (nm) diagonal matrix with diagonal elements vec(XI")),
vec(M) Övec(T*) Övec(R)T
(12)
Furthermore,
Övec(R) R)
Ovec(a I + 1T M)
Ōvec (MyT
Öven(M)T
da
05*
[l % In)) Dvcc(M)T + (Jm ? 1,)
Övec [(mm))
(13)
where the second identity is an application of Lemma Combine get
12 and 13 and then we
vec(M)T )vec(T")
dvec(M)T Avec(M o ovec(M@TT) T Övec(R) R)
Ovec M)T
Oa*
Xvec(M ( (T)) In) Ovec MYT
03*
+Xvec(M @ * (T @ 1n) Övec(M)T
-Avee(M C -Nvce(MoTT))
23

---

which is the first term of (10).
Using the identities Lemma] we have
vec(M (M o (1m
In)
[11m @ (n) vec(M o Vec oT))
vec(
[1 @ (
= [vec((M o T [mm! 11 = =,,
v o TIm & 1n) = (Im & 1n) T vec [ ④ TT) = [L @ VTM
|voc | (M o )]] = Pe-
Since we have set 9m O,
10
simplifies to
OS^ a, b) @vec(M)T
03 @vec(M)) + Me @vec(M)T vec(M @TT) + vec(TT)².
Oa
Hr
(14)
Let w² = ( and then the main challenge is to calculate Ow /0
First, note that the optimality condition for (a", , arg m [(a, B) is
VaL(o, (() = 0, VaL(a, ((a O.
(15)
Section C.2 has shown that Val(a, 8) = a - T1m and V VzL(a, 3) =b - b-TT1, TT,,- Moreover,
Vi, B) - Adiag(71m)
(z Adiagg(T 1..)) 1n)
V B (Va((a -XX-.
Define the function
F(w, M) = (vo ff
- T1m TT )
(16)
where (a= BT and then wi satisfies the equation F(w", M) 0, indicating that w" is implicitly function of M, written as w" = w(M) By the implicit function theorem,
Ow(M) Ōvec(M)T
@F M) Ow
0F M)
F F
dvec (M)T
W==""
Note that
Va v 8 (VaL(a, 8)]
V (Va((,,))))
)
F FT.
=-
A B
D ))
BT
Then by the inversion formula for block matrices, we have
= (
B
=->-1 ( A A BÂ BT A
-A -A-IB -1
)
BT
D
-A A-1
A-1
where Á D - ATAA B.
For g = (1T (TT)T, the vector 5 = (s hT = -XF ' has the
following expression:
8, A-1 A M4. Ι 1 Nr.
8, = A + A BÂ BT 1 pr A-1 Pe,
A U4 A B3,.
After some simplification, we obtain
Â diag( (T 1,,) T didggg( y7
8, = 1 i Á 1 3 Hr)
8a (Tlm)= - ① Hr (T1m)) c (Tsu).
24

---

G HM
Next, partition FMr as FM = definition,
where GM [ax(nm) and HM ( [((=-1x(am), By
aG, (veT(M)T
GM
G.
aG @vec(M)T
so
aG:
S0,
i k: i= k:
0M
(TRRI
This indicates that GM X (diag( (T1). T' Similarly, for Hij = Sin1 Tij.
diag(T (Tm)), where Ti...., T are the column vectors
aH, dvec(M))
0H;
[0,
i j=l'
HM
0MR
ATT::
aHm-1 (vec(
implying that
7f
HM =
TT
0
m-1
As result,
da- 1r (vec(M)T
03 'ee Övec(M)T
T
Ow"
- T. @vec(Mit
- T o. FF F
T
x G-
x HM ) ) (5T, 2T)
A G- G
XF
Ae
HM
((sg O T2) ".................. (su O Tm + (T
80 T 1 0
[vec (diag(su)T + (ddaag(s))))
(17)
Finally, substitute 11 into (4 and we have
2M,,Maan
(vec (diag(s.)T* + +T" vec(M T)) vec(T)
@vec(M)T
Transforming back to the matrix form, and we obtain
SM, a,
A(su + - -))
9M
Replacing T with T and noting that T 1m T get the stated result. The positive definiteness of the Á matrix is a direct consequence of Lemma
C.4 PROOF OF THEOREM 2
In the proof of Theorem
we have already shown that
2 0. f(a IVV -- (diag (1m) TT
T
diag( (
where T = exfa ① 3 Mj. Plugging a" (13) into L(a, 8). and then V V2 (9), (3) is the Schur
complement of the top left block of 2 a.d v2,,a(( C(a, 8), given by
Vacla = - diag(? (T 1,)) - 7th diag -17)
25

---

Since f(3) -L(a" B), by Lemma find that ()f( is positive definite, so f(3 is strictly convex on B, and hence 3 is unique.
The optimality conditions for "(") are 7 = a and 7 = b, where 7 == =ex(a" 4 8
M). Since Tij exp{\(a; +37 My)] } 2 0 and a, شØ= Tij by -i= Tij. have
Ty S minfa,,by b for all and j. implying that
af + 53 S Us = My + + min logg)
Since 3m by design, have a I Ua. too, i 1...., where Ua: Uim Min + X min (log(a,), lm This indicates that o is upper bounded.
Next, let arg max; Tim Since 1 b implies that bm. -i= Tim S n1Tmm we have
Tim exp(M(a Mmm)} 2 bm/nn,
and hence af 2 Mm + x log(Umm/)) n)). Again, since of + 83 U2 for all i and j. it holds that
83 $ U1 0 S Uı Mm pm X log(bm//))
M1 + x mui(log(ar)) log(()))))))))) (og())))))) Mm - x log(bm/nn)
M1j Mi + - min (loggooonnaana /bm)). /bm))
= Ug, < +00, j = I.....m. m.
(18)
On the other hand, 7 1, = implies that bj = - Túù = eBB
Σ e Ma;; - Miy) for any j. so
log by = A3 + log
Ma; M)))
5 A88 + log
Σ
((Un MMg)
1=1
It is well-known that
2n < log
Σ E
< max ,,, + logn
for any , En E R. so
Bj x log by - -x- log
A(Ua, - Mar)
x 1 log bj- max (Un; Mis) x-1 1logg
2 log(b;/n) max (Un, - Mis) = Lg, > -0, j= m.
(19)
Then 18 and
together show that 1851 < oc.
Similarly, T 1m,
log a; I Aa; + log Σ
MlUs Mar
a > - (a/gga//m)mmmm max (U2, Maj) = La, > -cc, i= 1,...,n
The trivial bounds L.. and U3g are obtained by removing the unknown index
The results above verify that |a Finally, plugging 33 to the objective function, and we immediately get f -
and |87| < and hence ||<x"| < xo and || < co.
C.5 PROOF OF THEOREM 3
Claims (a) and (b) direct consequences of the convergence property of the L-BFGS algorithm (Theorem 7.1, Liu & Nocedal 1989), and we only need to verify its three assumptions. The new results here are explicit expressions for the constants C1, C2, and r.
26

---

First, is twice continuously differentiable, Assumption 7.1(1) of Liu & Nocedal (9989) is verified. Second, f is closed convex function, and we define the level set of as Le {
Rm- ((3 I c). Theorem 2 has shown that f > -00, and when
f obviously Le (3)
is non-empty and bounded. Then Corollary 8.7.1 of |Rockaelllar (9700) shows that Le is bounded for every c. In particular, for a fixed initial value ((0) define L = {B f(3) S )(3(0)), and then
L a bounded, closed, and convex set, which verifies Assumption 7.1 of & Nocedal (1989). Third, let (3) == V f(3), and then in the proof of Theorem we have already shown that
H(3) = > [daag ( 1,,) TT diag (T '7)
where T =axla(B3) MM. Lemma 8 verifies that 2 mmin (H()) 2 X 2(n m 1)
min Tim: 1515n
m (Hnx) S X
max V/z
15/5m-1
where = T On the L set, Lemma 7 shows that max, Br S Ug and miny B3 2 Le with f(((1)). Therefore,
m-
a₁ == 17 logaa a, - 10
Σ 1=1
e (
2 nloga a, - 1log
-MMM +
e (Use-My)
j=1
m-1
= log a: - Uo 19 1l00
e MMMuntt
2Ms
Ai
j=1
Ton = Moi-Min) 2 e M(d--Mam),
On the other hand, Túj must satisfy Ti > and E5m Tij = a, for any and j. so for j
I,..., m 1, we have Tig S as Tim- Therefore,
1/ = WΣTH T S
Σ i=1
(a.
Tim) 1 Σ Tim S 1 Σ (MMg-Mam) j j= i......, 1=1
1=1
This implies that there exist constants M1, M2 such that
Mi ||| 2 < x (3)x S 22
for all x E 1m and B E L, with
Mi = X
n1 - 171
min 1895
Mom)
2(n m + 1)
M2 = x 1 Σ M(i--Min
This verifies Assumption 7.1(3) Liu & Nocedal 1989.
Next, we derive the explicit constants the theorem. Following the notations in equation (7.3) of Liu & Nocedal (1989), the BFGS matrix B₁ for the L-BFGS algorithm has the following expression
B = B(a), B(++)) = B sis B())
319
B
yT
where m min(k + 1, m.). my is a user-defined constant explained Section A and {m} and {si} are two sequences of vectors. We also choose ((0) = to be the identity matrix.
st ||: ? Th = 11 ///l////' sts and Fix I m. and let cos 0₂ = s7 ///ss/(()ss)) [[B((s I]) Pk /T
ge = .r (() 1/. Then it can be verified that
|IB"ssl
1| tr(
= tr(B")) tr()(++)) of B()sss
+
+ Tk.
(20)
of s1
cos2 0x
dct(B((+1)) = det(B(()) ppxp1)))//
(21)
27

---

Define (B) tr(B) - log det (B), and it is known that (B( > 0 for any positive definite matrix B. Equation (6.50) of Nocedal & Wright (2006) shows that
0 < (<(+(t)))= =r((B)) cos² 0k
9k
log del (Bl) log Pk log Qk
Th
1) + (1 (1- cos² 0g
=)) + (Th- log ph -1)
+ log cos? 0g
+ log cos 0s-
Under the assumptions verified above, equations (7.8) and (7.9) of Liu & Nocedal (989 show that
M1 y/" ///lss/////// || M2 and |/91 /// si M2 for every I, so M1 S Ph: M2 and To M2-
Since h(x)= = x + log(z) I 0 for all x > O, we have
(22) Now we show that (BB) can be upper bounded. First, Lemma iimlles that for I O, ..... m 1,
0 < ( + (M2 - log M₁ - 1) + log cos 0g = (Ba)) + M3 + log cos 0g.
O max (B(+1)) I Tmax + +0+ + 11aal1//1/ 5t S mmax + M2;
so maax (Bx) $ + M2 This also implies that q S 1 + ma Next, 20 shows that tr( B&) S m 1 + mg 22, and 21 gives
log det(B(++)) = log det (()) + log log 9b 2 log dct(( + log M1 log( 1 + mo)2),
implying that log det [x) 2 mo [log M1 log(1 1 + moM2)). As a result, we get
(B&) M4
Mo mo [log M1 - log(1 + mg
(23)
Combining
(2 and 23 have cos? 0g (M
Finally, using the argument in Byrd et al. 19887 we have ((++1) - S (r(p(h- f*), where
r= -1- M cos? 02 1 e1(1 -1-1-p0) M2) M4))
and C1. 02 are two constants for the Wolfe condition as explained in Section simply 2/M
A The constant C₁ is
For (c), we follow the analysis in Nocedal et al. 2002 Let g(3) = V VgJ(3), and then g(3²) 0.
By Taylor's theorem we have
((3) (6 (
(24)
where H1 (E) for some & in the line segment connecting 33 and Br. Also,
g( - g(8*) r(3) H2(3 - -3*), H2 [ H( + t130 3)
(25)
Combining (44 and (55) we get
(|o())112 = (3 10- H2 3T)) 2(3 1 304) 3-) H1 - 3r))
(f(B) - fr)
for any E L. Let 12 - B*). and then
[ 2n. T H 1/2 H| 1/2
(f(3) f²).
||0||2
is easy to find that
mmax H 10 min (H1) < MT!, Fmax Hz)) max H2))² < M2.
By Lemma we have x² 11, 1/2 11, 11, 1/2 x < M, M2 |0|| for any x E 1m- and hence 1|ا(3))12 < (2()3) - - for all 99 E L, where C2 = 2 M2. Since (8) & f(0) due to claim
(a), find that ((h)) L for all k 0. Therefore,
(x) (2( 5) s Car f*),
28

---

and claim (c) is proved.
Claims (d) and (e) can be verified as follows. For any & L, recall that T = @ M).
and then 71m a and
1n - -llo0 S 1T - bl = V (f((l)
indicating that
Σ Ti bj s < 1= ,..,..,, m -1
(26)
Then 0 < T5 S
Tú and 0 < Ti S E Ti I b + VV The gradient
V (B)] can be bounded using claim (c), (d) is also proved. On the other hand, 26 shows that
m
a.
Ti max Tij:
==1
and similarly we have
by - [V)((0) & Σ Ti S n max Tij-
Replacing |B)l by its upper bound, and claim (e) is verified.
C.6 PROOF OF THEOREM 4
For matrix Anxm = (asy), define || = maxi, la:yl. and the notation A 2 0 means azj 2 0 for all and j. First, is easy to show that ||A @ BlIF S || B since
| BlI
Σ (an)bn)) <
| = ||
Next, we show that if Baxm 2 O, Cpxn 2 0, and 2 0, where is a vector, then [IC(A Byul S ||A4| |CCB
Proof. Let un be the i-th element of (A @ B)v, and then
04
ludl
<
||
j=l
Consequently,
2
2
m
[[C(A o B)ul|
Ui
S
Σ
Chi |uss
S
|| o0
Σ
Czi
bigg
k=1
k=1
i=1
||-A| oo ||BBlll
Similarly, A 2 0 and e 2 0, then
2
Ui
S ||
Uijyi = ||
Moreover, for matrices Baxm 2 O, and Cmxp 2 0, let Cj be the j-th column of C. and then
[I(A @ B)C[|==
Σ (I(A @ B)C,11 S
|B,,12 = || |CC|
29

---

a = =f and B = +g
g for some perturbation vectors and and define T exla03-
Then it easy to find that
+ 03 -M + - T T exla"
TT.
Let Ep = ef 0g) 1,, S = o T". Since |e" -11 < 2/2r for |x|| < 1, we have
(ET) (e>(itt))) - 2 + g3l
as long as A + 9, < 1. This can be achieved by assuming & == 2X(11 Ilo + (1gll>e) < 1, since in 1/2. Then clearly Er <.
this case Alf + g31 X + =l
Consider S, = Su + S. and S, = S. + 8, for some perturbation vectors 6, and d,, and let
VMS T + M(5u s. M)
= T + A(3u 08, - 8,
-) o 06.
V ET 2 VMS+ A(x a (T Er o
Then we have
VMS- S - 0 S | (T Er @ l'
< || + 6o + Er @ T'
< e//V /
1TT 1
S ee/V
< [ [11]]]))1]1
X(1 @
+ (1) Therefore, just need to show proper bounds for [ and 10ol Consider P. = (M @ T) and fic- = (M o T and then
<
frr (M ④ (T* Fr T T* ) nr (Er O 1r dr;
+ Er Ae = (T
ue (B == == pe Sc.
It can be easily verified that
|a @ 6 = diag(( ) O M
T M
$ | a (M O Imll < <e/1/>
(a 2 0, |T (Er ④ M O 1m11) Erlls> diag a (M o 1mll
(Er
< <El|T (a o 11)
||&el| = [(Er 3 M M ( < E 1 Ellpell.
Define b. = Me T (a_ @ (4) and b, = i ( o (r). and we have
b, 6 th d. (T + Er o T
+
= He 5, T (a 2 (p, + +,)) (Er T (,))
= Ur T 1 (a 2 (,)) + b, ( 06,) - (E, oT o (p. + 6,))
As result,
1|, bull 10 TT (o
< |all + [a) |T أ,))
06
< ellleell + + 0 (
06,1
< <el|Mell ( o 0z) +e|TT o + -T eF]
<ellmell 33= o 0
30

---

On the other hand, let b.g be the (i,j) element of the matrix 7 diag(a yT, and bry be the (i, element of T'T diag (a Yr" Then
fig
The TasTa//on =
n + (Er)ss)(1 + (Er) /
k=1
1=1
14i - 44yl =
Σ
(ET)a: )ai + -(Er)as +
k=1
n - (2e + s²)
TTsTEy/ak < 35ty
This implies that
|| - D((r - diag(a (a yr -
"dagg(a )(|
Ihe - 1ig12 s 3ell DIle- Σ f
Then by Theorem 7.2 Higham (2002). if if |el < we have
|6.| |
8a
|Hell +31 @ اx)مl + +3|DIl))
S
35a D
||s||||
where '1oop 1/o (D). Assume that < min (/6|| then with slight simplification, we have || 1 |6+1| I CoE, where
C. = 20 (I1mell + T @ 0 + + [s,l))-
On the other hand,
|16u|100 $ ||6,|| $ ||a
@ 8, + | 3 (T8, - sw)l
$ =||m|||
||a
④ ( T T s,)l
< æ||æ~||
④ (T* @ ((Er 3 )8)]
S E|lA||
(T* diag (a (EE 1 s
S <||m~|
@ (T* - c|/diag( a )7"" sell
e||,~||
['',))] - =e||diig( y7 (8. + 6,)1l
$ E||A~||
+1++)) |a (TT + diag
| (
+ (1 (1 =e|la- @ (T's
S <[1/=11 + 2C, diag(a
||a @o@@or (7*s.)ll)
Combining the results together, we get
V Slr S e V M + 2 +C)C)
where C, || + 2C,, a
( | + ||a 2 (T*
Finally, Theorem b) shows that ||g||> - ||9|| = [|3(00) - 3rll S (ie(k) addition, by 6 we have
a (k) a (g(t)), log ny log
Σ -(0ftt) M,))
Then using Lemma 6 we obtain
m
Ja.0 (b)
-a---a1-- = 17 log
a='=1
AGalM) M.s) loog
200 M.s)
-B = =I[g|lo0
implying that |l ()ax(((())) a² a=1ooo I I|g|lo S (iis(h) As a result, e = 2A([ +
(|o])>>) 4XvCCe((H)
31


Reviews:
[{"review_id": "faBiyJkO7j9", "reviewer": "Reviewer_sDPq", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper focuses on the so-called 'sharp Sinkhorn distance', which can be read as between the Wasserstein distance (without regularization) and the vanilla Sinkhorn divergence. The main contributions are two-fold: (i) the authors propose an L-BFGS algorithm to optimize the dual objective function of the Sinkhorn divergence leading to a stable transport plan (ii) an analytic form of the gradient of the shape Sinkhorn with respect to the input matrix is given. Some numerical experiments to corroborate the theoretical finds are illustrated on synthetic and real datasets.", "summary_of_the_review": "Overall, the paper lacks some theoretical results or an extensive numerical study to prove the stability of the proposed LBFGS algorithm. \n\n**Superlative phrases**\n\nIn the core of the paper, I noticed that there are lots of superlative phrases. I am wondering if this fact is scientifically acceptable. \n- Abstract: \"superiror performance\"\n- Page 8: \"enjoys superior computational efficiency\"\n- Page 8: \"huge advantage of a closed-form\"\n- Page 8: \"proposed algorithm is the most efficient one\"\n- Page 9: \"highlights the value of the proposed method\"\n\n**Typos**\n- Page 2: \"The optimal solution to (1)\" --> \"An optimal solution to (1)\". Problem (1) could have many solutions.\n- Page 2: \"satisfies some suitable conditions\" --> Add a reference.\n- Page 3: \"diagonal matix\" --> \"diagonal matrix\".\n- Page 5: in Theorem 2: the norm $\\||{\\gamma^*}\\||$ is not specified; is it the $\\ell_2$-norm? or $\\ell_\\infty$-norm? etc.\n- Page 6: in Theorem 4: the matrix $D$ is not defined.\n- Page 9: \"methods achives\" --> \"methods achieves\"\n- References: \"wasserstein\" --> \"Wasserstein\"\n- References: \"sinkhorn\" --> \"Sinkhorn\"", "strength_and_weaknesses": "**Strength**\n- LBFGS algorithm to optimize the dual objective function of the Sinkhorn divergence.\n- An analytic form of the gradient of the sharp Sinkhorn with respect to the cost function.\n\n**Weakness**\n- The authors argue that the transport plan resulting from the L-BFGS algorithm is stable compared to ones from vanilla /stabilized Sinkhorn iterations. For me, this comparison is not fair, since it is only visually; it is mandatory to establish a theoretical result or an extensive empirical study showing this gain of stability. I mean one can show that the Frobenius norm (or any matrix norm) between the LBFGS plan and the unregularized plan is less than the vanilla Sinkhorn plan and the unregularized plan. Alternatively, one can investigate the difference between the objective dual function.\n- The motivation behind calculating the gradient cost should be highlighted. \n- The numerical experiments on real data are limited.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "dZoUW1efPYn", "reviewer": "Reviewer_eynF", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper studies the problem of entropy-regularized optimal transport (EOT). The authors advocate the usage of L-BFGS algorithm to solve the EOT problem instead of the popular Sinkhorn algorithm. They analyze the convergence of L-BFGS on EOT and demonstrate its better stability compared to the Sinkhorn algorithm. They also derive a closed form of the gradient of the Sinkhorn loss w.r.t the cost matrix, leading to efficient computation of the backward pass. Finally, we illustrate their results with numerical examples.", "summary_of_the_review": "The paper studies a relevant problem and provides some new insights. However, I have some concerns about the quality and novelty of the work as detailed in my major comments in the \"Strength and Weakness\" section. I thus recommend rejection.", "strength_and_weaknesses": "**Strength:**\n1. The problem studied is a highly relevant problem.\n1. They provided a convergence analysis of the L-BFGS method on the EOT problem.\n1. The advocated L-BFGS + analytic gradient is computationally more efficient than other methods with which they compare.\n1. The paper is well-written and easy to follow.\n\n**Major comments:**\n1. The novelty of the paper is not super clear to me.\n    1. As mentioned by the authors, using L-BFGS to solve the EOT problem is a known practice. One selling point of this paper is the convergence and stability analysis of this algorithm (Theorems 2 - 4). However, this seems to be a practice problem given that L-BFGS has been extensively studied in other (potentially broader) settings. Is the proof technique new? Moreover, why can you conclude $T^{(k)}$ does not underflow from (e)? Don't you need to prove a lower bound for $min T^{(k)}$?\n    1. Cuturi et al. (2019) and Cuturi et al. (2020) consider gradients w.r.t data points. Since the cost matrix is constructed from $c(X, g_\\theta(Z))$ in the generative modeling application, knowing the gradients w.r.t data points is enough to get the gradients w.r.t $\\theta$. Can you also compare your approach with theirs?\n    1. You missed an important reference on this topic [1].\n1. The numerical experiments are not convincing. In particular,\n    1. Fig. 2 does not provide a fair comparison. As mentioned by the authors themselves, one iteration of the L-BFGS algorithm is more time-consuming than one iteration of the Sinkhorn algorithm. Why not use a similar plot as in Fig. 3?\n    1. In Table 2, what are the stopping criteria for each of the methods?\n    1. In Sections 7.2 and 7.3, while the goal is to compute the Sinkhorn loss, the performance is measured in Wasserstein distance. This mismatch makes it difficult to evaluate the improvement. I would use the Sinkhorn divergence (Feydy et al. '19) instead since this value should decrease as well.\n1. I would like to see the full algorithm of the advocated forward and backward passes with complexity analysis.\n\n[1] Tianyi Lin, Nhat Ho, Michael I. Jordan. On the Efficiency of Entropic Regularized Algorithms for Optimal Transport. JMLR (2022).\n\n**Minor comments**\n1. The notation $\\mu^{-}$ is a bit confusing. I would use $\\mu^{-1}$ instead.\n1. It would be better to define $s_v$ directly rather than $\\tilde s_v$.\n1. For matrices $M$ and $T^\\star$, their tilde versions are not defined.\n1. It would be interesting to know how the marginal constraints of $T^{(k)}$ change with the number of iterations in Theorem 3.\n1. In Theorem 4, how large is $\\sigma$ typically?", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "OvMIDZ1uK8_", "reviewer": "Reviewer_toCu", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper provides an analytic overview of the ``sharp Sinkhorn loss'' and its derivatives with respect both the input weights and the cost matrix (that typically depends on the input locations). Formally, given two discretes probability measures $\\mu = \\sum_i a_i \\delta_{x_i}$ and $\\nu = \\sum_j b_j \\delta_{y_j}$, and a cost $M = (|x_i - y_j|^2)_{ij}$ (or actually any other expression), the considered loss reads\n\n$$ S(M,a,b) = \\braket{T^\\star, M} $$\n\nwhere\n\n$$T^\\star \\in \\mathrm{argmin} \\braket{T, M} - \\eta h(T), $$\n\nwhere $\\eta > 0$ is a regularization parameter, $h$ denotes the entropy function, and the minimization is performed over the transportation polytope between $\\mu$ and $\\nu$ (that, in this discrete setting, only depends on $a = (a_1,\\dots, a_n)$ and $b = (b_1,\\dots, b_n)$). A well-known fact is that $T^\\star$ can be computed (or at least, estimated) by running the Sinkhorn algorithm, an iterative fixed-point procedure. \n\nIn applications (especially in generative models), it is often needed to compute gradients of $S$ with respect to $a,b$ and/or $M$. However, the latter is typically obtained via ``unrolling'' the iterations of Sinkhorn (via automatic differentiation). In contrast, this work computes $T^\\star$ relying on a L-BFGS approach instead of the usual iterative loop. \n\nWhile using L-BFGS to get $T^\\star$ is not entirely new, this paper has the following main contributions : \n- They study precisely the analytical behavior of the L-BFGS algorithm. The main takeaway are that convergence rates are (asymptotically) the same (i.e. geometric convergence), but the L-BFGS approach can be provably stable. This is an advantage over the standard Sinkhorn loop, and while unstability of the Sinkhorn loop is typically handled by performing iteration in log-domain, this comes at the price of computational efficiency. Here, the proposed approach gets the best of both world, being stable while remaining efficient. \n- An analytic expression of $\\nabla_M S(M, a,b)$ has been derived, allowing \"direct\" differentiation rather than unrolling the iterations when---for instance---training a generative model.\n- Showcase their approach through some numerical experiments. ", "summary_of_the_review": "I think this is a solid work overall, with potential impact for computational OT practitioners. Nonetheless, few details remain to be clarified/discussed in my opinion before fully supporting it. ", "strength_and_weaknesses": "## Strengths\n\nThis is a solid work overall. I particularly appreciate that algorithmic aspects (convergence, stability) are supported by precise theoretical results (and not just empirical ones / informal intuition). \nI also think that the main idea conveyed by the paper may be impactful when it comes to use OT in ML tasks. \n\n## Weaknesses\n\n- My main concern is that, in some sense, focusing on the loss $\\braket{T^\\star, M}$ feels somewhat ``outdated''. While it is the seminal loss proposed by Cuturi in 2013, many variants/improvements have been considered since. In particular, I would think that the so-called _Sinkhorn divergence_ [1,2,3] have taken the lead when it comes to introduce OT-based losses in ML. While I guess that a vast majority of this work holds/adapts to that setting, I think it would have been arguably better to focus on it. If there is a precise reason to focus on $\\braket{T^\\star, M}$ rather than the Sinkhorn divergence, this should be discussed in the paper. \n- One may argue that, although fairly convincing, the experimental section is not ``groundbreaking'' \n\n## Complementary comments/suggestions/questions\n\n- I think that few statements in the introduction should be made more precise. For instance, \"the computational cost of OT quickly becomes formidable\" --> give the precise complexity; similarly, it would be nice to recall what is the best theoretical complexity known for the Sinkhorn algorithm. \n- It is slightly inaccurate to talk about \"**The** optimal solution to (1)\" as it may not be unique (in contrast to those of (2), as highlighted after in the paper). \n- [minor suggestion] Of course, that is a matter of taste, but I think that using $\\eta = \\lambda^{-1}$ would alleviate notations and has become more standard in computational OT literature. That would disambiguate claims like \"Sinkhorn-type algorithms may be slow to converge, especially for small regularization parameter\" (here, one may think that \"small regularization parameter\" refers to \"small $\\lambda$\", while if my understanding is correct, this should be understood as \"small $\\lambda^{-1}$\"). \n- Formally speaking, I think that (5) is a _concave maximization problem_ (but of course this is equivalent to a convex minimization problem). \n- In Theorem 1, $s_v$ has not been defined if I am correct. \n- \"exponentiation operation are more expensive than matrix multiplications\", is that correct? I would expect term-wise exponentiation to be faster as it is $O(n^2)$ vs $O(n^3)$ overall (unless \"matrix multiplications\" refers to \"matrix-vector multiplications\", but even there I would expect the complexity to be at least similar). I also think that both operations can be efficiently performed on a GPU. \n- [question] In Theorem 1 and 4, what guarantee that $D$ is invertible (and that $\\lambda_{\\min}(D) \\neq 0$)? I tried to quickly check the appendix but did not find any discussion on that topic (I may have missed it though). \n- [suggestion/question] In Theorem 4, it may be more \"reassuring\" to state/recall that the condition $\\epsilon^{k_0} < ... $ is always met since $\\epsilon^{(k)} \\to 0$ (and is decreasing); assuming the rhs term is not $0$. Related to the above question, do we have any control on how long it takes for this condition to be satisfied? \n- In numerical experiments, it is observed that when the regularization parameter $\\lambda^{-1}$ is too large, $S$ is a poor approximation of $W$ in sense that it does not go to zero in that setting. In my opinion, this is an additional argument in favor of using the Sinkhorn divergence $\\mathrm{SkDiv}$: while for large regularization parameters these divergences will naturally differ from exact OT in general, one still has $\\mathrm{SkDiv}(\\mu,\\nu) = 0 \\Leftrightarrow \\mu = \\nu$ (and more precisely, they metricize the same topology as $W$ for any $\\lambda$), so they may be more suited to study the expressiveness/convergence of generative models. \n- [typo] In the appendix, both section B.4 and B.5 are called \"proof of Theorem 3\". \n\n## References \n\n- [1] _On wasserstein two-sample testing and related families of nonparametric tests_, Ramdas et al.\n- [2] _Learning generative models with sinkhorn divergences_, Genevay et al.\n- [3] _Interpolating between optimal transport and mmd using sinkhorn divergences_, Feydy et al.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "1MNttDvJ-J", "reviewer": "Reviewer_7wxm", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposes to use L-BFGS algorithm to solve the OT problem by exploiting a loose condition in the Sinkhorn algorithm. They derive a formula to compute the derivative with respect to cost matrix of the transportation cost. Then they prove some results of the convergence of their proposed algorithm.\n", "summary_of_the_review": "Results of Theorem 1 is useful in application such as generative modeling, but other results are not. ", "strength_and_weaknesses": "Strength:\n\n1. The paper presents some useful results (Theorem 1) of computing derivatives of the OT cost between two distribution with respect to the cost matrix. The results could be used in some important settings, i.e generative modeling.\n\n2. The author also presents some theoretical results to support their arguments. More particular, the convergence rate of the objective function, bound of the derivative etc. Those results are presented in Theorem 2,3 and 4.\n\n3. Emperical results show that the method work well and faster than other methods. \n\nWeaknesses:\n\n1. Theorem 2,3 and 4 are helpful but not useful. In particular, index $I$ in Theorem 2 depends on the solution. The rate $r$ in Theorem 3 is unknown, it is supposed to be dependent on $\\lambda$.  The same situation in Theorem 4, where $\\sigma$, $C_S, C_v$ and $C_u$ are difficult to quantified as a function of the original parameters: the marginal distribution  $a$, $b$ and the cost matrix $M$. That makes it very difficult to compare the order of the complexities of this method and the Sinkhorn algorithm in order to answer the question if this method is faster because of its initialization or it has a larger step to the optimal value etc. Note that there are already some works [1] to show the linear convergence of the objective function of the Sinkhorn algorithm, but it does not improve the upper bound of its complexity [2]. \n\n2. Empirical results are limited to some toy examples and small data set (MINIST).\n\nReferences\n\n [1] Vladimir Kostic, Saverio Salzo and Massimiliano Pontil. Batch Greenkhorn Algo- rithm for Entropic Regularized Multimarginal Optimal Transport: Linear Rate of Con- vergence and Iteration Complexity. Proceedings of the 39th International Conference on Machine Learning, PMLR 162:11529-11558, 2022.\n\n[2] Dvurechensky, P., Gasnikov, A., and Kroshnin, A. Com- putational optimal trans- port: Complexity by accelerated gradient descent is better than by Sinkhorn’s algorithm. In International conference on machine learning, pp. 1367– 1376, 2018.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject