Safe and Sparse Newton Method for Entropic-Regularized Optimal Transport

Authors: zihao  tang
yixuan  qiu
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 17:20:10

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.

Main paper: 
 introduction: 
 In recent years, optimal transport (OT, [35] ) has received massive attentions from the deep learning community, and has become a fundamental modeling tool in modern statistical machine learning [34, 22] . One major challenge of applying OT to real-life problems, the computation of large-scale OT, has also gained many progresses in the direction of approximate OT methods, especially the entropic-regularized OT [6] . i=1 and {y j } m j=1 are two sets of data points, and δ x is the Dirac measure at position x. Without loss of generality we assume that n ≥ m, as their roles can be exchanged. Let M = (M ij ) ∈ R n×m be a cost matrix, whose entry M ij typically represents the distance between data points x i and y j . Also define Π(a, b) = {T ∈ R n×m : T 1 m = a, T T 1 n = b, T ≥ 0}, where the inequality sign applies elementwisely. Then OT between the two measures µ and ν can be characterized by the following linear programming problem,EQUATIONwhere ⟨A, B⟩ = tr(A T B). Assuming n = m, the computational complexity of solving (1) using standard linear programming solvers is typically at the order of O(n 3 log(n)) [28] , which can be difficult even for moderate n and m. One approach to approximating the solution is to add an entropic penalty term to the objective function, leading to the entropic-regularized OT problem [6] :EQUATIONwhereh(T ) = n i=1 m j=1 T ij (1 -log T ij )is the entropy term, and η > 0 is a regularization parameter. The objective function in (2) is η-strongly convex on Π(a, b), so (2) has a unique global solution, denoted as T * . The matrix T * is often referred to as the Sinkhorn transport plan. Problem (2) can be solved by the celebrated Sinkhorn algorithm [37, 32] , which is based on efficient matrix-vector multiplication operations. Due to its computational advantage, the Sinkhorn algorithm has become the de facto solution to large-scale OT problems for a 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-order methods, such as the Newton method, have attracted a growing number of researchers to rethink the solution to (2) [3, 33] . It has been shown that (2) has a dual problem (Proposition 4.4 of [29] ):EQUATIONLet α * = (α * 1 , . . . , α * n ) T and β * = (β * 1 , . . . , β * m ) T be one optimal solution to (3), and then the Sinkhorn transport plan T * can be recovered as T * = τ (α * , β * ), where for two vectors α = (α 1 , . . . , α n ) T and β = (β 1 , . . . , β m ) T , τ (α, β) is a matrix with entries [τ (α, β)] ij = exp η -1 (α i + β j -M ij ) .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-order 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 it typically requires much fewer iterations than the Sinkhorn algorithm. However, each iteration in the Newton method involves solving a dense linear system, resulting in a per-iteration cost of O(n 3 ), which is nearly of the same order as that of an unregularized OT problem. To this end, [33] proposes the Sinkhorn-Newton-Sparse (SNS) algorithm, which approximates the Hessian matrix of (3) by a sparse one, so that each iteration roughly has a computational cost of O(n 2 ), the same order as the Sinkhorn algorithm.The SNS algorithm demonstrates promising results for solving the entropic-regularized OT, yet there remain a number of open questions. First, the SNS algorithm relies on initial Sinkhorn iterations to achieve a 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 is no strong guarantee on the invertibility of the sparsified Hessian matrix. Third, [33] conjectures that SNS has a super-linear local convergence rate, but it is not yet proven by theoretical analysis.In this article, we make major progresses on the second-order method for entropic-regularized OT by resolving the three challenges above. First, we carefully design a 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 a global convergence guarantee, so no initial Sinkhorn iterations are needed. Most importantly, we provide solid theoretical analysis of the proposed method, and show that it achieves a 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 a novel Hessian sparsification scheme that promises a strict control of the approximation 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 to a 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 hyperparameter tuning. An efficient implementation is included in the RegOT Python package2 . 4. We provide rigorous global and local convergence analysis of the proposed algorithm, which is lacking in prior literature. 
background and related work: 
 Notation Throughout this article, we use f (x), g(x), and H(x) to represent the objective function, gradient, and Hessian matrix at point x, respectively. For a matrix A, let A i• be the vector of the i-th row of A, and A •j be the vector of the j-th column of A. For a vector v = (v 1 , . . . , v n ) T , let ṽ denote the vector (v 1 , . . . , v n-1 ) T , and for a matrixA = (A •1 , . . . , A •n ), let Ã represent the matrix (A •1 , . . . , A •,n-1). The symbol ∥ • ∥ stands for the Euclidean norm for a vector argument, and represents the operator norm for matrices. The ℓ 1 -norm and infinity norm are denote by ∥ • ∥ 1 and ∥ • ∥ ∞ , respectively, applying to both vectors and matrices.The main objective As has been explained in Section 1, the key to the computation of entropicregularized OT is to solve its smooth dual problem (3) , which is equivalent to the convex optimization problem min α,β -L(α, β). However, it is worth noting that the variables (α, β) have one redundant degree of freedom: if (α * , β * ) is one solution to (3), then so is (α * + c1 n , β * -c1 m ) for any c. Therefore, to ensure identifiability, we globally set β m = 0 without loss of generality. In what follows, we use the vector x = (α T , βT ) T ∈ R n+m-1 to collect all free dual variables. Clearly, we always have β = ( βT , β m ) T = ( βT , 0) T , so β and β will be used interchangeably when we consider functions of β. Then the main target of this article is to efficiently solve the optimization problemEQUATIONIt is known that f (x) is a strictly convex function, so if (4) has a solution, then it is unique. It 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 [20] ):EQUATIONH(x) = - ∇ 2 α L(α, β) ∇ β (∇ α L(α, β)) ∇ β (∇ α L(α, β)) T ∇ 2 β L(α, β) = η -1 diag(T 1 m ) T T T diag( T T 1 n ) .Computational OT There are a huge number of methods developed to solve (2) and ( 4). The Sinkhorn algorithm [37, 32, 6] can be interpreted as applying the block coordinate descent (BCD), also known as the alternating minimization (AM) 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 [1] . Accelerated first-order methods, such as the adaptive primal-dual accelerated gradient descent (APDAGD, [10] ) and adaptive primal-dual accelerated mirror descent (APDAMD, [18] ), 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-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method [19] , have also been used [7] . Second-order methods are less common in the literature to solve (4), mainly due to the high cost of solving the linear systems, and [33] proposes a practical sparsified Newton method to approximate the true Hessian by a sparse matrix.The Newton method is a classical second-order optimization algorithm that has been extensively studied. Assuming f is twice continuously differentiable, the Newton method generates a sequence of iterates {x k } using the updating formulax k+1 = x k -[H(x k )] -1 g(x k ).However, computing the linear system [H(x k )] -1 g(x k ) with a dense Hessian matrix generally costs O(n 3 ) of computation, so the sparsification of H(x k ) with a 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 k )] -1 g(x k ). 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 to characterize the transformation of statistical distributions, OT has wide applications in deep generative models [2, 11, 14] , domain adaptation and transfer learning [5, 4] , and fairness in machine learning [12, 24, 31] , among many others. Readers are referred to review articles such as [34, 22] for a summary of machine learning tasks and methods that utilize OT. 
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: a new scheme to sparsify the Hessian matrix, and a Newton-type algorithm to update iterates.Ideally, the sparsified Hessian matrix, denoted by H δ (x), should meet two criteria: first, it should be close to the true H(x) with a tunable approximation error; second, it needs to preserve the positive definiteness of the original Hessian matrix. The first point is to ensure that H δ (x) does not lose too much information of the true H(x). The second point is especially important in implementing the Newton method, since a linear system [H δ (x)]d = -g(x) needs to be solved to compute the search direction d, and one needs to make sure that H δ (x) 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 H(x), 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(v, δ), which takes a vector v = (v 1 , . . . , v n ) T and a scalar δ > 0 as inputs, and outputs a mask vector ϕ = (ϕ 1 , . . . , ϕ n ) T in the following way.Suppose that v π(1) ≤ v π(2) ≤ • • • ≤ v π(n)are the sorted values of v, where π(s) is the original index of the s-th smallest element in v. We also define π -1 (•) to be the inverse of π(•), i.e., π -1 (i) is the integer such that π(π -1 (i)) = i. Let S be the largest integer such that S s=1 v π(s) ≤ δ. Then we set ϕ i = 1 if π -1 (i) ≤ S, and ϕ i = 0 otherwise. For example, if v = (2, 1, 3, 5, 2) T and δ = 6, then we have ϕ = (1, 1, 0, 0, 1) T , since1 + 2 + 2 ≤ δ but 1 + 2 + 2 + 3 > δ.Next, define an operator apply_mask(v, ϕ), which inputs a vector v = (v 1 , . . . , v n ) T and a mask ϕ = (ϕ 1 , . . . , ϕ n ) T , and outputs a vector u = (u 1 , . . . , u n ) T with u i = ϕ i •v i . For the example above, we have apply_mask(v, ϕ) = (2, 1, 0, 0, 2) T . Then given the current dual variables x = (α T , βT ) T , Algorithm 1 outputs a sparse matrix H δ as an approximation to the true Hessian matrix H(x).It is worth noting that in the construction of H δ , the diagonal elements diag(T 1 m ) and diag( T T 1 n ) use the original T matrix, whereas the off-diagonal elements are the sparsified T δ 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 1 (Approximation error of sparsification). Let H δ be output by Algorithm 1 with a given vector x, and define D = H(x) -H δ . Then for any δ ≥ 0, we have D ≥ 0, where the inequality sign holds elementwisely, and for each i, j = 1, 2, . . . , n + m -1,∥D i• ∥ 1 ≤ η -1 δ, ∥D •j ∥ 1 ≤ η -1 δ.Moreover, ∥D∥ ≤ η -1 δ.Algorithm 1 Sparsifying the Hessian matrix.Input: Dual variable vector x = (α T , βT ) T , threshold parameter δ ≥ 0 Output: Sparsified Hessian matrix H δ 1: Initialize a zero matrix ∆ ∈ R n×m and compute T = τ (α, β) 2: for j = 1, 2, . . . , m -1 do3: ϕ ← select_small(T •j , δ), ∆ •j ← apply_mask(T •j , ϕ) 4: for i = 1, 2, . . . , n do 5: ϕ ← select_small(∆ i• , δ), ∆ i• ← apply_mask(∆ i• , ϕ) 6: T δ ← T -∆ 7: H δ ← η -1 diag(T 1 m ) Tδ T T δ diag( T T 1 n )Theorem 1 shows that all the row-wise and column-wise ℓ 1 -norms of the residual matrix D is bounded by the pre-specified threshold δ, up to a constant multiplier η -1 . More importantly, the overall operator norm of D also has the same upper bound, and this property plays a key role in analyzing the convergence behavior of the SSNS algorithm later in Section 3.2.Next, Theorem 2 validates our second criterion on H δ : the sparsified Hessian matrix by Algorithm 1 is guaranteed to be positive definite, regardless of the threshold parameter δ. Therefore, it can be safely used to compute the Newton search directions.Theorem 2 (Positive definiteness). Suppose that n ≥ m. For any α ∈ R n , β ∈ R m-1 , and M ∈ R n×m , let T = τ (α, β), r = T 1 m , and c = T T 1 n . Then for any δ ≥ 0, the matrix H δ output by Algorithm 1 is always positive definite. Specifically, we haveσ max (H δ ) ≤ 2η -1 • max{∥r∥ ∞ , ∥c∥ ∞ }, σ min (H δ ) ≥ η -1 • n -m + 1 2n • min i,j T ij ,where σ max (•) and σ min (•) represent the largest and smallest eigenvalues of a symmetric matrix, respectively.The two features presented by Theorem 1 and Theorem 2 are crucial in designing the optimization algorithm and analyzing its theoretical properties, which we introduce in details in Section 3.2.In Theorem 2, we have shown that the sparsified Hessian matrix H δ is always positive definite, thus invertible. However, in actual computation, H δ may be nearly singular, thus causing numerical instability. To further stabilize the optimization procedure, we propose a safe and sparse Newton method inspired by the regularized Newton method [17] and the Levenberg-Marquardt algorithm [16, 21] .One of the key ingredients of SSNS is to introduce a shift parameter λ in solving the Newton linear system, leading to a search direction of the form p = -(H δ + λI) -1 g(x) in each iteration. The shift parameter λ has multiple functions. First, it stabilizes the vanilla Newton linear system H -1 δ g(x), as the matrix H δ + λI has a smaller condition number than H δ , which potentially accelerates iterative linear solves such as the CG method. Second, if we allow λ to vary along iterations, λ has the effect of adjusting the search direction when -H -1 δ g(x) is not a good candidate. Third, as we will later show in the convergence analysis, an adaptive λ plays an important role in both the global and local convergence of the algorithm.On the other side, the sparsification threshold δ is also essential to the design of the algorithm, as we need to control the approximation error of H δ during the iterations. Intuitively, when the current iterate x k is close to the optimum, H δ should be also close to the true Hessian matrix H(x k ) to achieve the fast convergence rate, whereas when x k is far away, it may be beneficial to use a large δ, since H δ is then more sparse and leads to a faster computation of the linear system -(H δ + λI) -1 g(x).Overall, in SSNS we allow δ and λ to vary in each iteration, and use the notation δ k and λ k to reflect this adaptivity. Intuitively, both δ k and λ k cause approximation errors to the true Hessian matrix in determining the Newton search direction, so they should tend to zero when x k is sufficiently close to the optimum. On the other hand, we do not want δ k and λ k to decay too fast, since otherwise they make H δ too dense and the linear system may be ill-conditioned. To this end, we set their values according to the current gradient norm ∥g(x k )∥, which can be viewed as an indicator of the distance between x k and the optimum x * . Specifically, we set λ k = µ k ∥g(x k )∥ and δ k = ν 0 ∥g k ∥ γ , where µ k is a parameter that dynamically changes in each iteration, and ν 0 > 0 and γ ≥ 1 are two constants. These settings are justified by the convergence analysis later in Theorems 3 and 5.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 a quantity ρ k ,ρ k = f (x k ) -f (x k + ξ k p k ) m k (0) -m k (ξ k p k ) ,where ξ k > 0 is a step size parameter, and the m k (•) function is a local quadratic approximation to the objective function:EQUATIONThe quantity ρ k has two functions: first, it adjusts the µ k parameter, and second, it also determines whether one should accept the proposed new iterate x k + ξ k p k . A negative value of ρ k means that the objective function actually increases at the new iterate, so one should reject the proposed point, and increase λ k in the next iteration for a potentially better search direction. The overall SSNS algorithm is then given in Algorithm 2.Algorithm 2 Safe and sparse Newton method for Sinkhorn-type optimal transport.Input: Initial point x 0 , parameters {µ 0 , ν 0 , c l , c u , κ} > 0, γ ≥ 1, ρ 0 ∈ (0, 1 2 ), ε tol > 0 Default values: µ 0 = 1, ν 0 = 0.01, c l = 0.1, c u = 1, κ = 0.001, γ = 1, ρ 0 = 1 4 Output: x k 1: for k = 0, 1, 2, . . . do 2: Compute g k = g(x k ), δ k = ν 0 ∥g k ∥ γ 3: if ∥g k ∥ < ε tol then 4: return x k 5:Compute H δ k according to Algorithm 1 with x ← x k 6:Computep k = -(H δ k + µ k ∥g k ∥I) -1 g k 7:Select anyEQUATIONUpdateµ k+1 =    4µ k , if ρ k < ρ 0 max{µ k /2, κ}, if ρ k ≥ 1 -ρ 0 µ k , otherwise 10:if ρ k > 0 then 11:x k+1 = x k + ξ k p k 12:else 13:x k+1 = x kIt is worth noting that in Algorithm 2, the step size ξ k can be taken an arbitrary value from a fixed interval [c l , c u ]. This is quite different from typical line search methods, which require the new iterate x k + ξ k p k to satisfy certain conditions, such as the sufficient decrease condition and the curvature condition [23] . In other words, it is one of the advantages of SSNS that the cost of step size selection is predictable: we can consider a fixed number of candidates, and select one of them according to some criterion. At one extreme, it is completely acceptable to always set ξ k = 1. However, in practice, we advocate Algorithm 3 to 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 a decrease in function value is 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 (ξ [0] , ξ [1] , ξ [2] , ξ [3] ) = (1, 0.5, 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 1 = ξ [0] > ξ [1] > • • • ξ [N ] > 0, current x k , p k , objective function f (•) 1: Initialize ξ * = 1, f * = +∞ 2: for i = 0, 1, . . . , N do 3: Compute x trial = x k + ξ [i] p k , f trial = f (x trial ) 4: if f trial < f * then 5: ξ * = ξ [i] , f * = f trial 6: if f * < f (x k ) then 7: return ξ * , f * 8: return ξ * , f * 
theoretical guarantees on convergence: 
 In this section, we present our major theoretical contributions to the sparsified Newton method for entropic-regularized OT. We first show that starting from an arbitrary initial point x 0 , the iterates generated by Algorithm 2 eventually converge to the unique global optimum x * . Therefore, SSNS does not require a warm initialization using the Sinkhorn algorithm as in [33] . Theorem 3 (Global convergence guarantee). Let {x k } be generated by Algorithm 2, and x * is an optimal point of (4). Then either Algorithm 2 terminates in finite iterations, or x k satisfieslim k→∞ ∥g(x k )∥ = 0, lim k→∞ ∥x k -x * ∥ = 0.The next theorem characterizes the behavior of SSNS when the iterates are sufficiently close to the optimum. Since in Algorithm 2, the proposed new iterate x k + ξ k p k is rejected if ρ k ≤ 0, it is important to show that such rejections will be very rare as the algorithm proceeds. Theorem 4. There exists an integer K > 0 such that for all k ≥ K, ρ k ≥ 1 -ρ 0 , µ k+1 ≤ κ, andx k+1 = x k + ξ k p k .Theorem 4 indicates that eventually the parameter µ k is upper bounded by the pre-specified value κ, making the shift parameter of the Newton method, namely, µ k ∥g k ∥, tend to zero. This is consistent with the intuition we have described in Section 3.2. More importantly, the quantity ρ k will also be greater than the threshold 1 -ρ 0 , so the new step proposed by the Newton direction will always be accepted, regardless of the bounds c l and c u . 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 a quadratic local convergence rate that matches the Newton method based on genuine and dense Hessian matrices. Theorem 5 (Quadratic local convergence rate). Fix ξ k ≡ 1. Then there exists an integer K ′ > 0 and a constant L > 0 such that for all k ≥ K ′ ,∥x k+1 -x * ∥ ≤ L∥x k -x * ∥ 2 .Overall, Theorems 3 to 5 provide solid theoretical guarantees on the proposed SSNS algorithm, which are a major complement to the existing literature on solving entropic-regularized OT using realistic second-order methods. More importantly, they validate that SSNS is indeed a "safe" algorithm that is quite robust to the choice of initial value, step sizes, strength of sparsification, and hyperparameters.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it 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 not include experiments.• If the paper includes experiments, a 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 a 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 a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to 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, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example (a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm. (b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully. (c) If the contribution is a new model (e.g., a large language model), then there should either be a way to access this model for reproducing the results or a way to 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 some 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 it should be possible for other researchers to have some path to reproducing or verifying the results. 
numerical experiments: 
 In this section, we test the performance of the proposed SSNS algorithm on various numerical experiments 3 . There are a huge number of algorithms developed for entropic-regularized OT, and we focus on representative ones from each category of optimization methods: 1. the Sinkhorn algorithm as the default option for entropic-regularized OT, interpreted as a BCD method; 2. the APDAGD 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.• (Fashion-)MNIST: OT between a pair of images from the MNIST [15] or Fashion-MNIST [36] dataset. The a and b vectors are flattened and normalized pixel values, and the cost matrix holds the ℓ 1 -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 4 , 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 network, resulting in a 512-dimensional vector, then followed by a dimension reduction by principal component analysis. Let x i ∈ R 30 be the feature vector of an image in the first category, i = 1, . . . , n, and y j ∈ R 30 be the feature vector of an image in the second category, j = 1, . . . , m. Thena = n -1 1 n , b = m -1 1 m , and the cost matrix is M ij = ∥x i -y j ∥ 1 or M ij = ∥x i -y j ∥ 2 . The problem size is n ≈ m ≈ 1000.To make the regularization parameter η comparable in different settings, we normalize all cost matrices to have unit infinity norms, namely, M ← M/∥M ∥ ∞ . Then we consider two settings of the regularization parameter, η = 0.01 and η = 0.001.For entropic-regularized OT, a commonly-used criterion to evaluate optimality is the marginal error of the estimated transport plan. Let x k = (α T k , βT k ) T be the current iterate for dual variables, and setT k = τ (α k , β k ). Then the marginal error is given by ∥T k 1 m -a∥ 2 + ∥T T k 1 n -b∥ 2 .Coincidentally, this is exactly the current gradient norm ∥g(x k )∥ as indicated by (5) . Figure 1 shows the plots of marginal errors against iteration number or run time for different algorithms on the three benchmark datasets. From the top row of Figure 1 , it is clear that second-order methods have much faster convergence speed compared to 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 as shown in the bottom row of Figure 1 . 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 Fashion-MNIST data. The run time for SSNS is relatively longer in the ImageNet dataset, mostly because the transport plan is more dense under η = 0.01. We show in Figure 2 that as η becomes smaller, the Hessian matrix can be better approximated by a sparse matrix, thus enlarging the advantage of SSNS. This point is consistent with Theorem 1 of [33] , 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. We also note that with η = 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(x) is nearly singular, -[H(x)] -1 g(x) may not be a good search direction, and the introduction of the shift parameter λ in SSNS may give a better candidate -(H δ + λI) -1 g(x).In the ImageNet experiments above, we map each image to a 30-dimensional feature vector to compute the cost matrix. To study the impact of the feature dimension, in Section A.4, we conduct additional experiments with d = 60, 90, 200, 300, 500, under the same setting as 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 A.5 in the appendix.Finally, to study the scalability of SSNS, we consider the following synthetic OT problem that can generate data with arbitrary dimensions.• Large-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 • N (1, 0.2) + 0.8 • N (3, 0.5). We discretize the problem in the following way: let x i = 5(i -1)/(n -1), i = 1, . . . , n, and y j = 5(j -1)/(m -1), j = 1, . . . , 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 ãi = f 1 (x i ), bj = f 2 (y j ), a i = ãi / ( Similar to the previous experiment setting, we normalize the cost matrix, M ← M/∥M ∥ ∞ , and set η = 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 time-consuming to proceed. The results are visualized in Figure 4 , 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. Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: All datasets used are publicly available, and we have provided the programming code. Justification: This article is mostly a theoretical and computational work. Guidelines:• The answer NA means that there is no 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 or unintended uses (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 a direct path to 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 is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is 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.• If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
conclusion: 
 In this article, we have carefully designed and analyzed a new variant of the sparsified Newton method for solving large-scale entropic-regularized OT problems. As demonstrated by [33] , 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 on sparsified Newton method on the following aspects. First, a 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 [33] 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.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. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in 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 is 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:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL.• The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of 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, paperswithcode.com/datasets 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 it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
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/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • 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 an anonymized URL or include an anonymized zip file. 
 
 Appendix: All experiments in this article are conducted on a personal computer with an Intel i9-13900K CPU, 32 GB memory, and a Ubuntu 24.10 operating system.The experiments in Section 5 are designed to reflect two typical uses of OT. In the (Fashion-)MNIST experiment, images are vectorized as density vectors, and OT is a used as a tool for image morphing and interpolation [26] . We follow the experiment settings in the existing literature such as [27] and [33] to define a, b, and M : given two images, let a ∈ R 784 be the vectorized and normalized pixel values of the first image, and similarly define b ∈ R 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 j isM ij = |w -w ′ | + |h -h ′ | based on the ℓ 1 -distance, or M ij = (w -w ′ ) 2 + (h -h ′ ) 2base on the squared Euclidean distance. The whole cost matrix M ∈ R 784×784 then collects all pairwise cost values M ij .The ImageNet 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. Since each image is mapped to a 30-dimensional feature vector, the entries of the cost matrix are simply the ℓ 1 -distances or squared Euclidean distances between the feature vectors of two images.In Figures 1 and 2 , we have shown that the regularization parameter η may affect the performance of optimization algorithms. To further study this effect, we consider a series of equally-spaced η 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 in Table 1 . 1 we can find that BCD is very sensitive to the value of η. When η 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 η becomes smaller, the number of iterations also increases, but the overall runtime of SSNS may even decrease. This is because smaller η 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.In the ImageNet experiment in Section 5, we map each image to a 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. 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 ℓ 1 -distances to construct cost matrix, and the second one, demonstrated in Figure 7 , is based on squared Euclidean distances. Both test cases use the regularization parameter η = 0.001.In this proof all inequality signs hold elementwisely when applied to vectors and matrices.By Algorithm (1), T δ = T -∆, where each element ∆ ij of ∆ is either zero, or the corresponding element T ij of T . Clearly, we have ∆ ≥ 0. Due to the selection scheme in Algorithm (1), each row sum and each column sum of ∆ is upper bounded by the threshold δ, i.e., ∆1 m ≤ δ1 n and ∆ T 1 n ≤ δ1 m . Since the last column of ∆ is always a zero vector by the design of the algorithm, we also easily get ∆1 m-1 = ∆1 m ≤ δ1 n and ∆T 1 n = d ≤ δ1 m-1 , where d = ∆ T 1 n . D = H(x) -H δ = η -1 O ∆ ∆T O := D 1 D 2 ,so D ≥ 0, D 1 1 n+m-1 = η -1 ∆1 m-1 ≤ η -1 δ1 n and D 2 1 n+m-1 = η -1 ∆T 1 n ≤ η -1 δ1 m-1 . Overall, we have D1 n+m-1 ≤ η -1 δ1 n+m-1 . Since D ≥ 0, this is equivalent to ∥D i• ∥ 1 ≤ η -1 δ, i = 1, 2, . . . , n + m -1.Since D is symmetric, we also have∥D •j ∥ 1 = ∥D j• ∥ 1 ≤ η -1 δ, j = 1, 2, . . . , n + m -1. Let P k = j̸ =k |D kj |, k = 1, . . . , n + m -1,and note that all diagonal elements of D are zero. So by the Gershgorin circle theorem, every eigenvalue of D must be smaller than or equal to the maximum of D kk + P k = ∥D k• ∥ 1 , which is upper bounded by η -1 δ. Therefore, ∥D∥ = σ max (D) ≤ η -1 δ.Consider the matrix 1) ,A = diag(r) Tδ T T δ diag(c) ∈ R (n+m-1)×(n+m-and suppose that T = (T ij ), T δ = (S ij ). Clearly, S ij = T ij -∆ ij ≤ T ij . Let P k = j̸ =k |A kj |, k = 1, . . . , n + m -1,and then it is easy to find thatA kk + P k = r k + m-1 j=1 S kj ≤ 2r k , k = 1, . . . , n c k-n + n i=1 S i,k-n ≤ 2c k-n , k = n + 1, . . . , n + m -1 .Overall, we have A kk + P k ≤ U := 2 max{∥r∥ ∞ , ∥c∥ ∞ }. By the Gershgorin circle theorem, every eigenvalue of A must be smaller than or equal to the maximum value of A kk + P k , which is upper bounded by U < ∞.For some columns of ∆, say, column k, it may be the case that ∆ •k = 0, indicating that there is no thresholding on T •k . Note that if we simultaneously permute the columns of T and T δ , then the eigenvalues of A do not change. Therefore, without loss of generality we can assume that the first d columns and the last column of ∆ are exactly zero, 0 ≤ d ≤ m -1 (∆ •m = 0 always holds by the design of Algorithm 1). Then we can partition T δ as T δ = [T (1) , Tδ , T •m ], where T (1) ∈ R n×d and T •m ∈ R n have strictly positive elements. Accordingly, A can be partitioned asA =   diag(T 1 m ) T (1) T (2) δ T (1)T diag(T (1)T 1 n ) O T (2)T δ O diag(T (2)T δ 1 n )   .Consider the matrix B = A -svv T , where v ∈ R n+m-1 is a vector whose first (n + d) elements are ones and the remaining elements are zeros, and s is a positive scalar. DefineR k = j̸ =k |B kj |, k = 1, . . . , n + m -1,and suppose that s ≤ min 1≤i≤n,1≤j≤d T ij . Then for k = 1, . . . , n, it is easy to find thatEQUATIONThe inequality sign in ( 7) is strict, since we have assumed that at least one element of ∆ •,k-n is nonzero. Moreover,B kk =    m j=1 T kj -s, k = 1, . . . , n n i=1 T i,k-n -s, k = n + 1, . . . , n + d n i=1 T i,k-n , k = n + d + 1, . . . , n + m -1, so overall we can obtain thatB kk -R k =      m-1 j=1 ∆ kj + T km -(n -d)s, k = 1, . . . , n (n -d)s, k = n + 1, . . . , n + d n i=1 ∆ i,k-n , k = n + d + 1, . . . , n + m -1 . Now take s = min 1≤i≤n,1≤j≤m T ij 2n ,and then we can verify that s satisfies the condition s ≤ min 1≤i≤n,1≤j≤d T ij . Moreover, for k = 1, . . . , n, we haveB kk -R k ≥ T km -ns = T km - 1 2 • min i,j T ij ≥ 1 2 • min i,j T ij . For k = n + 1, . . . , n + d, B kk -R k = (n -d)s ≥ (n -m + 1)s = n -m + 1 2n • min i,jT ij .For k = n + d + 1, . . . , n + m -1, since ∆ i,k-n is either zero or T i,k-n, and at least one element of ∆ •,k-n is nonzero, we then haveB kk -R k = n i=1 ∆ i,k-n ≥ min 1≤i≤n T i,k-n ≥ min i,j T ij .Combining all the possibilities, we getB kk -R k ≥ L := n -m + 1 2n • min i,j T ij > 0for all k and d.By the Gershgorin circle theorem, every eigenvalue of B must be greater than or equal to the minimum value of B kk -R k , which is lower bounded by L > 0. Since A = B + svv T and vv T is positive semi-definite, we have σ min (A) ≥ L > 0, implying that A is positive definite.We first introduce some necessary notation. Define the level set as L(x 0 ) = {x : f (x) ≤ f (x 0 )}, and letf * = inf{f (x) : x ∈ L(x 0 )}, β 1 = sup{∥g(x)∥ : x ∈ L(x 0 )}, β 2 = sup{∥H δ (x)∥ : x ∈ L(x 0 ), δ > 0}. Define the function Φ(x, d) = f (x + d) -f (x) -g(x) T d.We then present a few technical lemmas, followed by the proof of the main theorem. Lemma 6. L(x 0 ) is a bounded and closed convex set, f * > -∞, β 1 < ∞, and β 2 < ∞.Proof. Clearly, f (x) is a continuously differentiable and strictly convex function. Define L c = {x : f (x) ≤ c}, and then L c is a closed convex set.Lemma 1 of Dvurechensky et al. [10] shows that if (α * , β * ) is an optimal solution to (3), thenmax j β * j -min j β * j ≤ R,where R = ∥M ∥ ∞ /η -log(min i,j {a i , b j }). Since we have fixed β * m = 0, we immediately obtain that ∥β * ∥ < ∞. Given β * , α * has a closed form,α * i = η log a i -η log   m j=1 e η -1 (β * j -Mij )   ,so ∥α * ∥ is also bounded.Let x * = (α * T , β * T ) T and f * = f (x * ), and then clearly f * > -∞. Take c = c * = {f * }, and obviously, L c * is non-empty and bounded. Then Corollary 8.7.1 of [30] shows that L c is bounded for every c, implying that L(x 0 ) = L f (x0) is bounded. Since g(x) is continuous on L(x 0 ) and L(x 0 ) is compact, it is easy to show that β 1 < ∞.By Theorem 2, ∥H δ ∥ = σ max (H δ ) ≤ 2η -1 • max{∥r∥ ∞ , ∥c∥ ∞ }, where r = T 1 m and c = T T 1 n .Clearly, every element of T is bounded if x is bounded, making ∥r∥ ∞ and ∥c∥ ∞ also bounded. As a result, ∥H δ ∥ is bounded on L(x 0 ), implying that β 2 < ∞.Lemma 7. For any ε > 0, there exists a constant δ 0 > 0 such that for all η ∈ [c l , c u ], x ∈ L(x 0 ), and ∥p∥ ≤ δ 0 satisfying x + ηp ∈ L(x 0 ), we have|Φ(x, ηp)| ≤ ρ 0 c l ε 4 ∥p∥.Proof. Since f is continuously differentiable on L(x 0 ), by the mean value theorem for multivariate functions, there exists some ψ ∈ (0, 1) such that for all x ∈ L(x 0 ) and x + ηp ∈ L(x 0 ),f (x + ηp) -f (x) = g(x + ψηp) T (ηp).Therefore,Φ(x, ηp) = g(x + ψηp) T (ηp) -g(x) T (ηp) = [g(x + ψηp) -g(x)] T (ηp).By Lemma 6, g is continuous on L(x 0 ) and L(x 0 ) is a compact set, so g is uniformly continuous on L(x 0 ). Therefore, for any ε > 0, there exists a constant δ 0 > 0 such that whenever ∥p∥ ≤ δ 0 , we have∥g(x + ψηp) -g(x)∥ ≤ ρ 0 c l ε 4c u . As a result, |Φ(x, ηp)| ≤ ∥g(x + ψηp) -g(x)∥ • ∥ηp∥ ≤ ρ 0 c l ε 4 ∥p∥.Lemma 8. Let {x k } be generated by Algorithm 2, and thenm k (0) -m k (ξ k p k ) ≥ 1 2 ∥g k ∥ • min ∥ξ k p k ∥, ∥g k ∥ ∥H δ k ∥ . Proof. Let f k = f (x k ) and r k = ∥ξ k p k ∥. If we take (p, λ) = (ξ k p k , µ k ∥g k ∥), and then we can verify that (p, λ) meets the following relations:       g k + H δ k p + λp = 0 λ ≥ 0 r k -∥p∥ ≥ 0 λ • (r k -∥p∥) = 0.As m k (•) is a convex function, we have that ξ k p k is a KKT point and an optimal solution to the constrained optimization problemEQUATIONNext, consider the Cauchy point defined in Chapter 4 of Nocedal and Wright [23] ,p k c = -τ k r k ∥g k ∥ g k ,whereτ k = 1, if g T k H δ k g k ≤ 0 min 1, ∥g k ∥ 3 r k g T k V k g k , if g T k H δ k g k > 0.Then we can obtain the following inequality by Lemma 4.3 of Nocedal and Wright [23] :m k (0) -m k (p k c ) ≥ 1 2 ∥g k ∥ • min r k , ∥g k ∥ ∥H δ k ∥ .Meanwhile, the Cauchy point p k c is a feasible point for the constrained problem (8), and ξ k p k is an optimal solution for this problem. Therefore, m k (p k c ) ≥ m k (ξ k p k ), and we havem k (0) -m k (ξ k p k ) ≥ 1 2 ∥g k ∥ • min r k , ∥g k ∥ ∥H δ k ∥ .Lemma 9. Let g k and µ k be generated by Algorithm 2. If there exists a constant ε > 0 and an integer K such that ∥g k ∥ ≥ ε for all k ≥ K, then there must exist a sufficiently large constant μ > 0, such that µ k+1 ≤ 4μ for all k ≥ K.Proof. For a given μ, defineI 1 = {k : k ≥ K, µ k < μ} and I 2 = {k : k ≥ K, µ k ≥ μ}.If for some μ, I 2 is finite, then it trivially holds that all µ k has a global upper bound, which leads to the desired conclusion. Therefore, we only consider the case that I 2 is infinite for all μ > 0. Since μ can be chosen arbitrarily, {µ k } k∈I2 must be unbounded.Then we will estimate the bound of the following quantity,|ρ k -1| = m k (ξ k p k ) -f (x k + ξ k p k ) m k (0) -m k (ξ k p k ) .According to Lemma 6, one has ∥g k ∥ ≥ ε and∥g k ∥ ≤ β 1 for k ≥ K. Theorem 2 shows that H δ k is positive definite, so σ min (H δ k + µ k ∥g k ∥) ≥ µ k ε, (H δ k + µ k ∥g k ∥I) -1 ≤ µ -1 k ε -1 . By Algorithm 2, we have ξ k p k = ξ k (H δ k + µ k ∥g k ∥I) -1 g k ,and hence for all k ∈ I 2 ,EQUATIONEquation ( 9) indicates that ∥ξ k p k ∥, k ∈ I 2 can be made arbitrarily small with a sufficiently large μ. Therefore, we can choose some μ such thatEQUATIONwhere δ 0 is defined in Lemma 7. Meanwhile, we know|m k (ξ k p k ) -f (x k + ξ k p k )| = Φ(x k , ξ k p k ) - 1 2 (ξ k p k ) T H δ k (ξ k p k ) ≤ |Φ(x k , ξ k p k )| + 1 2 ∥H δ k (ξ k p k )∥ • ∥ξ k p k ∥ ≤ |Φ(x k , ξ k p k )| + c u 2 ∥H δ k ∥ • ∥ξ k p k ∥ • ∥p k ∥.Lemma 7 indicates that when ∥ξ k p k ∥ ≤ δ 0 , we have|Φ(x k , ξ k p k )| ≤ ρ 0 c l ε 4 ∥p k ∥.Therefore, with p k that satisfies (10) , we obtain|m k (ξ k p k ) -f (x k + ξ k p k )| ≤ ρ 0 c l ε 4 ∥p k ∥ + c u 2 • β 2 • ρ 0 c l ε 2c u β 2 • ∥p k ∥ = ρ 0 c l ε 2 ∥p k ∥.On the other hand, Lemma 8 indicates thatm k (0) -m k (ξ k p k ) ≥ 1 2 ∥g k ∥ • min ∥ξ k p k ∥, ∥g k ∥ ∥H δ k ∥ ≥ 1 2 ε • min ∥ξ k p k ∥, ε β 2 = 1 2 ε • ∥ξ k p k ∥ ≥ c l ε 2 • ∥p k ∥,so we have|ρ k -1| = m k (ξ k p k ) -f (x k + ξ k p k ) m k (0) -m k (ξ k p k ) ≤ ρ0c l ε 2 ∥p k ∥ c l ε 2 • ∥p k ∥ = ρ 0 , which means that ρ k ≥ 1 -ρ 0 for all k ∈ I 2 .Then by the design of Algorithm 2, we haveµ k+1 ≤ µ k for all k ∈ I 2 .Now we can show that µ k+1 ≤ 4μ for all k ≥ K by induction. First, enlarge μ when necessary to ensure that μ > µ K . Then by Algorithm 2, we must have µ K+1 ≤ 4μ. Now suppose that µ l+1 ≤ 4μ for some l ≥ K. If µ l+1 ∈ I 1 , then clearly µ l+2 ≤ 4μ immediately holds. Otherwise, µ l+1 ∈ I 2 , so by the argument above, we have µ l+2 ≤ µ l+1 ≤ 4μ. In both cases, the conclusion holds.Lemma 10. Under the same conditions as in Lemma 9, defineK = {k : k ≥ K, ρ k ≥ ρ 0 } .Then K must be a finite set.Proof. We use proof by contradiction. Suppose that K is an infinite set. Sinceρ k = f (x k ) -f (x k + ξ k p k ) m k (0) -m k (ξ k p k ) ,and note thatx k+1 = x k if ρ k < 0, we have f (x k+1 ) ≤ f (x k + ξ k p k ).Therefore, by Lemma 8, it holds thatf (x k ) -f (x k+1 ) ≥ ρ 0 [m k (0) -m k (ξ k p k )] ≥ ρ 0 2 ∥g k ∥ • min ∥ξ k p k ∥, ∥g k ∥ ∥H δ k ∥ ≥ ρ 0 2 ε • min ∥ξ k p k ∥, ε β 2for all k ∈ K. Combined with Lemma 6, we find that f (x k ) is monotonically non-increasing with a lower bound, so f (x k ) has a limit, and hence lim k∈K,k→∞∥ξ k p k ∥ = 0.On the other hand,p k = -(H δ k + µ k ∥g k ∥I) -1 g k , which means that g k = -(H δ k + µ k ∥g k ∥)p k . Therefore, ε ≤ ∥g k ∥ = ∥(H δ k + µ k ∥g k ∥I)p k ∥ ≤ ∥H δ k + µ k ∥g k ∥∥ • ∥p k ∥ ≤ [∥H δ k ∥ + µ k ∥g k ∥] • ∥p k ∥ ≤ (β 2 + µ k β 1 ) • 1 c l • ∥ξ k p k ∥.In other words,0 ≤ c l ε β 2 + µ k β 1 ≤ ∥ξ k p k ∥ → 0, which implies that µ k → ∞ for k ∈ K, k → ∞.This contradicts with the fact that µ k ≤ 4μ as shown in Lemma 9. So to conclude, K must be a finite set.Lemma 11. Let {x k } be generated by Algorithm 2. Then either Algorithm 2 terminates in finite iterations, or g(x k ) satisfies lim inf k→∞ ∥g(x k )∥ = 0.Proof. We prove this lemma by contradiction. Suppose that there exist some ε > 0 and an integer K such that ∥g k ∥ > ε, ∀k ≥ K. Then Lemma 10 implies that the index setK = {k : k ≥ K, ρ k ≥ ρ 0 }is infinite. This means that there is a sufficiently large integer K ′ such that ρ k < ρ 0 for all k ≥ K ′ . According to Algorithm 2, we must haveµ k+1 = 4µ k , ∀k ≥ K ′ ,which means that µ k → ∞. However, this contradicts with the fact that µ k ≤ 4μ for some μ > 0 as shown in Lemma 9. Therefore, we must havelim inf k→∞ ∥g k ∥ = 0.Then we are ready to prove Theorem 3. Let x * be the unique global optimum of (4). Clearly, f (x) and g(x) are Lipschitz continuous on L(x 0 ), so there exist constantsEQUATIONEQUATION)On the other hand, by taking δ = 0, Theorem 2 shows that there is a constant c 1 > 0 such that σ min (H(x)) ≥ c 1 for all x ∈ L(x 0 ), so there is a constant c 2 > 0 such thatEQUATIONAlso by Taylor's theorem,f (x k ) = f (x * ) + [g(x * )] T (x k -x * ) + 1 2 (x k -x * ) T H(z 1 )(x k -x * ),where z is some point between x k and x * . Since g(x * ) = 0, we obtainEQUATION)By the design of Algorithm 2, f (x k ) is non-increasing and is lower bounded by f * = f (x * ), so f (x k ) must have a limit. Suppose thatlim k→∞ |f (x k ) -f * | = ε ≥ 0,and then (11) and (13) indicate that for sufficiently large k,1 2 C -1 1 ε ≤ C -1 1 |f (x k ) -f (x * )| ≤ ∥x k -x * ∥ ≤ c -1 2 ∥g(x k )∥.Since in Lemma (11) we have shown that lim inf k→∞ ∥g k ∥ = 0, ε cannot be any positive value. Therefore, ε = 0. Combining ( 12) and ( 14), we havec 1 2C 2 2 ∥g k ∥ 2 ≤ c 1 2 ∥x k -x * ∥ 2 ≤ f (x k ) -f (x * ) → ε = 0,which implies that lim k→∞ ∥g k ∥ = 0 and lim k→∞ ∥x k -x * ∥ = 0.For convenience, letH k = H(x k ) and B k = H δ k + µ k ∥g k ∥I.Theorem 2 shows that for any δ,σ min (H δ ) ≥ η -1 • n -m + 1 2n • min i,j T ij . Note that T = τ (α, β) = e η -1 (αi+βj -Mij ) ,and each (α i + β j ) is bounded on L(x 0 ), so each element of T must be bounded below from zero on L(x 0 ). Therefore, there exists a constant 0 < C 1 < ∞ such that σ min (H δ ) ≥ C -1 1 for all x ∈ L(x 0 ).By the design of Algorithm 2, x k ∈ L(x 0 ) for all k. As a result,EQUATIONand hence ∥p k ∥ = ∥B -1 k g k ∥ ≤ C 1 ∥g k ∥ for all k. Moreover, Theorem 2 shows that ∥H δ k ∥ = σ max (H δ k ) ≤ C 2 for some C 2 > 0, so by Lemma 8, we havem k (0) -m k (ξ k p k ) ≥ 1 2 ∥g k ∥ • min ∥ξ k p k ∥, ∥g k ∥ ∥H δ k ∥ ≥ 1 2C 1 ∥p k ∥ • min c l ∥p k ∥, ∥p k ∥ C 1 C 2 ≥ min{c l , C -1 1 C -1 2 } 2C 1 • ∥p k ∥ 2 := C 3 ∥p k ∥ 2 .On the other hand,EQUATIONBy Taylor's theorem,f (x k + ξ k p k ) -f (x k ) -g T k (ξ k p k ) - 1 2 (ξ k p k ) T H k (ξ k p k ) T = o(∥ξ k p k ∥ 2 ).For the second term of ( 16), let D k = H k -H δ k , and then Theorem 1 shows that∥D k ∥ ≤ η -1 δ k . By Algorithm 2, δ k ≤ ν 0 ∥g k ∥ γ , so ∥H k -H δ k ∥ ≤ η -1 ν 0 ∥g k ∥ γ . Then 1 2 (ξ k p k ) T (H k -H δ k )(ξ k p k ) T ≤ ν 0 2η ∥g k ∥ γ • ∥ξ k p k ∥ 2 = o(∥ξ k p k ∥ 2 ). As a result, |m k (ξ k p k ) -f (x k + ξ k p k )| = o(∥ξ k p k ∥ 2 ), and hence|ρ k -1| = m k (ξ k p k ) -f (x k + ξ k p k ) m k (0) -m k (ξ k p k ) ≤ o(∥ξ k p k ∥ 2 ) C 3 ∥p k ∥ 2 → 0.This implies that there is an integer K > 0 such that for all k ≥ K, ρ k ≥ 1 -ρ 0 . By the design of Algorithm (2), we have µ k+1 ≤ κ and x k+1 = x k + ξ k p k for all k ≥ K.We first present a classical result derived from the mean value theorem of vector-valued functions. Lemma 12 (Theorem 3.2.12 of [25] ). Let g : D ⊂ R n → R m be continuously differentiable on a convex set D 0 ⊂ D and suppose that for constants α ≥ 0 and p ≥ 0, ∇g satisfies∥∇g(u) -∇g(v)∥ ≤ α∥u -v∥ p , ∀u, v ∈ D 0 .Then, for any x, y ∈ D 0 ,EQUATION)Then we are ready to prove Theorem 5. Let x * be the unique global optimum of (4), and then g(x * ) = 0. Clearly, g(x) and H(x) are Lipschitz continuous on L(x 0 ), so there exist constants L 1 , L 2 > 0 such that for all x, y ∈ L(x 0 ),∥g(x) -g(y)∥ ≤ L 1 ∥x -y∥, ∥H(x) -H(y)∥ ≤ L 2 ∥x -y∥.This implies that g satisfies the conditions in Lemma 12 with α = L 2 and p = 1. By substituting x ← x k and y ← x * , we have∥g(x k )∥ = ∥g(x k ) -g(x * )∥ ≤ L 1 ∥x k -x * ∥, ∥g(x * ) -g(x k ) -H(x k )(x * -x k )∥ = ∥g(x k ) -g(x * ) -H(x k )(x k -x * )∥ ≤ 1 2 L 2 ∥x k -x * ∥ 2 .Without loss of generality let K ′ > K, and then for all k ≥ K ′ , µ k ≤ κ and x k+1 = x k + p k . In the proof of Theorem 5, we have shown that∥H k -H δ k ∥ ≤ η -1 ν 0 ∥g k ∥ γ , so for ∥g k ∥ ≤ 1, ∥H k -B k ∥ ≤ ∥H k -H δ k ∥ + µ k ∥g k ∥ ≤ (η -1 ν 0 + κ)∥g k ∥ ≤ L 1 (η -1 ν 0 + κ)∥x k -x * ∥.Also note that ∥B -1 k ∥ ≤ C 1 from ( 15), and then we have∥x k+1 -x * ∥ = ∥x k + p k -x * ∥ = ∥x k -x * -B -1 k g k ∥ = ∥B -1 k (B k (x k -x * ) -g k )∥ ≤ ∥B -1 k ∥ • ∥g k -g(x * ) -H k (x k -x * ) + H k (x k -x * ) -B k (x k -x * )∥ ≤ C 1 (∥g k -g(x * ) -H k (x k -x * )∥ + ∥(H k -B k )(x k -x * )∥) ≤ C 1 1 2 L 2 ∥x k -x * ∥ 2 + ∥H k -B k ∥ • ∥x k -x * ∥ ≤ C 1 1 2 L 2 ∥x k -x * ∥ 2 + L 1 (η -1 ν 0 + κ)∥x k -x * ∥ 2 = C 1 1 2 L 2 + L 1 (η -1 ν 0 + κ) ∥x k -x * ∥ 2 .• 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 crossreferenced. • All assumptions should be clearly stated or 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.https://github.comm/yixuan/regot-python.The programming code to reproduce the results is available at https://github.com/TangZihao1997/ SSNS.https://github.com/fastai/imagenette. 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: I have not found any major weaknesses. Minor comments/typos:
1.	The 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?
2.	The 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. 
3.	Line 67: “proved” should be “proven.”
4.	Line 110: The full name of “BFGS” is missing.
5.	What 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.
questions: 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: limitations: The authors have mentioned the limitation of their work.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.
questions: questions: See weaknesses.

- Following the above, is it possible to provide a runtime in the style of Sinkhorn given a prescribed tolerance level?
- Can this approach be adapted for unbalanced entropic OT?
limitations: limitations: N/A
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

2. 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.

3. 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.
questions: questions: 1. What is the computational complexity of Algorithm 1 Sparsification?
2. What is the computational complexity of each iteration in Algorithm 2?
3. Can you compare with the SNS method and report results in tables during rebuttal?
4. 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?
5. Can you also replicate the Table 3 and Table 4 in the SNS paper? Basically conduct perturbation study on the entropy regularization parameters.
6. 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: limitations: NA
paper_summary: paper_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.
strength_weakness: strength_weakness: (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.
questions: 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.

I have the following comments and questions:

1. 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".

2. 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?

3. 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.

4. In all images, Newton method seems to be very slow. Is this all due to the costly process of solving the dense linear system?

5. 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: limitations: The authors are clear about their limitations.

Comments: comment: Dear reviewer,

We really appreciate your comments and totally understand your concern. To avoid violating the double blind policy, we have communicated with the Area Chairs to handle this issue. What we can comment here is that we have tried our best to enhance the reproducibility and integrity of this work, and the comparison with [31] is subject to some practical difficulties. Given this situation, we have added similar experiments in [31] per your suggestions, which are contained in the global rebuttal.

Thanks for your understanding.
Thank you very much for your reply.

Most of my questions have answered, except the requested comparison with [31], which is the following work:

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

This might sound a little bit too harsh, but I have some leftover question on whether there are any author overlaps between this submission and [31].
If there is an author overlap, then you defintely have access to the code base in [31].
Some research groups do this so that they can publish another paper with impressive results, without having to beat against their previous method.
I apologize if this sounds too mean.

However, I wouldn't know this because of NeurIPS's double blind policy. I will give you the benefit of doubt and trust you that there is no author overlap and you do not have access to the codebase of [31]. I do want to leave this comment here, just in case this possibility happens.

As for now, I will maintain my score since it is already above the acceptance threshold.
Thanks for the comments. I will still keep my score as-is.
Thanks for the feedback!
I thank the authors for their rebuttal. I decide to keep my score as I think this is a good work and the authors do clarify most of my concerns.
Thank you for the clarifications!