Safe and Sparse Newton Method for Entropic-Regularized Optimal Transport
Authors: Zihao Tang
Yixuan Qiu
Venue: NeurIPS 2024 poster (2024)
Keywords: optimal transport
Newton method
sparsified Hessian matrix
global convergence
quadratic convergence rate
Primary area: optimization
TL;DR: A sparsified Newton method to solve entropic-regularized optimal transport, with low per-iteration computational cost, high numerical stability, provable global convergence, and quadratic local convergence rate.

Abstract: Computational optimal transport (OT) has received massive interests in the machine learning community, and great advances have been gained in the direction of entropic-regularized OT. The Sinkhorn algorithm, as well as its many improved versions, has become the *de facto* solution to large-scale OT problems. However, most of the existing methods behave like first-order methods, which typically require a large number of iterations to converge. More recently, Newton-type methods using sparsified Hessian matrices have demonstrated promising results on OT computation, but there still remain a lot of unresolved open questions. In this article, we make major new progresses towards this direction: first, we propose a novel Hessian sparsification scheme that promises a strict control of the approximation error; second, based on this sparsification scheme, we develop a *safe* Newton-type method that is guaranteed to avoid singularity in computing the search directions; third, the developed algorithm has a clear implementation for practical use, avoiding most hyperparameter tuning; and remarkably, we provide rigorous global and local convergence analysis of the proposed algorithm, which is lacking in the prior literature. Various numerical experiments are conducted to demonstrate the effectiveness of the proposed algorithm in solving large-scale OT problems.

Safe and Sparse Newton Method for Entropic-Regularized Optimal Transport
Zihao Tang
School of Statistics and Data Science
Shanghai University of Finance and Economics tangzihaodstu, sufe. .du.eedue.. cn
Yixuan Qiu"
School Statistics and Data Science Shanghai University of Finance and Economics giuyixuan@sute edu. cn
Abstract
Computational optimal transport (OT) has received massive interests in the machine learning community, and great advances have been gained the direction of entropic-regulariied OT. The Sinkhorn algorithm, as well as its many improved versions, has become the de facto solution to large-scale OT problems. However, most of the existing methods behave like first-order methods, which typically require a large number of iterations to converge. More recently. Newton-type methods using sparsified Hessian matrices have demonstrated promising results on OT computation, but there still remain lot of unresolved open questions. this article, we make major new progresses towards this direction: first, we propose novel Hessian sparsification scheme that promises a strict control of the approximation error; second, based on this sparsification scheme, we develop a safe Newton-type method that is guaranteed to avoid singularity in computing the search directions; third, the developed algorithm has a clear implementation for practical use, avoiding most hyperparameter tuning; and remarkably, provide rigorous global and local convergence analysis of the proposed algorithm, which is lacking in the prior literature. Various numerical experiments conducted demonstrate the effectiveness of the proposed algorithm in solving large-scale OT problems.
Introduction
In recent years, optimal transport (OT, (36]) has received massive attentions from the deep learning community, and has become fundamental modeling tool in modern statistical machine learning 35 23). One major challenge of applying OT to real-life problems, the computation of large-scale OT, has also gained many progresses the direction of approximate OT methods, especially the entropic-regularized OT [6].
Consider two discrete probability measures Σ u,ō.; and Σド= bydy,- where (01,...,a,)) AT and (b.....,b,.) bm.) are two vectors satisfying Σ a; プm by 1, a; > 0. b, > 0, L,.., 1,.... m. { and {y, } are two sets of data points, and d is the Dirac measure position æ. Without loss of generality we assume that n 2 m, their roles can be exchanged. Let M (M,y) mx be cost matrix, whose entry Msj typically represents the
"Conrespondigg author.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

distance between data points and yj. Also define II( b) = {T E Ruxm Tim a, 1", T' 2 0}, where the inequality sign applies elementwisely. Then OT between the two measures j and D can be characterized by the following linear programming problem,
min P)ll(a
(P,M),
where
Assuming n m. the computational complexity of solving (I) using
standard linear programming solvers is typically at the order of O( log(n)) [2 which can be difficult even for moderate and m. One approach approximating the solution is to add an entropic
penalty term to the objective function, leading to the entropic-regularize OT problem G
(T. (T(T), Tell(a,b)
min
(2)
where ((T) Ei Tiy Log Tip) is the entropy term, and v > is regularization
parameter. The objective function in (2) is ry-strongly convex on [T b), so 2 has a unique global
solution, denoted as T The matrix T often referred to as the Sinkhorn transport plan. Problem ( can be solved by the celebrated Sinkhorn algorithm multiplication operations. Due to its computational advantage, the Sinkhorn algorithm has become
38 (3], which is based on efficient matrix vecctor
the de facto solution to large-scale OT problems for long time.
However, from a practical point of view, the Sinkhorn algorithm typically demonstrates a sub-linear convergence speed, thus requiring a large number of iterations for a moderately high precision. More recently, second-orde methods, such as the Newton method, have attracted growing number of researchers rethink the solution to (Proposition 4.4 of (30)::
34) It has been shown that @ has a dual problem
max L(a, B) max a²a a + b 1 Σ Σ e" P(a.t)) a a RER", (3) 0,8
a.f
Let a = a.,....,,,"""" a and 8 = ........ )m)) be one optimal solution to (3) and then the Sinkhorn transport plan T be recovered as T 7(a", _,,), where for two vectors a (ai...., ax))" and (BI...., Bm)) T ra, (a,9 3) is a matrix with entries
[r(a,
exp (1 (a; + 33 -,;))-
Remarkably, 3) is equivalent to a smooth and unconstrained convex optimization problem, so it can be solved using first-order methods such as gradient descent, or second- methods including the Newton method. It is well known that under some smoothness assumptions on the objective function, the Newton method has a fast quadratic local convergence rate, and hence typically requires much fewer iterations than the Sinkhorn algorithm. However, each iteration in the Newton method involves solving dense linear system, resulting in a per-iteration cost of O which is nearly of the same order as that of unregularized OT problem. To this end, [34] proposes the Sinkhorn Sinkhorn-Newton Sparse (SNS) algorithm, which approximates the Hcssian matrix of 33 by a sparse one, so that each iteration roughly has a computational cost of O(r²), the same order the Sinkhorn algorithm.
The SNS algorithm demonstrates promising results for solving the entropic-regulariin OT, yet there remain a number of open questions. First, the SNS algorithm relies on initial Sinkhorn iterations to achieve moderately small optimality gap, but in practice, the cost for this warm initialization also needs to be taken into account. Second, though SNS advocates using the conjugate gradient (CG) method to solve the sparse linear systems, there no strong guarantee on the invertibility of the sparsified Hessian matrix. Third, (34] conjectures that SNS has super-linear local convergence rate, but it is not yet proven by thcoretical analysis.
In this article, we make major progresses on the second-orded method for entropic-regularized OT by resolving the three challenges above. First, we carefully design new sparsification algorithm that has a well-controlled approximation error, and then propose a safe Newton-type algorithm, in the sense that the linear systems for computing the search directions are always positive definite. This property addresses the invertibility issues of existing sparsified Newton methods, and is crucial for practical implementation. We also show that the proposed algorithm has global convergence guarantee, so no initial Sinkhorn iterations are needed. Most importantly, we provide solid theoretical analysis of the proposed method, and show that achieves quadratic local convergence rate similar to the vanilla Newton method.

---

Contribution Our main contribution compared to prior art is summarized as follows:
1. We propose novel Hessian sparsification scheme that promises strict control of the approxi- mation error.
2. Based on this sparsification scheme, we prove that the sparsified Hessian matrix is always positive definite regardless of the sparsification parameter. This property then leads safe Newton-type method that is guaranteed to avoid singularity in computing the search directions.
3. The developed algorithm has a clear implementation for practical use, avoiding most hyperpa- rameter tuning. An efficient implementation included in the RegOT Python package]
4. We provide rigorous global and local convergence analysis of the proposed algorithm, which is lacking in prior literature.
2 Background and Related Work
Notation Throughout this article, we use f(x), g(x), and (x)) to represent the objective function, gradient, and Hessian matrix at point æ, respectively. For a matrix A, let AG. be the vector of the i-th row of A, and A. be the vector of the j-th column of A. For a vector e (uI...., )T let denote the vector (UI...., en-1)", and for a matrix A A= = (A......, A.). let represent the matrix (A....... A,,nn.). The symbol || stands for the Euclidean norm for vector argument, and represents the operator norm for matrices. The fn-mrrm and infinity norm are denote by || || and ||. oc: respectively, applying both vectors and matrices.
The main objective As has been explained in Section [ the key to the computation of entropic- regularized OT is to solve its smooth dual problem ( which is equivalent to the convex optimization problem min... - L(a. 3). However, it is worth noting that the variables (a, 3) have one redundant degree of freedom: if (a",3* is solution to 33. then so is (a" + cla, 8 clm) for any Therefore, to ensure identifiability, we globally set Bm = 0 without loss of generality. In what follows, we use the vector x = (a yT E IR=+=-1 to collect all free dual variables. Clearly, we )T always have 3 = , = (े 0)T, so 3 and /3 will be used interchangeably when we consider functions of 3. Then the main target of this article is to efficiently solve the optimization problem
min f(x) = C(a, 3) = 71 T(, 3)1,. - a - 3Tb.
is known that (x) is a strictly convex function, so if (4) has solution, then it is unique. has also been proven in the existing literature that the gradient and Hessian matrix of f(x) have the following closed-form expressions (see for example (21):
V Val B)] =
T1m a
g(x) =
V (()
TT
T T=,
(5)
VoL(a, 3)
V (V (VaL(a,)))
[diag [diag[T1m) (1'm)
T diag(TT1) 1.))
H(x)
-
[ (V ((,((a,))) B))]
C(, 3)
T
Computational OT There are a huge number of methods developed to solve and (4). The 33 o can be interpreted as applying the block coordinate descent (BCD),
Sinkhorn algorithm [38,
also known as the alternating minimization (A.) method, to (4). Along this direction, many extensions of the Sinkhorn algorithm have been proposed. For example, [13] develops an accelerated version of AM, and a greedy variant of the Sinkhorn algorithm, named Greenkhorn, is developed in
M. Accelerated first-order methods, such as the adaptive primal-dual accelerated gradient descent (APDAGD, (00) and adaptive primal-dual accelerated mirror descent (APDAMD, (19)
) have been proposed to solve the constrained problem (2). Moreover, since the dual problem 4 is smooth, various quasi-Newton methods, such as the Limited-memon Broyden-FFetccher Goldfarrb--od-bbbboooooo Saadnoo (L-BFGS) method 20), have also been used 77]. Second-order methods are less common in the literature to solve 44. mainly due to the high cost of solving the linear systems, and (34] proposes practical sparsified Newton method to approximate the true Hessian by a sparse matrix. Another direction to utilize the sparsity is the importance sparsification method [18], which replaces the dense kernel matrix in the Sinkhorn algorithm by a sparsified one through random sampling.
htttp::////e thub. www. com/yixuan/regot python
3

---

Newton-type methods The Newton method is classical second-order optimization algorithm that has been extensively studied. Assuming f is twice continuously differentiable, the Newton method generates sequence of iterates { using the updating formula 3k++1
(*x)) (=*)))=~9(aa ).
However, computing the linear system (Tx,) ) g'g()))) with a dense Hessian matrix generally costs O(n²) of computation, so the sparsification of H (axx) with solid convergence guarantee is the main focus of this article. The Newton method also has many extensions, for example the regularized Newton method to solve problems with non-isolated solutions (17), and the inexact Newton method 8 that allows inexact solutions to the linear system [H (x+)] (xxx). Our proposed method use a combination of these ideas to reduce the computational cost of the classical Newton method.
OT in machine learning OT is a blooming research topic in modern machine learning research. As a powerful tool characterize the transformation of statistical distributions, OT has wide applications in deep generative models 22 14] domain adaptation and transfer learning [5] 4] and fairness in machine learning
12 25 32 among many others. Readers are referred to review articles such as for a summary of machine learning tasks and methods that utilize OT.
35
23
3
The Proposed Algorithm
To complement existing second-order methods for entropic-regularized OT, in this article we propose the SSNS algorithm, short for safe and sparse Newton method for Sinkhorn-type OT. It has two central components: new scheme to sparsify the Hessian matrix, and a Newton-type algorithm to update iterates.
3.1 Sparsifying the Hessian matrix
Ideally, the sparsified Hessian matrix, denoted by Ha(x), should meet two criteria: first, it should be close to the true (x) with a tunable approximation error; second, it needs to preserve the positive definiteness of the original Hessian matrix. The first point to ensure that Hz(x) does not lose too much information of the true (x). The second point is especially important in implementing the Newton method, since a linear system [Hs(<)]d
g(x) needs to be solved to compute the search direction d, and one needs to make sure that H3 is invertible. These two points are briefly mentioned in 33 as an implementation practice, but without a strong theoretical guarantee.
To this end, we introduce a new adaptive method to sparsify Hx),, which takes the two points above as first principles in both theoretical analysis and practical implementation. To describe the whole algorithm, first define an operator select_ small(m. (v. o), which takes a vector v = (01,-.., )² and scalar 8 > 0 as inputs, and outputs mask vector Ø = (#1,-... Rn)"" in the following way. Suppose that (1) I (2) S S Us(n)) are the sorted values of c. where T(s) is the original index of the s-th smallest element in v. We also define π 7-) to be the inverse of (-), i.e., IT T is the integer such that π (- (i)) i. Let S be the largest integer such that ΣS Us(s) 8. Then we set s 1 π -i) < S, and Φ; otherwise. For example, if = (2, 1, 5, 2)² and 5 then have 6 = (1, 1, 0, 0, 1)",, since 1 + 2 + 2 < 8 but 1 + 2 + 2 + 3 > 5.
(U1 (Vi..... )² and a mask é = ( .r))² and outputs a vector u = (u....., ) with u, = Φ; u, For the example above, we have apply (y Φ) = (2, 1, 0, 0, 2) Then given the current dual variables r = (a" 3r Algorithm outputs a sparse matrix HA as an approximation the true Hessian matrix H (xr)
Next, define an operator apply mask(u, Φ), which inputs vector
It is worth noting that in the construction of Hs, the diagonal elements dTag(Tmm)) and diag (TT) 1,) use the original 7 matrix, whereas the off-diagonal elements are the sparsified Ta matrix, with the last column removed. In this way, Theorem 2 shows that Hş indeed has a well-controlled approximation error as requested by our first criterion.
Theorem (Approximation error of sparsification). Let H5 be output by Algorithm 7 with a given vector x, and define D H(x) H3. Then for any 5 2 O. we have D 2 0, where the inequality sign holds elementwisely, and for each i.j = 2......, m
|D+11 I 17
[D =
Moreover. ||a|| S -88

---

Algorithm Sparsifying the Hessian matrix.
Input: Dual variable vector z = (a".,"") 3"""", threshold parameter ò 2
Output: Sparsified Hessian matrix Hg
I: Initialize a zero matrix E _ and compute T =(a,
do 3: Φ f select (T,
2: for j = 1,2.....
A t apply_pppaaee (T.).
4: for == 1........" n do
5: Φ 1 select_smalle (A,, 5),
A f apply mmply_mask// A (),
6: Ta T -
diag (T T5
Tg diag(77 ( 1n)
7: Hs + 1
Theorem T shows that all the row-wise and column-wise 4--norms of the residual matrix D is bounded by the pre-specified threshold 6, up to constant multiplier 17 More importantly, the overall operator norm of D also has the same upper bound, and this property plays key role in analyzing the convergence behavior of the SSNS algorithm later in Section 3.2
Next, Theorem 2 validates our second criterion Hg: the sparsified Hessian matrix by Algorithm T is guaranteed to be positive definite, regardless of the threshold parameter 6. Therefore, it can be safely used to compute the Newton search directions.
Theorem 2 (Positive definiteness). Suppose that n 2 m. For any a E R", B E R= and M E R=xm, let T T(a, T1mm and = T" 1,, Then for any 5 2 0. the matrix Hş output
by Algorithm 7 is always positive definite. Specifically, we have
+
mmax (Hs) 2n max maax(1rlles [|rl|l> [cc]bo0h
Tmin (Hi) 21
min Tajs
2n
[.)
where muax ()) and m-in(-) represent the largest and smallest eigenvalues of symmetric matrix, respectively.
The two features presented by Theorem D and Theorem 2 are crucial designing the optimization algorithm and analyzing its theoretical properties, which we introduce in details in Section 3.2
3.2 The SSNS algorithm
In Theorem 2,
we have shown that the sparsified Hessian matrix H& always positive definite, thus invertible. However, in actual computation, Hls may be nearly singular, thus causing numerical instability. To further stabilize the optimization procedure, we propose safe and sparse Newton method inspired by the regularized Newton method [17 and the Levenberg-Marquar algorithm
16 22
One of the key ingredients of SSNS is to introduce shift parameter À in solving the Newton linear system, leading to a search direction of the form p = -HH + Al) (4)-'g(aa)) "((x)) in each iteration. The shift parameter À has multiple functions. First, it stabilizes the vanilla Newton linear system H1, 115'gg( "9(x), as the matrix H₃ + AI has smaller condition number than H3. which potentially accelerates iterative linear solves such as the CG method. Second, if we allow X to vary along iterations, has the effect of adjusting the search direction when - H "=')(a)) is not a good candidate. Third, show in the convergence analysis, an adaptive À plays an important role in both the global and local convergence of the algorithm.
will later
On the other side, the sparsification threshold 8 is also essential to the design of the algorithm, as we need to control the approximation error of H,s during the iterations. Intuitively, when the current iterate Iy is close to the optimum, Hş should be also close to the true Hessian matrix (xxx) to achieve the fast convergence rate, whereas when I& is far away, it may be beneficial use large 5, since Hg is then more sparse and leads to a faster computation
the linear system -(H, + AI) g(x).
Overall, in SSNS we allow 5 and À to vary in each iteration, and use the notation d5 and X to reflect this adaptivity. Intuitively, both òg and d cause approximation errors to the true Hessian matrix in determining the Newton search direction, they should tend to zeto when zg is sufficiently close to
5

---

the optimum. On the other hand, we do not want ō5 and d to decay too fast, since otherwise they make Hş too dense and the linear system may be ill-conditioned. To this end. we their values according to the current gradient norm [[g( Il, which can be viewed as an indicator of the distance between FyK and the optimum z². Specifically, we set A = Us (++]] and 5y = Vo |9 where μi is a parameter that dynamically changes in each iteration, and Vo > 0 and y 2 1 are two constants.
These settings are justified by the convergence analysis later in Theorems 3 and
Before presenting the full details, we first introduce a few functions and quantities that play important roles in the algorithm. In each iteration, SSNS computes quantity Pk.
f(xx) f( + CAPR)
Pk
my (0) mx(EpPN)
where & > 0 is step size parameter, and the m²(-) function is local quadratic approximation to the objective function:
m& (p) f(xx) + [g([))))>>> --- - Ha.p.
(6)
The quantity Ak has two functions: first, it adjusts the th parameter, and second, also determines whether one should accept the proposed new iterate 24 + E&PE. A negative value of Pk: means that the objective function actually increases at the new iterate, so one should reject the proposed point, and increase A in the next iteration for potentially better search direction. The overall SSNS algorithm then given in Algorithm 2
Algorithm 2 Safe and sparse Newton method for Sinkhorn-type optimal transport.
Input: Initial point Fo. parameters {No, Vo, Ct2 Gus K3 > 0, r 2 1, Po E (0, 1)) Etot > 0 Default values: fuu = 1, Vu = 0.01, q = 0.1, ca = 1, K = 0.001, 7 = 1, Pu =
Output: 2%
1: for k == 0, 2, do
2: Compute fk g(rx), 6% vollgu|1
3: if ||9x|| < Eod then
4: return F%
5:
Compute Hš, according to Algorithm D
with x + Tk
6:
Compute pk (1Hôr p&
gh
7:
Select any & E [q,, Col
f(zk) f(zk + ENPk)
8:
Compute PI: (0) (m(PpP))
m+(-)) is defined (6)
if Pk: < Pa
9:
Update p6++1
max(uu/2, k]]
if px Po otherwise
14::
10:
if Pk: > 0 then
11:
Tk++1 = Ig + E1Pk
12:
else
13:
Tk+1 = Tk
3.3
Step size selection
It is worth noting that in Algorithm the step size εx can be taken arbitrary value from fixed interval [q, Gu]- This is quite different from typical line search methods, which require the new iterate I ExPE to satisfy certain conditions, such the sufficient decrease condition and the curvature condition [24). In other words, it is one of the advantages of SSNS that the cost of step size selection is predictable: we can consider fixed number of candidates, and select one of them according to some criterion. At one extreme, it is completely acceptable to always set ÇA 1. However, practice, we advocate Algorithm 3 heuristically choose εk. which considers a fixed number of candidates, and computes the objective function value for each candidate. The algorithm will early stop if decrease in function value found. If all candidates increase the function value, then return the step that results in the smallest function value. Based on our empirical results, setting (Εjo)- Eny E1a)) = (1, 0.25, 0.1) gives reasonably good performance in most numerical experiments.

---

Algorithm 3 A practical step size selection method for Algorithm 2
Input: Candidate step sizes Initialize Ç =
S10] > E11)
SiN 0, current 2i,, Pk, objective function f(()
+000
2: for i i= 0, 1,.,,.... N do
3: Compute Tirial = 2k + PH Sirial ((trial)
4: if fuial f then
5: <= Eti. f = furil
6: if < f(x)) then
7: return 5"
8:
return {"
4 Theoretical Guarantees on Convergence
In this section, present major theoretical contributions to the sparsified Newton method for entropie-regularized OT. We first show that starting from an arbitrary initial point Tg. the iterates generated by Algorithm eventually Algoritmmeeveenuull converge to the unique global optimum x². Therefore, SSNS does not require a warm initialization using the Sinkhorn algorithm as [4].
Theorem (Global convergence guarantee). Let { be generated by Algorithm 2 and is an optimal point of (4) Then either Algorithm terminates in finite iterations, I& satisfies
lim [[9( O, lim || x x=11=0.
The next theorem characterizes the behavior of SSNS when the iterates sufficiently close the optimum. Since Algorithm 2 the proposed new iterate 3y + ExPx is rejected Pk S 0, it is important to show that such rejections will be very rare the algorithm proceeds.
Theorem 4. There exists an integer K such that for all k 2 K, pk 2 Po. ck+11 K, and Ik+1 = Ik + EkPK-
Theorem indicates that eventually the parameter 14 upper bounded by the pre-specified value s, making the shift parameter of the Newton method, namely, A4 119111 tend to zero. This is consistent with the intuition we have described in Section 3.2 More importantly, the quantity Ps. will also be greater than the threshold 1 Po. so the new step proposed by the Newton direction will always be accepted, regardless of the bounds e cu. This makes SSNS quite flexible in picking the step sizes, unlike classical line search algorithms that pose various conditions that need to be satisfied.
Finally, Theorem 5 shows that SSNS indeed has quadratic local convergence rate that matches the Newton method based on genuine and dense Hessian matrices.
Theorem 5 (Quadratic local convergence rate). Fix & = 1. Then there exists an integer K' >> O and a constant L 0 such that for all K',
| x11l S L| Fr x²1
Overall, Theorems to provide solid theoretical guarantees the proposed SSNS algorithm, which are a major complement to the existing literature on solving entropic-reguulriieed OT using realistic second-order methods. More importantly, they validate that SSNS is indeed "safe" algorithm that is quite robust to the choice of initial value, step sizes, strength of sparsification, and hyperparameters.
5 Numerical Experiments
In this section, we test the performance the proposed SSNS algorithm on various numerical experiments]] There are huge number of algorithms developed for entropic-regullrizee OT, and focus on representative ones from each category of optimization methods: 1. the Sinkhorn algorithm the default option for entropic-regulariized OT, interpreted
BCD method; 2. the APDAGI
'The programming code to reproduce the results is available at hithtts: :/ thub. hht..cccom/aagg///aa199977//
SSNS
7

---

algorithm for accelerated first-order method; 3. the L-BFGS algorithm for quasi-Newton method; 4.
the globalized Newton method with line search; 5. the proposed SSNS algorithm.
We first consider three benchmark datasets to define the OT problem, each with two ways of constructing the cost matrices. More explanations of the experiment setting are given in Section A.2
JMNIS:: OT between pair of images from the MNIST [15] or Fashion-MNIST [37 dataset. The a and b vectors are flattened and normalized pixel values, and the cost matrix holds the fi-distances or squared Euclidean distances between individual pixels. The problem size is
n = m = 784.
ImageNet: OT between two categories of images from the ImageNet dataset [9]. We use a subset of ImageNet from the Imagenette Github repository which contains ten classes of ImageNet images. Approximately 1000 images per category are selected. We map each image to a 30-dimensional feature vector by first passing the image to a ResNet18 88 network, resulting in a 512-dimensional vector, then followed by a dimension reduction by principal component analysis. Let 2; R300 be the feature vector of an image in the first category, i ...... n, and Vj E I 30 be the feature vector of an image in the second category, J = -I...., m. Then a = n 1,, b = m "m, and the cost matrix is Mi = I|z ys la or My = || - yy 1l? The problem size is n = m R 1000.
To make the regularization parameter n, comparable in different settings, we normalize all cost matrices to have unit infinity norms, namely, M M| Then we consider settings of the regularization parameter, = 0.01 and 0.001.
For entropic-regullaizee OT, commonly-used criterion to evaluate optimality is the marginal error of the estimated transport plan. Let Iı (af, )²) be the current iterate for dual variables,
and set Ta T(aks 3a)) Then the marginal error is given by / VI al1² + -17 1n bl(2.
Coincidentally, this is exactly the current gradient norm [g(x)))l as indicated by
5 Figure I shows
the plots of marginal errors against iteration number or run time for different algorithms on the three benchmark datasets.
VRIST (101-2, 1D0 43400)
PssticnMWi15T (D2-22, (D2-14009)
log (N00Ngg Erree
Marginal
b
Aroou Nevie
Affono
-8
50
102
155
200
250
300
50
150
350
200
250
300
50
100
150
200
250
Boration
HorationNumm
Ticration
VRIST (102=2, 102 +34000)
PashicnMhiST 1D2 148981
Image/ett TTech Ceretto playert
Margin tron
Ertrr
BC3
Log30
-20
Mun 2me Secconosl
Run Smm Geconisl
Run Goconds)
Figure 1: Top: Marginal error vs. iteration number for different algorithms on three datasets. Bottom: Marginal error vs. run time. The cost matrix based on the f disttttae,, and 17 0.01.
From the top row of Figure I it is clear that second-order methods have much faster convergence speed compared first-order and quasi-Newton methods. However, for the vanilla Newton method, this advantage is weakened by its high per-iteration cost, resulting in less competitive run time performance shown in the bottom row of Figure [ The SSNS algorithm avoids this issue by using sparse matrix operations, and hence for the run time results, it still shows an order of magnitude speedup for MNIST and Mashion-MNIST data. The run time for SSNS relatively longer in the ImageNet dataset, mostly because the transport plan is more dense under 7 0.01. We show in
8

---

Figure that as becomes smaller, the Hessian matrix can be better approximated by sparse matrix, thus enlarging the advantage of SSNS. This point consistent with Theorem of [44, which states that smaller in general results in better sparse approximation. More discussions on the impact of regularization parameter is given in Section A.3
MNIST (IB]= I02 -54600)
PeshicoMRiST ODI - 505588
Imapeeeer (Dhrch Consette payer
FL
Log Margon mmox
Log30 Kurgin
Noviie
50
103
155
200
300
3.00
200
250
300
150
150
250
Beretion
VNIST (103=2, 102 (54660)
ReshioeWRSST =2. 102 546-900
Cassette
1og Margha Emox
Logt Emox
44
B20 Pool
Run Sm (seconce)
Run Siee (encente1
Run time (secends)
Figure 2: Top: Marginal error vs. iteration number for different algorithms on three datasets. Bottom:
Marginal error vs. run time. The cost matrix based on the é didistance, and 17 0.001.
We also note that with 17 0.001, the Newton method encounters numerical issues in solving the linear systems: the Hessian matrix is nearly singular, and the line search procedure fails along the computed search direction. Therefore, the Newton method does not generate any valid results in Figure 2 This issue further validates our intuition in Section 3.2 when H (æ) is nearly singular, ([H(x) g(xr) may not be a good search direction, and the introduction of the shift parameter À in SSNS may give a better candidate -(H + 9f)=+
In the ImageNet experiments above, map each image to 30-dimensional feature vector to compute the cost matrix. To study the impact of the feature dimension, in Scction [A.4, we conduct additional experiments with d 60, 90, 200, 300, 500, under the same setting in Figure 2
Next, we show in Figure 3 the results for cost matrices based on squared Euclidean distances. The implications are similar: the Newton method fails in the ImageNet dataset, and costs too much run time in the other two datasets. In contrast, SSNS shows great computational advantage. The experiment results for more test cases are given in Section An5 the appendix.
Finally, study the scalability of SSNS, we consider the following synthetic OT problem that can generate data with arbitrary dimensions.
Lcrge-scale synthetic data: The basic setting is to approximate the OT between two continuous distributions: the source is an exponential distribution with mean one, and the target is a normal
mixture distribution 0.2 0.2) + 0.8 (3, We discretize the problem the following way: let x, = 5(i - 1)/(rz (m 1), i = 1,
and yj = 5(j 1m - 1), j = 1...., m, which are equally-spaced points on [0, 5J. Define the cost matrix as My = (x; y )² Let fi and f2 be the density functions of the source and target distributions, respectively. Then we set ā, fi.)). b, fa(u,). = / (> k= āx), and bj /,,/ (E by:
Similar to the previous experiment setting, we normalize the cost matrix, M + M/IJM loc. and set 17 0.001. We then solve the problem at the scales of n m 1000, 5000, and 10000, but only considering BCD and SSNS, as other methods are too tome-consuming to proceed. The results are visualized in Figure 4 whose pattern is clear: BCD demonstrates a linear-like convergence rate, and SSNS has fast convergence speed consistent with the theoretical quadratic rate. Thanks to the Hessian sparsification, SSNS does not suffer from a high per-iteration cost, so overall it provides an efficient solver for entropic-regularriee OT even on very large problems.

---

VAIST 1122-2 -34000)
PashicoMeiST (02-2, (02-54890)
Image (Tanch Cassetie player
Lo5 Marg Troe
Froo
Amouch
50
100
150
200
250
200
30
200
350
200
250
300
50
120
350
200
250
300
MNIST (101m22 102 (5466)
Pishico/RIST (10 (02-506991
Imagetetr (Tench vs Cassette payyr
15010 to tror
Kuro
6
4
ATTACD blouiie
ATRAAD benci
0.5
1.5
2.5
2,0
3.0
2.5
1.0
2.0
2.5
1.0 2.5
2.5
Aun 5mee (SRCOnES)
Run time joscontst
time jseconde
Figure 3: Top: Marginal error vs. iteration number for different algorithms on three datasets. Bottom: Marginal error vs. run time. The cost matrix is based the squared Euclidean distance, and 77 0.001.
Syrthetle 1000, 10000
Synthetic data One 2220
Syrthhtic 10000 5=100001
con
500
Log10 Maggra treee
Lagio
130
2000
700
420
500
100
200
300
400
500
0
100
200
200
410
500
Symthetr dats on 5000, 50300
Syrthetic daca (mo 10000, tull (0000)
4.
(05 Mavgral Erors
4
Log
25
50
100
125
150
175
130
200
200
400
500
700
Run (eecends]
jeecende)
Figure 4: Comparing BCD and SSNS on large OT problems. Top: Marginal error vs. iteration number for different problem sizes. Bottom: Marginal error vs. rum time.
Conclusion
In this article, we have carefully designed and analyzed variant of the sparsified Newton method for solving large-scale entropic-regularized OT problems. As demonstrated by [34], the sparsified Newton method typically enjoys visibly faster convergence speed than first-order methods, and meanwhile maintains a comparable computational cost per iteration. This paper substantially complements the existing literature sparsified Newton method the following aspects. First, new Hessian sparsification method is developed, which nicely interacts with the convergence theory. Second, based on this scheme, we prove that the sparsified Hessian matrix is always positive definite, thus promising a safe Newton direction computation. Third, we have a thorough theoretical analysis of the proposed SSNS algorithm. which affirmatively answers a previous conjecture made in (34] on the convergence speed. Finally, our algorithm features an easy and clear implementation that does not require warm initialization or heavy hyperparameter tuning, thus being user-friendly for practical use.
L Limitations A potential limitation of the proposed method is that for some specific OT problems, the transport plan may not be well approximated by a sparse matrix. In such cases, the Hessian matrix may be too dense even after the adaptive sparsification step.
10

---

Acknowledgements
Yixuan Qiu's work was supported in part by National Natural Science Foundation of China (12101389), Shanghai Pujiang Program (21PJC056), MOE Project of Key Research Institute of Humanities and Social Sciences (22JJD110001), and Shanghai Research Center for Data Science and Decision Technology.
References
[1] Altschuler, J., Niles-Weed, J., and Rigollet, P. (2017). Nenr-linear time approximation algorithms for optimal transport via Sinkhorn iteration. In Advances Neural Information Processing Systems, volume 30.
[2] Arjovsky, M., Chintala, S., and Bottou, L. (2017). Wasserstein generative adversarial networks. Proceedings of the 34th International Conference on Machine Learning, pages 214-223.
[3] Brauer, C., Clason, C., Lorenz, D., and Wirth, B. (2017). Sinkhorn-Newton method for entropic optimal transport. arXiv preprint arXiv:1710.06633.
[4] Courty, N., Flamary, R., Habrard, A., and Rakotomamoniy. A. (2017a). Joint distribution optimal transportation for domain adaptation. In Advances in Neural Information Processing Systems, volume 30.
[5] Courty, N., Flamary, R., Tuia, D., and Rakotomamonjy, A. (2017b). Optimal transport for domain adaptation. IEEE Transactions Pattern Analysis and Machine Intelligence, (853))185-18888
[6] Cuturi, M. (2013). Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in Neural Information Processing Systems, volume 26.
[7] Cuturi, M. and Peyré, G. (2018). Semidual regularized optimal transport. SIAM Review.
[8] Dembo, R. S., Eisenstat, S. C., and Steihaug, T. (1982). Inexact Newton methods. SIAM Journal on Numerical Analysis, 19(2):400-408.
[9] Deng, J., Dong, W., Socher, R., Li, L.JJ,, Li, K., and Li, F.F.. (2009). ImageNet: A large- scale hierarchical image database. In 2009 IEEE Conference Computer Vision and Pattern Recognition, pages 248-255.
[10] Dvurechensky, P., Gasnikov, A., and Kroshnin, A. (2018). Computational optimal transport: Complexity by accelerated gradient descent is better than by Sinkhorn's 's algorithm. In Proceedings of the 35th International Conference on Machine Learning, pages 1367-1376.
[11] Genevay, A., Peyre, G., and Cuturi, M. (2018). Learning generative models with Sinkhorn divergences. In Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, pages (608-1611
[12] Gordaliza, Barrio, E. D., Fabrice, G., and Loubes, J.-M. (2019). Obtaining fairness using optimal transport theory. In Proceedings of the 36th International Conference on Machine Learning, pages 2357-2365.
[13] Guminov, S., Dvurechensky, P., Tupitsa, N., and Gasnikov, (2021). On combination of alternating minimization and Nesterov's momentum. In Proceedings of the 38th International Conference on Machine Learning, pages 3886-3898.
[14] Huynh, V., Phung, D., and Zhao, H. (2021). Optimal transport for deep generative models: State of the art and research challenges. In Twenty-Nitth International Joint Conference Artificial Intelligence, volume 5, pages 4450-
[15] Lecun, Y., Bottou, L., Bengio,
Y., and Haffner, P. (1998). Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):278-22244
[16] Levenberg, K. (1944). method for the solution of certain non-linear problems in least squares. Quarterly of Applied Mathematics, 2(2):164-168

---

[17] Li, D.-.., Fukushima, M., Qi, L., and Yamashita, N. (2004). Regularized Newton methods for convex minimization problems with singular solutions. Computational Optimization and Applications, 28(2):131-147.
(18) Li, M., Yu, J., Li, T., and Meng, C. (2023). Importance sparsification for Sinkhorn algorithm. Journal of Machine Learning Research, 24(247):1-44.
[19] Lin, T., Ho, N., and Jordan, M. I. (2022). On the efficiency of entropic regularized algorithms for optimal transport. Journal of Machine Learning Research, 23(137):1-42.
[20] Liu, D. C. and Nocedal, J. (1989). On the limited memory BFGS method for large scale optimization. Mathematical Programmingg 45(1)5503-528.
[21] Luise, G., Rudi, A., Pontil, M., and Ciliberto, C. (2018). Differential properties of Sinkhorn ap- proximation for learning with Wasserstein distance. In Advances in Neural Information Processing Systems, volume 31.
[22] Marquardt, D. W. (1963). An algorithm for least- squares estimation of nonlinear parameters. Journal of the Society for Industrial and Applied Mathmmatics, 11(2):33-441.
[23] Montesuma, E. F., Mboula, F. N., and Souloumiac, A. (2023). Recent advances in optimal transport for machine learning, arXiv:2306.16156
[24] Nocedal, J. and Wright, J. J. (2006). Numerical Optimization. Springer Series in Operation Research and Financial Engineering. Springer, New York, NY, second edition.
[25] Oneto, L., Donini, M., Luise, G., Ciliberto, C., Maurer, A., and Pontil, M. (2020). Exploiting MMD and Sinkhorn divergences for fair and transferable representation learning. In Advances in Neural Information Processing Systems, volume 33.
[26] Ortega, M. and Rheinboldt, W. C. (1970). Iterative Solution of Nonlinear Equations in Several Variables. SIAM.
[27] Papadakis, N. (2015). Optimal Transport for Image Processing. Thesis, Université de Bordeaux Habilitation thesis.
[28] Pasechnyuk, D. A., Persiianov, M., Dvurechensky, P., and Gasnikov, A. (2023). Algorithms for poolidean-rreularisse optimal transport. In Olenev, N., Evtushenko, Y., Jaćimović, M., Khachay, M., and Malkova, V., editors, Optimization and Applications, Lecture Notes in Computer Science, pages 84-98. Springer Nature Switzerland.
[29] Pele, O. and Werman, M. (2009). Fast and robust earth mover's distances. In 2009 IEEE 12th International Conference on Computer Vision, pages 460-467.
[30] Peyré, G. and Cuturi, M. (2019). Computational optimal transport: With applications to data science. Foundations and Trends@ in Machine Learning, 115-6))355-607
[31] Rockafellar, R. T. (1970). Convex Analysis. Princeton University Press.
[32] Silvia, C., Ray, J., Tom, S., Aldo, P., Heinrich, J., and John, A. (2020). A general approach to fairness with optimal transport. Proceedings of the AAAI Conference on Artificial Intelligence, 34(04))3633-3640.
[33] Sinkhorn, R. (1964). A relationship between arbitrary positive matrices and doubly stochastic matrices. The Annals of Mathematical Statistics, 35(2):876-877 879.
[34] Tang, X., Shavlovsky, M., Rahmanian, H., Tardini, E., Thekumparampil, K. K., Xiao, T., and Ying, L. (2024). Accelerating Sinkhorn algorithm with sparse Newton iterations. In The Twelfth International Conference on Learning Representations.
[35] Torres, L. C., Pereira, L. M., and Amini, M. H. (2021). A survey on optimal transport for machine learning: Theory and applications. arXiv:2 206.:200.000 01..0/96..
[36] Villani, C. (2009). Optimal Transport: Old and New, volume 338 of Grundlehren Der Mathe- matischen Wissenschaften. Springer, Berlin, Heidelberg.
12

---

[37] Xiao, H., Rasul, K., and Vollgraf, R. (2017). Fashion-MNIST: A novel image dataset for benchmarking machine learning algorithms.
[38] Yule, G. U. (1912). On the methods of measuring association between two attributes. Journal of the Royal Statistical Society, 75(6):579-66 642..
13

---

A Additional Experiment Details
A.1
Computing environment
All experiments in this article are conducted on personal computer with an Intel 19- 19-13900K CPU, 32 GB memory, and a Ubuntu 24.10 operating system.
A.2
Motivation
The experiments in Section are designed to reflect two typical uses of OT. In the (Fashion PaashionJMMNIST experiment, images are vectorized as density vectors, and OT used as tool for image morphing and interpolation [27 We follow the experiment settings in the existing literature such as [28] and T8744 be the vectorized and normalized pixel
34 to define a, b, and M: given two images, values of the first image, and similarly define b E IR 784 for the second image. For pixel i in the first image and pixel j in the second image, let (w, h) and (w', h() be their original coordinates in the image, respectively. Then the cost value between pixel i and pixel is
Mij wl- w'1 + 1h- h'1
based on the didisttace, or
Mi (w w w²)² + (h h²)²
base on the squared Euclidean distance. The whole cost matrix M E 78447884 then collects all pairwise cost values Mij-
The ImageNet experiment uses OT statistical distance to measure the difference between two distributions. In this setting, each image is one observation of distribution, and we use OT to compute the (approximate) Wasserstein distance between two classes images. Since each image is mapped to 30-dimensional feature vector, the entries of the cost matrix are simply the iddistances or squared Euclidean distances between the feature vectors two images.
A.3 Impact the regularization parameter
In Figures I and 2 we have shown that the regularization parameter v may affect the performance of optimization algorithms. To further study this effect, we consider a series of equally-spaced r values in the logarithmic domain, and display the run times and number of iterations of each algorithm. For brevity, we focus on the comparison between BCD and SSNS, and the results are given Table T
Table I: Performance comparison between BCD and SSNS under different regularization parameters for the ImageNet experiment in Section The convergence tolerance set to Etol 10-8 Left: cost matrix based on the & disttance Right: cost matrix based on the squared Euclidean distance.
log
Method
Time (s)
Iterations
log (1)
Method
Time (s)
Iterations
BCD
1.628
217
BCD
0.438
59
-2
SSNS
1.523
13
SSNS
2.235
11
BCD
3.765
500
BCD
0.853
114
-2.25
-2.55
SSNS
0.960
20
SSNS
1.066
15
BCD
3.765
500
BCD
3.327
443
-2.5
-2.5
SSNS
0.461
30
SSNS
0.997
BCD
> 3.766
500
BCD
> 3.773
> 500
-2.75
-2.75
SSNS
0.383
57
SSNS
0.529
35
BCD
3.767
500
BCD
3.773
500
-3
-3
SSNS
0.771
120
SSNS
0,458
68
From Table we can find that BCD is very sensitive to the value of 1). When 1) is large, BCD may demonstrate some computational advantages, but when is small, BCD typically fails to meet the error tolerance within 500 iterations. The pattern of SSNS shows some interesting points: when n becomes smaller, the number of iterations also increases, but the overall runtime of SSNS may even
14

---

decrease. This is because smaller values typically result in more sparse Hessian approximations, thus leading faster sparse linear system solving. These findings are consistent with our explanations in Section 5
A.4 Impact of the feature dimension
In the ImageNet experiment in Section 5 we map each image to 30-dimensional feature vector to compute the cost matrix. To study the impact of the feature dimension, we conduct five more experiments with d = 60, 90, 200, 300, 500, under the same setting as in Figure 2 The results are shown in Figure 5 The plots demonstrate similar patterns, implying that the convergence property of SSNS is robust to the feature dimension of input images.
Imspeliet (Trnch Caceette
Imagehitt (Bench Cacorte
ImageWet (ondhis Cassetic playo
Logio Huighal Eror
Lo0) Manorra Emees
BD
X3
Dull
500
0055
Run Brr Sencondel
Tuy 3m (ORCOTcN)
Runtine (ocends)
(a) d = 30
(h) d = 60
(c) = 90
(tench
Imagehit (lench
Imageleer (Thcch Cesettl (laye)
No Enos
100 Mox di Emoss
Deal
LENOS
14153
-10
5055
Run time (erezrdee)
Run 1mre (ONTETOO)
Rintim (ovcend))
(d) = 200
(c) d = 300
(f) = 500
Figure 5: ImageNet experiment with different feature dimensions. The cost matrix is based on the l disttnnce and 1 = 0.001.
A.5 Additional test cases
Besides the numerical experiments in Section 5 here we include two more test cases from each
dataset. The first test case shown in Figure 6 uses fi-distances construct cost matrix, and the second one, demonstrated in Figure 7 is based on squared Euclidean distances. Both test cases use
the regularization parameter 7) = 0.001.
B Proof of Theorems
B.1 Proof of Theorem
In this proof all inequality signs hold elementwisely when applied vectors and matrices.
By Algorithm (T. T3 = A, where each element Ai of A is either zero, or the corresponding have A O. Due to the selection scheme in Algorithm ( each
element Tú of T. Clearly, row sum and each column sum of A is upper bounded by the threshold 8, i.e., Alm < 01, and 1 < d1m, Since the last column of A always zero vector by the design of the algorithm, we also easily get Å1m-1 A1m < 81, and A d < 81mm11 1, where d A'1,-
Note that
A
D1 D2
D H(x) HA =
A
15

---

-37372 102-16090)
102-540281
Imagehet Churchi
0
log Margin Tmor
Erooe
amouch
50
100
150
200
250
200
30
100
350
200
250
300
30
120
250
200
250
300
WAST 0D3 -37372, 546980
RaSS (BI 57372 002-5009811
IrrggsNet Clench vs Churchi
150 Kor to troet
WY
Lng
ATRAAD
0.5
2,0
30
2.5
1.0
2,0
0.5
1.0
2,5
Bun (ERCOnSS)
Run (secontst
jsecende
Figure 6: Additional test cases based on disiistanc and 0.001. Top: Marginal error vs. iteration number for different algorithms three datasets. Bottom: Marginal error time.
101-235, (b2422011)
Bun (D1=239, (02-43991)
ImageNet (lench toni
Log 1511-Mbriinll Erren
Marginal
can
50
109
150
200
250
300
50
100
50
200
250
300
150
150
200
250
300
MNBT 101-220) (D-420011)
PostioretesT 121-239, (02-43991)
ImageNet (Teech va french noml
Lo5 Morgoe trro
0.5
2,5
3.0
0.5
1.
2.0
2.5
1.0
2.5
2.5
Figure 7: Additional test cases based on squared Euclidean distances and r/ = 0.001. Top: Marginal error iteration number for different algorithms three datasets. Bottom: Marginal error time.
so D > O, Dj D11ntm-1=100 11-1mm I 01, and Overall, have D1,+m-1 S - Since D 0, this equivalent
- S 1-110 201mm
|| 25 i= 1,2, ............ # m-
Since D is symmetric, we also have
||D, = |D sn 1= ........+......
Let
Ps Dash
A 1, -m
J
and note that all diagonal elements of D are zero. So by the Gershgorin circle theorem, every
eigenvalue of must be smaller than or equal the maximum of DEs + Ps = |D II, which is upper bounded by -18. Therefore, || maxx (D) < 17 -8. 118
16

---

B.2 Proof of Theorem 2
Consider the matrix
A (rag(r) Ts Tog
E -R(ntm-)x(n-mm-)
diag (c)
and suppose that T = (Tiz). To = (Sig). Clearly, Sij = Tij Aij S Tij. Let
Pk - Σl
Akyl:
k == 1=,....,4 + m - 1,
3yk
and then it is easy to find that
rk + Σ Skj S 2rk
h = 1, ....,, h = n +
Ask + Pk =
Ch-n Ei= Si,k-n S 2ck-n
m m-1
Overall, we have Ars + Pk S U == max l|r [| ). By the Gershgorin circle theorem, every eigenvalue A must smaller than or equal to the maximum value of Auk + Pr. which is upper bounded by U < co.
For some columns of A, say, column k may be the case that A 0, indicating that there is no thresholding on Tx. Note that if we simultaneously permute the columns of I and Tg. then the cigenvalues of A do not change. Therefore, without loss of generality we can assume that the first columns and the last column of A are exactly zero, 0 < d S m (A m = 0 always holds by the design of Algorithm lp. Then can partition T.s as T3 [T(1), T T.m.. where T Raxd and T.m E IR" have strictly positive elements. Accordingly, A can be partitioned as
diag(T1 TÚT T (2)T
(())
T diag
A
diag 1
(T T
Consider the matrix B A sont, where e E R=+m-1 a vector whose first (n + d) elements are ones and the remaining elements are zeros, and is positive scalar. Define
R = Σ
Bayl.
k = 1..... n m
jth
and suppose that S minj Sainssisn.lsss/ Tig.
Then for h = I,.... n, it is easy to find that
m-1
R₂ (n- 1)s
Σ
(Sk)
8) +
Σ
Sky (n d 1)8
Skj
f=1
j=d+1
1=1
For h: == 1,..., d.
R =
1)s =
Σ
T (n d 1)s.
1=1
For A = n d n m. 1,
R4 Σ Skn-n = Σ
Tuk-n
Σ
A A
Tik-n
1=1
i=1
The inequality sign in is strict, since we have assumed that at least one element of A A is nonzero. Moreover,
j=1 Tki
k = 1,
Bu
Tik-n
k + 1, d
Tik-nn
k n + n +
17

---

so overall we can obtain that
[1 j=1 Auj + Tkm - -(n - d)s,
m-
,..... n k + 1, ,.... + d
BEk RR
djs,
Aik-n:
n + m
Now take
minnsisnss/sm Tij
2n
and then we can verify that 8 satisfies the condition k 1 n, we have
S minggisnssssa Moreover, for
Bus - Rr 2 Tkm - 118 Tkm
min 4d
Túj 2
min Téj- [j
For k = = I,.... + d.
BER Rs = (n d)s 2 (n 1)s
m
min T
ij.
2n
For = n d ,,..., n + m 1, since Ai,k-n is either zero Tik--n and at least one element of A k-n is nonzero, we then have
B&s R4
Σ A
min
Tik-n 2 Ty...
18i5n
1.
i=1
Combining all the possibilities, get
m
B&s Rs
2
L:
min Tig
2n
13
for all & and d.
By the Gershgorin circle theorem, every eigenvalue of B must be greater than equal to the minimum value of B& Rx, which is lower bounded by > 0. Since A B + svoT and vu" is positive semi-definite, have m (A) L > 0, implying that A is positive definite.
B.3 Proof of Theorem 3
We first introduce some necessary notation. Define the level set L(xo) = {x f()))() s ((a)),
and let
f* inf( : e L(z)),
81 sup(]g((z)] : L(zo)),
B. = (
E L(xo), > 0}.
Define the function (xx, d) = f(x + d) f(x) -())(((x)))))))))) d. We then present a few technical lemmas, followed by the proof of the main theorem.
Lemma 6. L(xo) is bounded and closed convex set, > oo, BI < oo, and 6₂ < cc.
Proof. Clearly, (x) is continuously differentiable and strictly convex function. Define T., {a: f(x) E c), and then T., is closed convex set.
Lemma of Dvurechensky al. 10 shows that if (a*,8*) an optimal solution 3 then
max 3 - min 37 < Rt
where |M |/o//// log( ( {a, b,))) }). Since have fixed 8 we immediately obtain that || < 0o. Given 8",, a" has a closed form,
at = 17 log a; 7 log
(33 M,y)
18

---

so ||a" is also bounded.
Let x² = (a* "" and f f(=)), and then clearly f > -OC. Take с = c* = {f*) J. and obviously, Le non-empty and bounded. Then Corollary 8.7.1 of [31] shows that Le is bounded for every implying that L(zo)= = ((x)) is bounded. Since g(x) is continuous on L(xp) and L(x)) is compact, it easy to show that $1 < 00.
By Theorem 2 |||||| = mmax (Hs) I 2n maxi||r]lee ||c||<>) where r = Tim and c = 7 In Clearly, every element of T' is bounded if bounded, making ||||| and ||e]|oc also bounded. As result, H&|l is bounded on [.oo), implying that $2 < co.
IT
Lemma For any & > O, there exists a constant 50 such that for n E lc, x e L(xo). and ||p|| s ón satisfying x + Tp E L(ro), we have
Poots |ppl-
(0 (p) <
Proof. Since f is continuously differentiable on [. by the mean value theorem for multivariate functions, there exists some tb E (0, 1) such that for all x ((x))) and x up E [(xo).
f(r mp) f(x) = g + vmp) (mp)
Therefore,
(2x np) g(x + vmp)) (mp) - g(x)) (mp) = [g(x + vmp)) g(x)]² (πp).
By Lemma is continuous [(x0) and L(x²) is compact set, so g is uniformly continuous Lx²)) Therefore, for any = > 0, there exists a constant do > 0 such that whenever l|pl < 50. we have
119( + emp)) - g(x))1 POGE
4cn
As result,
[(x, 19(( + shmp)) ((x) ||pl|| < PaGlE
||p|-
Lemma 8. Let {ox} be generated by Algorithm 2 and then
11gx my (0) mn (EuPR) 2 2 |g min [ISaaaall ll,
[
Proof. Let fx f(xx) and rą [& If we take (p, Mellxl))) and then we can verify that (p. A) meets the following relations:
Ok + HopP P + Ap = 0
A 2 0
Tk ||||| 2 0
A (rk [lpl]) 0.
As (+) is convex function, we have that EkPr is a KKT point and an optimal solution to the constrained optimization problem
1
min mp(p) = fa + of + 70" HoLP
|lll| $ Tk.
(8)
Next, consider the Cauchy point defined in Chapter 4 of Nocedal and Wright 24
Fkk
-h = - Th
|19h6l1
gk.
where
1,
if 9h Ho. 9k S 0
Tk
min
1,
Ige
if of H 9k > 0
FLo VeGr
19

---

Then we can obtain the following inequality by Lemma 4.3 of Nocedal and Wright (24
1
(0) (F.) |lahll min
Th:
[Hs.
Meanwhile, the Cauchy point ph is feasible point for the constrained problem (). and ExPE is an optimal solution for this problem. Therefore, mp(p.)) 2 mnx (CEPR), and we have
1 m& (0) (E(PE)) 2 =19x|11 min
1xll
rk:
|
Let gh and pa be generated by Algorithm If there exists constant e >> and an integer such that [[ 2 for all 2 K, then there must exist sufficiently large constant jr > 0, such that Pk-1 I 4p for all k 2 K.
Lemma 9.
Proof. For a given j. define I, = k K, < p. and 12 : K, uu 2 p... for some In is finite, then it trivially holds that all μız has global upper bound, which leads to the desired conclusion. Therefore, we only consider the case that 1 is infinite for all ş >> O. Since ji can be chosen arbitrarily, {ue]RET, must be unbounded.
Then we will estimate the bound of the following quantity,
(E(Ep))) f(zx + EkPR)
Pk 1 =
mxx(0) - m)((EpP))
According to Lemma one has ||gx|| 2 and |1 S. B1 for k 2 K. Theorem 2 shows that Hag is positive definite, so
O (Ho, + (xlllxl1) 2 MNE,
[HH,, + +x|1gx>11 T) S 1/ 1
By Algorithm
2. we have
EkPr Ex(Ho, + +xl1gg gk.
and hence for all T2,
[IεxPиl S Ea (H + Mell 11 I E
Equation indicates that [[εxPи]|, At E le can be made arbitrarily small with a sufficiently large f. Therefore, can choose some z such that
[k<PPl| S min do, 8 Pacee
Yk Jz, where dn is defined in Lemma 7 Meanwhile, we know
B2 2c S
(10)
(mx(E)P)) fxx <xPx)) ( EuPk)
5 (ExPx) H (EAPH)
1
< ( EaPE))
+
|H (Expk)]l ||EEpEll
2
< ( EuPR))
+ 2
HHs [REePlll ||pk|l-
Lemma indicates that when ||Expæ|| ||
< 6o. we have
(xx, EaPk)l <
PaC4E
l
Therefore, with that satisfies (10). we obtain
(Epp))))) f(xx + EaPx)) 5 POGIE
Cu.
POGE
POCEE
|
||pk||
2
2c,
2
20

---

On the other hand, Lemma indicates that
1
m& (0) my mz (EuPs)
2 19 min [Ksapall
|1
2 2 min [IEapull,
E $32
2 l[Expall 2
CGE
||pxkll
so we have
Ipk 1| mx(Epp))) - f(xx + EaPR)
PuczE
||pell
m²(0) (ExPP))
Po:
which means that Pk: Pn for all T2. Then by the design of Algorithm have μ1++1 S P4 for all & E 12.
Now
show that 16++1 I 4ü for all k 2 K by induction. First, enlarge when necessary to ensure that j > μКk- Then by Algorithm 2 we must have μ I 4ñ. Now suppose that P1+1 S 4 for some 1 K. If 14+11 I1, then clearly M1+2 I 4ü immediately holds. Otherwise, μ4+1 E Iz, so
by the argument above, we have pL42 I pt+1 S 4ü. In both cases, the conclusion holds.
0
Lemma 10. Under the same conditions as in Lemma 9 define
K = fk k 2 K, 2 Po]
Then K must be finite set.
Proof. We use proof by contradiction. Suppose that K. is an infinite set. Since
Pk: [(xx) - s(x + EaPA).
(0 mx(EkPk)
and note that rx+1 -z if Pk 0, we have f(Ix+1) I f(æx EaPR). Therefore, by Lemma 8 holds that
1
f(xR) f(x241) po(m=(() mu(Exp)))
|l min |IEaPS |l.
H
min
[EspPall: B2
2
for all & E K. Combined with Lemma 6 find that f(xxx is monotonically non-increasing with lower bound, so ((xx has a limit, and hence
lim koo
[EPPHll = 0.
On the other hand, = -(Hs, + Ma 11 -'ggs. which means that Øt - -(H5. + 1xll9a 11) Therefore.
E S || [1(Hs Helloul4))pll S | + ps lall ||
1
< [| || 1xl19x ||px||
(B2 + m1331) ||EpPull-
In other words,
0
||Expell 0,
B2 + 33.31
which implies that пı + for & E K,, k + This contradicts with the fact that μ4 < 4 as shown in Lemma 9. So to conclude, X. must a finite set.
Lemma 11. Let { be generated by Algorithml 2 Then either Algorithm terminates in finite iterations, or g (x*) satisfies
lim ininnll inf]]]]g[()]]] 0.
k0000
21

---

Proof. We prove this lemma by contradiction. Suppose that there exist some E > 0 and an integer K such that
|9 > e,
vk 2 K.
Then Lemma 10 implies that the index set
K fk: K, Pk Po)
infinite. This means that there is a sufficiently large integer K' such that Ps < Po for all k 2 K" According to Algorithm 2 must have
Mk++ Apk-
4k 2 K"
which means that μiu - 00. However, this contradicts with the fact that μa S 4 for some p > 0 as shown in Lemma Therefore, we must have
lim inf gi k-400
O.
Then we are ready to prove Theorem 3 Let xt be the unique global optimum of (4). Clearly, f(x) and g(x) are Lipschitz continuous on L (zo), so there exist constants C1, C2 such that
If(*)) f(=")(())))))))))))) S Cill - 2011,
(11)
[19(x*)) - g(x I S C2ll x='
(12)
On the other hand, by taking 5 0, Theorem 2 shows that there is constant C₁ > 0 such that amin (x)) 2 ei for x E L(ro), so there a constant 62 > 0 such that
||o(xx))| ((g(2x)) s(x*)l 2 c2 10 x=
(13)
Also by Taylor's theorem,
f(xx) = [(x)) + [a(x=))"((x> - x²)
(ak xx) H())(((x - x²),
where = some point between uy and xt. Since g(x") = 0. we obtain
((x)) - s(x)) 2 C1 10 - x 11
(14)
2
By the design of Algorithm 2 f(xx) is non-increasing and is lower bounded by f f(x*), so f(xx) must have a limit. Suppose that
lim [f(x)) l"1 =2 k
and then and 13) indicate that for sufficiently large ki,
CT's S CT (T' - f(x*(( S || - x¹11 S cg (lxxx1)
Since in Lemma
11 we have shown that lim inf ||gx|l = 0, cannot be any positive value.
Therefore, = 0. Combining (12) and (14), we have
C1 |
S
C1 |
x 112 S f(xx) - -f(x") + =
203
2
which implies that Dom Iiga|
0 and limp-soo || - x²1 =0 0.
B.4 Proof of Theorem 4
For convenience, let Hy (II)) and B₁ = H8g + ux vull II. Theorem shows that for any 6,
m (Hs) 2 n1 2n
min Tej- 4.
Note that
T= T(a,3) 2 "(a
22

---

and each (a;; + By) is bounded on Lzzo), so each element of T must be bounded below from zero on [(ra). Therefore, there exists constant O < C1 < co such that m (Hs) z CT for all L(xo).
By the design of Algorithm E L(xo) for all k. As a result,
B || min (Bx) }
-1 < Omin (Hs,)l < Ci:
(15)
and hence ||p*|| = w gul I @xll || for all k. Moreover, Theorem shows that |Hs = Tmax (H S C2 for some C2 > 0, so by Lemma we have
mx mx(Ekpk) 2 || min [IEaPall:
1
119x11 }
2
|Ha.
1
|| minn
(alllplli Il=
||px1] C1O2
201
min(cr,
-1lpx|1 = C311px
2C1
On the other hand,
(m<(E)Pk) f(ax + ExPx)]
fxx + EAPA) fxxx) - g (EMPH)
(EAPA)) Ho,
( (ENPK)"
< fxx + CAPA) f(xx) gf (ExPh)
2 (EAPL)" (((PP)))
+ 1 (EEPA)" (H& Ho, (EAPH))
(16)
By Taylor's theorem,
f(Ex + Eupx) f(xx)= - 9% (ExPK) 2 (ExPx)) HA (ExPx)²
= o(llEAPN11?)
For the second term of (6]. let Dr. 11 H5. and then Theorem shows that | S 1 28
By Algorithm
2 5x I
vo VoTT/11 so 1| - Hos || S 1 'vollge 10llgx1' Then
(ELPk)) (Hs- - (EEPE)) -
Po 2n
|| || = of SAPP
2
As result, mmm (ExP))) f(* + [aPH)] = 11 and hence
lpk 1| = mx(ExPk) - f( + (APK)
< <[lEApk|1 Csllpu|12
- 0.
m0) m& (CaPA)
This implies that there is an integer K > 0 such that for all & 2 K. 2 1 - Po. By the design of Algorithm (2, we have p4++1 S and F&+1 2& + EuPR for all l: 2 K.
B.5 Proof of Theorem 5
We first present a classical result derived from the mean value theorem vector-valued functions.
Lemma 12 (Theorem 3.2.12 of (26)). Letg:
g: D C R + Rm be continuously differentiable convex set Do C D and suppose that for constants a 2 and p 2 0, Vg satisfies
[Vg(u) - Vg(o)) = alu - =ll", Vu, E Do-
Then, for any x, E Du.
llq(v) - g (xx) - Vg(z) x) s Ja/(p + 1)] lly x1|
(17)
Then we are ready to prove Theorem Let x be the unique global optimum of (4p. and then g(x*) O. Clearly, g(z) and (xx are Lipschitz continuous on L(xo), so there exist constants L1, La > such that for all v, E [xxo),
(1g(x)) - [(y)] < LLIl --ll--vl y/b || [x)) H(v)l < L2 =v1-
23

---

This implies that g satisfies the conditions in Lemma 12 with a L2 and 1. By substituting t- 24 and y t x². have
||9(xx)l (1(())) g(x")) - L1 Lilla x011;
llo(x)) - g(x)) (xx)((== xx)Il = (1((*)) -s(x") (I x²)ه s 2 2Z x² x=1
Without loss of generality let K' > K, and then for all 2 K', ui S K and 0k+1 2x + Pk. In the proof of Theorem 5 have shown that
|1 S
so for 11 - 1,
| - Bxll S H2 H5. || ue 1 (n (n-'v + k)1l9ul1 S Li + +n)l x²11-
Also note that B
C1 from 15 and then we have
| - x²1 = |2 + Pk - x²|l = || - x² - Bx gkll
= k (Ixx - x²) - -gh)ll
< [ I1| 119h g(x²) - -H&(II - xx")== + H (2 x²) - BB(Lx x²) }
< C (1lgs - g(x") Hv Hx( - xx x²)]l + [|H - B1)( - x² (1)
<
T 1 x. x²112 H B&ll || - x²1l
Le 12 x x=112 Ln (17 2oo 1:) 111 x=" 112)
SQ
2
(2/2 1n2 + 1n(1itup + xi) 1 - x x²11².
G1
24

---

NeurIPS Paper Checklist
1. Claims
Question:
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We have provided detailed theoretical analyses and numerical experiments to support the claims made the abstract.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contri- butions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes]
Justification: We have mentioned in the conclusion part that our proposed method relies on the approximate sparsity of the transport plan.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these ussumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach only tested on few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play important role developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: We have clearly listed the assumptions of each theorem, and have provided the proof in the appendix.
25

---

Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification:
We have provided the code to reproduce the experiments in the article.
Guidelines:
The answer NA means that the paper does include experiments.
If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, may be necessary to either make it possible for others replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, other means that appropriate the research performed.
While NeurIPS does not require releasing code, the conference does require all submissions provide reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., large language model), then there should either be a way to access this model for reproducing the results way reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but should be possible for other researchers to have some path reproducing verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, described in supplemental material? Answer: [Yes]
Justification: All datasets used are publicly available, and we have provided the programming code.
26

---

Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nipss.cc/ppu))
guides/ codeeuubisssionoooee
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this central the contribution (e.g.. for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: Sorniides/ooeeSubii Po for more details.
//nips.cc/publl
The authors should provide instructions on data access and preparation, including how to access the data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only a subset of experiments are reproducible, they should state which ones omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Does the paper specify all the training and test details (e.g., data splits, hyperparame- ters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
Justification:
The detailed settings are reflected in the main article as well as the programming code.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [NA]
Justification: This article is mostly concerned with the optimization convergence properties instead of statistical uncertainties.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main
claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should given (c.g., Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error the mean.
It is OK to report 1-sigma emror bars, but one should state it. The authors should preferably report 2-sigma error bar than state that they have a 96% CI, if the hypothesis Normality of errors is not verified.
27

---

For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are of range (e.g. negative error rates). If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Answer: [Yes]
Justification:
We have provided the details in Appendix A.1.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't tt make it into the paper).
9. Code Of Ethics
Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics hattps://neurre
Question:
Answer: [Yes]
Justification: This research does not involve human subjects, and there is no clear harmful consequence the output of this work.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e-g., if there is special consideration due to laws regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
This article is mostly a theoretical and computational work.
Guidelines:
The answer NA means that there is societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it not needed to point out that generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
28

---

If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g. gated release of models. providing defenses in addition attacks, mecha- nisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Question:
Does the paper describe safeguards that have been put place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification:
This paper does not pose such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We have cited the original papers that provide the datasets.
Guidelines:
existing assets. The authors should cite the original paper that produced the code package dataset.
The answer NA means that the paper does
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from particular source (e.g., website), the copyright and terms service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithee wwwthcode...oooo cwwc /dttasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors encouraged to reach out to the asset's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
We do not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets. Researchers should communicate the details of the dataset/code/moddl part of their sub- missions via structured templates. This includes details about training, license, limitations,
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include anonymized zip file.

---

14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: This article does involve crowdsoureing nor research with human subjects. Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This article does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere the NeurlPS Code of Ethies and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
30


Reviews:
[{"review_id": "VKO5bMjBNw", "reviewer": "Reviewer_vqB6", "summary": "The paper proposes a Newton method to solve the entropic-regularized optimal transform (OT) problem.  The method includes a novel strategy for the sparse approximation of the Hessian to reduce computational complexity compared to the classical Newton method and applying a diagonal shift on the sparse Hessian to avoid singularity. The sparse approximation of the Hessian is a simple procedure based on removing the small off-diagonal block elements;  the authors have theoretically proven that the approximation error is bounded and that the sparse approximate is always positive definite. The authors have demonstrated the efficacy of their method on two OT problems.", "questions": "I find the first experiment setting a bit confusing. Could the authors give the motivation behind this experiment? Why are the pixel values used as probabilities? Please see Weaknesses for some other questions on this experiment.", "limitations": "The authors have mentioned the limitation of their work.", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "1.\tThe paper is very well-written, organized and clear. I congratulate the authors for writing a paper that is this easy to read while they make significant theoretical contributions.\n2.\tThe proposed method is practical with only few tuning parameters.\n3.\tThe theoretical guarantees of the proposed method are impressive and this is important in eliminating the limitations of prior art.\n4.\tThe numerical experiments clearly demonstrate the advantages of the proposed method in the error performance and computational complexity.", "weaknesses": "I have not found any major weaknesses. Minor comments/typos:\n1.\tThe authors have not explained how they chose the pair of images for their first experiment. Is it random? Does the algorithm have similar performance for other images?\n2.\tThe authors have usen $K$ for both the number of small elements and the number of iterations. Please use a different symbol for either one to avoid confusion. \n3.\tLine 67: “proved” should be “proven.”\n4.\tLine 110: The full name of “BFGS” is missing.\n5.\tWhat the authors mean by “safe” only becomes clear at line 169. I would suggest making this clarification earlier in the paper. Maybe, Contribution #2 could be rephrased so that is clear that “safe” refers to avoiding singularity.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the comments. Below are our point-by-point responses for the questions.\n\n### **Weaknesses**\n\n1. To enhance the reproducibility, we do not pick the pairs of images. Instead, the randomly selected image IDs are taken from the prior literature [R1] that studies quadratically regularized OT, and we follow their experiment setting to study entropic-regularized OT. In fact, we have already included other image pairs in Figures 4 and 5 in Appendix A.3, and the IDs of image pairs can be seen from the titles of images. In the rebuttal PDF we have added more test cases in Figure R1.\n\n2-4. Thanks for the suggestions. We have fixed these issues in our local revision.\n\n5. Thanks for the suggestion. We call our proposed algorithm a safe and sparse Newton method, in the sense that the linear systems for computing the search directions are always positive definite. This property addresses the invertibility issues of existing sparsified Newton methods, and is crucial for practical implementation. In our local manuscript, we have added explanations in the introduction, above the contribution section.\n\n### **Questions**\n\nWe are trying to demonstrate two typical uses of OT. In the first experiment, images are vectorized as density vectors, and OT is a used as a tool for image morphing and interpolation [R2]. The experiment setting is already used in previous literature such as [R1] and [R3]. The second experiment uses OT as a statistical distance to measure the difference between two distributions. In this setting, each image is one observation of a distribution, and we use OT to compute the (approximate) Wasserstein distance between two classes of images. In our local revision, we have added a section in the appendix to explain the motivation of these experiments.\n\n[R1] Pasechnyuk, D. A., Persiianov, M., Dvurechensky, P., & Gasnikov, A. (2023). Algorithms for Euclidean-regularised optimal transport.\n\n[R2] Papadakis, N. (2015). Optimal transport for image processing.\n\n[R3] Tang, X., Shavlovsky, M., Rahmanian, H., Tardini, E., Thekumparampil, K. K., Xiao, T., & Ying, L. (2024). Accelerating Sinkhorn algorithm with sparse newton iterations.", "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": "", "extra_scores": {}}, {"review_id": "JnkHP3gD2M", "reviewer": "Reviewer_DpgA", "summary": "This paper proposes a Newton-based algorithm to solve the entropic optimal transport problem on the basis of samples. The approach hinges on a \"sparsification\" scheme for the Hessian (which is explained in Algorithm 1) that retains many favorable properties such as an approximation error due to the sparsification (Theorem 1) and positive definite-ness (Theorem 2). The conclude with numerical experiments on image-data, where they demonstrate that their algorithm significantly outperforms other algorithms for EOT on the examples they considered.", "questions": "See weaknesses.\n\n- Following the above, is it possible to provide a runtime in the style of Sinkhorn given a prescribed tolerance level?\n- Can this approach be adapted for unbalanced entropic OT?", "limitations": "N/A", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The paper is well-written (with very few typos), with the overall story and methodology quite clear.", "weaknesses": "The experiments are mildly underwhelming, and also a bit inconsistent. In the first set of experiments, the authors are running their new (E)OT algorithm between two images, where the pixel frequency denotes the weights for the EOT problem. In the second set of experiments, where they perform (E)OT between two classes of images, there they flatten the images, performing some pre-processing so they lie in d=30, and use $n\\simeq m \\simeq 1000$ samples to constitute the weight vectors ($1/n\\bm{1}$ and $1/m\\bm{1}$). These are relatively small-scale in the number of samples $n$, and I am inclined to believe that it does not scale well when $n=10000$ or more, which the block-coordinate-descent approach (i.e., the vanilla Sinkhorn algorithm) does allow for. I could be wrong, but this is the more realistic/modern use-case for EOT, which is not covered in the experiments.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### **Weaknesses**\n\nThanks for the comments. As we have explained in the global rebuttal and other threads, the experiments are intentionally designed to reflect two typical uses of OT, one for image morphing and interpolation, and the other for computing statistical distances. In our local revision, we have added a section in the appendix to explain the motivation of these experiments.\n\nAs for the scalability, we have added a new set of experiments at the scales of $n=m=1000,5000,10000$, as explained in the global rebuttal. The results show that SSNS is capable of solving entropic-regularized OT on very large problems, and is more efficient than BCD in terms of run time and number of iterations.\n\n### **Questions**\n\n1. We suppose the \"runtime\" here refers to the theoretical computational complexity based on the global convergence speed. First of all, to the best of our knowledge, most convergence rate results for Newton-type methods are local, and there is few on the global convergence rate. However, there is a simple \"ensemble\" method to safeguard the global convergence speed. Note that our Algorithm 2 never increases the objective function value in any iteration, so we can append a Sinkhorn iteration after every $N$ Newton-type steps, where $N$ is a fixed integer. In this way, we obtain at least a Sinkhorn-like convergence speed, but are likely to make larger progresses during the Newton steps. This idea is also mentioned in Section 3.2, page 41 of [R4].\n2. Following the formulation in [R5], entropic-regularized unbalanced OT also has a smooth dual objective function (equation (4) of [R5]). Its difference with the balanced version is that the linear term in equation (4) of our article,\n$$-\\alpha^T a-\\beta^T b=-\\sum_{i=1}^n \\alpha_i a_i-\\sum_{j=1}^m \\beta_j b_j,$$\nis replaced by sums of exponentials,\n$$\\tau \\sum_{i=1}^n e^{-\\alpha_i/\\tau} a_i+\\tau \\sum_{j=1}^m e^{-\\beta_j/\\tau} b_j.$$\nClearly, this means that only the diagonal elements of the Hessian matrix would be different from the setting in our article, and we think that the proposed sparsification scheme and the Newton-type algorithm still apply.\n\n[R4] Nocedal, J., & Wright, S. J. (2006). Numerical optimization.\n\n[R5] Pham, K., Le, K., Ho, N., Pham, T., & Bui, H. (2020). On unbalanced optimal transport: An analysis of Sinkhorn algorithm.", "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": "", "extra_scores": {}}, {"review_id": "kl5odKchUv", "reviewer": "Reviewer_gzRs", "summary": "This paper proposes to solve the entropy-regularized OT problem by proposing a customized Newton method. More specifically, the authors propose a new way to sparsify the Hessian matrix as well as adaptive mechanisms to adjust the hyper-parameters in each iteration. Experimental results are also encouraging.", "questions": "1. What is the computational complexity of Algorithm 1 Sparsification?\n2. What is the computational complexity of each iteration in Algorithm 2?\n3. Can you compare with the SNS method and report results in tables during rebuttal?\n4. The SNS paper uses the metric Log of Optimality Gap while you use the metric Log of Marginal Errors? What is the difference? Can you report your results in terms of Log of Optimality Gap as well during rebuttal?\n5. Can you also replicate the Table 3 and Table 4 in the SNS paper? Basically conduct perturbation study on the entropy regularization parameters.\n6. For the imageNet dataset, how well does your proposed method perform if you change the final feature dimension from 30 to 60 and 90, respectively?", "limitations": "NA", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The new sparsification method is novel and different from [31]. I like the fact that you can prove that this method positive definite matrices.\n\n2. The self-adaptive mechanism to adjust Newton method hyper-parameters look interesting. I am not 100% positive (considering that there are already line search methods to guarantee convergence of Newton methods), but this part seems novel to me. The authors only cite very old papers [16, 17, 21] as motivations.", "weaknesses": "1. The major part is the scalability of the proposed method. Newton-variants methods are second-order methods. The largest scale the authors have tried is on the transformed features of ImageNet. The feature dimension is only 30. If we use the Sinkhorn algorithm, I believe we can handle much large scales.\n\n2. The experimental results are only sample results based on very few samples from the dataset. I am concerned that the results could be cherry-picked. It would be much more comprehensive if the authors can present training efficiency improvement on the overall dataset.\n\n3. There is no experimental comparison with SNS [31] even though this work takes a step further on top of SNS. Looking at the SNS paper, I believe it is very easy to implement.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for the comments. Below are our point-by-point responses.\n\n### **Weaknesses**\n\n1. The dimension of the features does not impact the scale the problem. The extracted features are only used to compute the cost matrix, which has a size of $n\\times m$. We do not make the feature dimension $d$ too large mainly because it is commonly accepted that the Euclidean distance suffers from the curse of dimensionality. In the global rebuttal, we have added experiments to study the impact of feature dimension and the scalability of algorithms on very large problems.\n2. As we have mentioned in the global rebuttal and other threads, to enhance the reproducibility, we do not pick the pairs of images. Instead, the randomly selected image IDs are taken from the prior literature [R1] that studies quadratically regularized OT, and we follow their experiment setting to study entropic-regularized OT. We have included different examples in Figures 4 and 5 in Appendix A.3. Also, to investigate the scalability of algorithms on very large problems, we have added experiments to test the performance for $n=m=1000,5000,10000$.\n3. We agree that the overall framework of SNS [R3] is easy to understand, but the major challenges for directly comparing with SNS are: (a) there is no publicly available code of SNS to reproduce the experiments; (b) in [R3], the algorithm has many hyperparameters, and we have no clear scheme to set those hyperparameters. For example, how many Sinkhorn iterations to run before switching to the Newton step, how large the sparsification parameter should be, what the fallback strategy should be when the linear system is not invertible, etc. In fact, one of the major motivations of this article is to make the SNS framework more adaptive and practical, and we position our SSNS algorithm as a concrete and practical variant of the original SNS method. \n\n### **Questions**\n\n1. A trivial upper bound of the computational complexity of Algorithm 1 is $O(n^2\\log(n))$, assuming $n=m$. This is done by sorting the values of each column or row of $T$. However, if we have an estimated upper bound $\\varrho$ of the density of the sparsified $T$, meaning that each column or row of $T$ has at most $\\varrho n$ nonzero elements after sparsification, then the cost can be reduced to $O(n^2+\\varrho n^2\\log(\\varrho n))$. This is done by selecting the largest $\\varrho n$ elements in each column or row of $T$, and only sort values within these $\\varrho n$ elements. In a typical setting, $\\varrho$ is very small due to the approximate sparsity of $T$.\n2. Each iteration of Algorithm 2 has a computational cost of $O(n^2)$ to compute the gradient and objective function value, plus the cost in Algorithm 1, plus the cost in solving the sparse linear system. When $T$ is very sparse, the number of nonzero elements of $T$ may be as small as $O(n)$, so the overall computation is dominated by the first and second parts, *i.e.*, $O(n^2)$.\n3. As we have explained above, there are certain practical difficulties in directly comparing with SNS. But we have tried our best to provide more experiment results similar to those in the SNS paper, as introduced in the points below.\n4. A typical definition of the optimality gap is the value $f(x_k)-f(x^*)$, where $f(\\cdot)$ is the objective function, and $x^*$ is the optimal point. The marginal error in entropic-regularized OT coincides with the gradient norm $\\Vert g(x_k) \\Vert$. We do not use the optimality gap because the ground-truth optimal point $x^*$ is in general unknown, and needs to be approximated using existing solvers. But this process also introduces rounding errors, and may be very sensitive when $x_k$ is close to convergence. Instead, $\\Vert g(x_k) \\Vert$ is easy to compute, and has an absolute lower bound of zero. Overall, we think that the marginal error is a more objective and easy-to-compute criterion to evaluate convergence.\n5. Yes, we have added such experiments in the global rebuttal.\n6. Per the suggestion, we have added such experiments in the global rebuttal (Figure R2).\n\n[R1] Pasechnyuk, D. A., Persiianov, M., Dvurechensky, P., & Gasnikov, A. (2023). Algorithms for Euclidean-regularised optimal transport.\n\n[R3] Tang, X., Shavlovsky, M., Rahmanian, H., Tardini, E., Thekumparampil, K. K., Xiao, T., & Ying, L. (2024). Accelerating Sinkhorn algorithm with sparse newton iterations.", "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": "", "extra_scores": {}}, {"review_id": "5Ok3PRThz3", "reviewer": "Reviewer_bgcB", "summary": "This paper proposed a new Newton-type algorithm for entropic-regularized optimal transport (OT) which utlizes sparsification and safeguard techniques and achieves global convergence and local quadratic convergence for the entropic-OT problem. Numerical experiments are provided to verify the effectiveness of the proposed method.", "questions": "In general, I think this is a technically solid paper. The authors identify the issue of applying second order methods to entropic-OT and use sparsification and safeguard to make it applicable. I didn't check the detail of the proof, but believe it to be correct since these all follow the standard techniques. I'd recommend this work to be accepted.\n\nI have the following comments and questions:\n\n1. Regarding Theorem 2 which studies the positive-definiteness of the sparsified matrix $H_{\\delta}$: what is the largest and smallest eignvalues if we don't conduct Algorithm 1, meaning we use the original Hessian in (5) directly? The author should add some discussion on this to consolidate their point of \"safely used to compute the Newton search directions\".\n\n2. When the authors talk about \"safe\", does it mean Theorem 2 or the trust region type update on line 7-13 in Algorithm 2, or both?\n\n3. In the experiments for MNIST and ImageNet, the authors should provide the dimension information corresponds to the theory part, such as what $n$ and $m$ are.\n\n4. In all images, Newton method seems to be very slow. Is this all due to the costly process of solving the dense linear system?\n\n5. What is $\\tilde{T}$ in (5)? I feel that it's not clearly defined. Is it corresponds to $\\tilde{\\beta}$ which is $T$ having one row/column removed?", "limitations": "The authors are clear about their limitations.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well-rounded and well-presented. The paper addresses the computational issue caused by a dense Hessian by sparsification, and addresses the singularity issue by a safeguard mechanism. Both the local and global convergences are provided. Numerical experiments are provided to show the effectiveness of the proposed method.", "weaknesses": "(Please answer the Questions section directly) The safeguard mechanism in this paper is actually a standard trust region technique; The numerical experiment settings could be improved.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### **Weaknesses**\n\nThanks for the comments. First of all, we shall clarify that SSNS is directly inspired by the regularized Newton method [R6] and the Levenberg–Marquardt algorithm [R7,R8], which have close relations to the trust-region methods but with some subtle differences. Specifically, standard trust-region methods first set a trust-region radius $\\Delta_k$ in each iteration, and then (approximately) solve a subproblem to determine the search direction. In SSNS, we update the shift parameter $\\lambda_k$ instead, and the search direction has a closed form. The formula of search directions in trust-region methods has a similar form, but its shift $\\lambda_k$ is implicitly determined by $\\Delta_k$, and typically requires some root-finding techniques. Also, SSNS allows for a step size selection procedure as in Algorithm 3, thus enhancing the flexibility.\n\nAlso per the suggestions of reviewers, we have improved the numerical experiments as explained in the global rebuttal.\n\n### **Questions**\n\n1. Recall that in Theorem 2 we allow $\\delta=0$, which corresponds to the case that we use the genuine Hessian matrix. In this setting, Algorithm 2 reduces to a conventional regularized Newton method, but its per-iteration cost would be huge. Theorem 2 is an improvement to existing sparsifed Newton methods such as [R3], since previous works do not guarantee that the sparsified Hessian is invertible (but note that the true Hessian matrix is always invertible by Theorem 2).\n2. In our original design of the algorithm, \"safe\" mainly refers to the positive definiteness of the sparsified Hessian matrix. But from the experiments, we also find that the shift parameter (which has a close relation to trust-region techniques) is very useful, as the classical Newton method fails in some examples due to numerical instability.\n3. Thanks for the suggestion. We have added such information in our local revision. Basically it is $n=m=784$ for (Fashion-)MNIST and $n\\approx m\\approx 1000$ for ImageNet data.\n4. Yes. This can be inferred from the fact that the number of iterations of the Newton method is very small, but its runtime is large. The huge per-iteration cost is the consequence of a dense linear system.\n5. Yes. We have defined this notation in the beginning of Section 2, meaning removing the last column of $A$.\n\n[R3] Tang, X., Shavlovsky, M., Rahmanian, H., Tardini, E., Thekumparampil, K. K., Xiao, T., & Ying, L. (2024). Accelerating Sinkhorn algorithm with sparse Newton iterations.\n\n[R6] Li, D. H., Fukushima, M., Qi, L., & Yamashita, N. (2004). Regularized Newton methods for convex minimization problems with singular solutions.\n\n[R7] Levenberg, K. (1944). A method for the solution of certain non-linear problems in least squares.\n\n[R8] Marquardt, D. W. (1963). An algorithm for least-squares estimation of nonlinear parameters.", "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": "", "extra_scores": {}}]

Author rebuttal:
# To All Reviewers

Thank you all reviewers for the encouraging and insightful comments. We appreciate the time and effort the reviewers have dedicated to providing valuable feedback on our manuscript. In this round, we have made every effort to address the comments of the reviewers. **The point-by-point responses are provided in the reply**.

Additionally, we would like to take this opportunity to make some global clarifications as well as improvements we have made during this period.

### **Explaining what "safe" means**

We call our proposed algorithm a safe and sparse Newton method, in the sense that the linear systems for computing the search directions are always positive definite. This property addresses the invertibility issues of existing sparsified Newton methods, and is crucial for practical implementation.

### **Motivations on the experiment setting**

We are trying to demonstrate two typical uses of OT when designing the numerical experiments. In the first experiment, images are vectorized as density vectors, and OT is a used as a tool for image morphing and interpolation [R2]. The experiment setting is already used in previous literature such as [R1] and [R3]. The second experiment uses OT as a statistical distance to measure the difference between two distributions. In this setting, each image is one observation of a distribution, and we use OT to compute the (approximate) Wasserstein distance between two classes of images.

### **Improvements on experiments**

Per the suggestions of reviewers, we have improved the numerical experiments on the following aspects:

1. We make it clear that to enhance **the reproducibility**, we do not pick the pairs of images in the numerical experiments. Instead, the randomly selected image IDs are taken from the prior literature [R1] that studies quadratically regularized OT, and we follow their experiment setting to study entropic-regularized OT. In the article, we have included other image pairs in Figures 4 and 5, and in the rebuttal PDF we have added more test cases in Figure R1.
2. We have studied **the impact of feature dimension** $d$ in the ImageNet experiment (Figure R2, rebuttal PDF). The plots imply that the convergence property of SSNS is robust to the feature dimension of input images.
3. We have studied **the impact of regularization parameters** on the performance of optimization algorithms. The table below shows the performance comparison between BCD and SSNS under different regularization parameters for the ImageNet experiment in Section 5. The convergence tolerance is set to $\varepsilon_{tol}=10^{-8}$, and the cost matrix is based on the $\ell_{1}$-distance. 

|$\log_{10}(\eta)$|BCD Time (s)|BCD Iterations|SSNS Time (s)|SSNS Iterations|
|-|-|-|-|-|
|-2   | 1.628| 217|1.523|13 |
|-2.25|>3.765|>500|0.960|20 |
|-2.5 |>3.765|>500|0.461|30 |
|-2.75|>3.766|>500|0.383|57 |
|-3   |>3.767|>500|0.771|120|

The table below shows the case using cost matrices based on the squared Euclidean distances.

|$\log_{10}(\eta)$|BCD Time (s)|BCD Iterations|SSNS Time (s)|SSNS Iterations|
|-|-|-|-|-|
|-2   |0.438 |59  |2.235|11|
|-2.25|0.853 |114 |1.066|15|
|-2.5 |3.327 |443 |0.997|23|
|-2.75|>3.773|>500|0.529|35|
|-3   |>3.773|>500|0.458|68|

The results show that BCD is very sensitive to the value of $\eta$. When $\eta$ is large, BCD may demonstrate some computational advantages, but when $\eta$ is small, BCD typically fails to meet the error tolerance within 500 iterations. The pattern of SSNS shows some interesting points: when $\eta$ becomes smaller, the number of iterations also increases, but the overall runtime of SSNS may even decrease. This is because smaller $\eta$ values typically result in more sparse Hessian approximations, thus leading to faster sparse linear system solving. These findings are consistent with our explanations in Section 5.

4. We have investigated **the scalability** of SSNS on very large OT problems. We consider a synthetic OT problem that can generate data with arbitrary dimensions. The basic setting is to approximate the OT between two continuous distributions: the source is an exponential distribution with mean one, and the target is a normal mixture distribution $0.2\cdot N(1,0.2)+0.8\cdot N(3,0.5)$. We discretize the problem in the following way: let $x_i=5(i-1)/(n-1)$, $i=1,\ldots,n$, and $y_j=5(j-1)/(m-1)$, $j=1,\ldots,m$, which are equally-spaced points on [0, 5]. Define the cost matrix as $M_{ij}=(x_i-y_j)^{2}$. Let $f_1$ and $f_2$ be the density functions of the source and target distributions, respectively. Then we set $\tilde{a}_i=f_1(x_i)$, $\tilde{b}_j=f_2(y_j)$, $a\_i=\tilde{a}\_i/\left(\sum\_{k=1}^n\tilde{a}\_k\right)$, and $b\_j=\tilde{b}\_j/\left(\sum\_{k=1}^m\tilde{b}\_k\right)$. Similar to the experiment setting in Section 5, we normalize the cost matrix and set $\eta=0.001$. We then solve the entropic-regularized OT problem using BCD and SSNS at the scales of $n=m=1000,5000$, and 10000. The results are visualized in Figure R3 in the rebuttal PDF, whose pattern is clear: BCD demonstrates a linear-like convergence rate, and SSNS has a fast convergence speed consistent with the theoretical quadratic rate. Thanks to the Hessian sparsification, SSNS does not suffer from a high per-iteration cost, so overall it provides an efficient solver for entropic-regularized OT even on very large problems.

[R1] Pasechnyuk, D. A., Persiianov, M., Dvurechensky, P., & Gasnikov, A. (2023). Algorithms for Euclidean-regularised optimal transport.

[R2] Papadakis, N. (2015). Optimal transport for image processing.

[R3] Tang, X., Shavlovsky, M., Rahmanian, H., Tardini, E., Thekumparampil, K. K., Xiao, T., & Ying, L. (2024). Accelerating Sinkhorn algorithm with sparse newton iterations.

Meta-review: From the reviews, the reviewers reach a consensus that the paper contains interesting result and novel contribution to the community. The Author-Reviewer discussions were also substantial, the authors should incorporate the important ones to the final version of the paper or supplementary to futher strengthen the paper. During the Reviewer-AC discussions, there are two points raised regarding the numerical experiments, if time allows, it would be great if the authors conduct the experiments
 - From Reviewer DpgA: in the comparison against BCD, what the result would be if the tolerance is relatively larger (still acceptable)? As for this setting, BCD would be much faster.
 - From Reviewer gzRs: For showing the scalability of the method, it is applied to solve OT on large-scale synthetic datasets. However, on the real-world datasets such as imageNet (Figure R2 in the rebuttal pdf), the authors have only conducted on dimension of 30, 60, and 90. For revision, it would be great if the authors could solve OT on the imageNet for dimension of 500 and 1000, respectively. If they can show that the same computational advantage can hold on real-world datasets.
Decision: Accept (poster)