The Sparse-Plus-Low-Rank Quasi-Newton Method for Entropic-Regularized Optimal Transport
Authors: Chenrui Wang
Yixuan Qiu
Venue: ICML 2025 poster (2025)
Keywords: optimal transport
entropic regularization
quasi-Newton method
Hessian sparsification
low-rank approximation
convergence analysis
Primary area: optimization->convex
TL;DR: A new quasi-Newton method to solve entropic-regularized optimal transport, which combines the advantages of Hessian sparsification and low-rank approximation, with provable convergence properties.

Abstract: The entropic-regularized optimal transport (OT) has gained massive attention in machine learning due to its ability to provide scalable solutions for OT-based tasks. However, most of the existing algorithms, including the Sinkhorn algorithm and its extensions, suffer from relatively slow convergence in many cases. More recently, some second-order methods have been proposed based on the idea of Hessian sparsification. Despite their promising results, they have two major issues: first, there is limited theoretical understanding on the effect of sparsification; second, in cases where the transport plan is dense, Hessian sparsification does not perform well. In this paper, we propose a new quasi-Newton method to address these problems. First, we develop new theoretical analyses to understand the benefits of Hessian sparsification, which lays the foundation for highly flexible sparsification schemes. Then we introduce an additional low-rank term in the approximate Hessian to better handle the dense case. Finally, the convergence properties of the proposed algorithm are rigorously analyzed, and various numerical experiments are conducted to demonstrate its improved performance in solving large-scale OT problems.

The Sparse-Plus-Low-Rame Quasi-Newton Method for Entropic-Regularize Optimal Transport
Chenrui Wang' Yixuan Qiu
Abstract
with statistical modeling and machine learning tasks. In particular, OT provides a principled way to measure the similarity between probability distributions by considering the cost of transporting the mass between them (Villani et al., 2009), making it highly relevant in applications involving structured data or distributions with geometric properties. See Torres al. (2021); Montesuma al. (2024) for an overview of its applications in modern machine learning.
The entropic-regularee optimal transport (OT) has gained massive attention in machine learn- ing due to its ability to provide scalable solutions for OT-based tasks. However, most of the ex- isting algorithms, including the Sinkhorn algo- rithm and its extensions, suffer from relatively slow convergence in many cases. More recently, some second-orde order methods have been proposed based on the idea of Hessian sparsification. De- spite their promising results, they have two major issues: first, there limited theoretical under- standing on the effect of sparsification; second, in cases where the transport plan is dense, Hessian sparsification does not perform well. In this pa- per, we propose new quasi-Newton method to address these problems. First, we develop new theoretical analyses to understand the benefits of Hessian sparsification, which lays the foundation for highly flexible sparsification schemes. Then we introduce an additional low-rank term in the approximate Hessian to better handle the dense case. Finally, the convergence properties of the proposed algorithm are rigorously analyzed, and various numerical experiments are conducted to demonstrate its improved performance in solving large-scale OT problems.
The discrete OT problem can be characterized by the fol- lowing linear programming problem:
min (P,M), P
[Ta b) = (P R Raxm Plm = a, 1T, = b, D2P200 0}:
where E R" , E Rm are two vectors satisfying a a²', f²1, 1, a > 0, b 0, M is given cost matrix, and all incquality signs are elementwise.
More recently, the application of OT has experienced a significant boost due to the development of approximate solvers such as the entropic-regularized OT. Popularized by Cuturi (2013), the entropic-regularized OT incorporates an entropic regularization h(P) ELEE Psy
log P.j) into the OT problem:
(P. nh(P). PELI(a,b)
min
(I)
This regularization significantly reduces the computational cost of OT based on the Sinkhorn Sinkhorn-Knopp algorithm (Yule, 1912; Sinkhorn, 1964), thus unlocking its potential in large- scale problems. Along this direction, many extensions of the Sinkhom algorithm have also been proposed (Altschuler et al., 2017: Dvurechensky et al., 2018; Guminov et al., 2021; Lin et al., 2022). As a result, OT is increasingly applied to solve a wide range of challenges in fields such as image processing, graphics, and machine learning (Peyré et al., 2019).
1. Introduction
Optimal transport (OT), as originally described by Gaspard Monge, addresses the problem moving a distribution (e.g.. s.. pile of sand) match target configuration (e.g., prescribed shape) while minimizing cost, such the total distance effort required. In recent years, OT has received significant attention, due to its strong connections
'School of Statistics and Data Science, Shanghai University of Finance and Economics, Shanghai, China. Correspondence to: Yixuan Qiu <qiuyixuan sduc.e.
However, the computation of problem (1) is still major challenge. For example, the Sinkhorn algorithm generally requires large number of iterations to converge, especially when the regularization parameter small. Other first-order methods that extend the Sinkhorn algorithm also show rela- tively slow convergence. To this end, another approach to
Proceedings of the nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

The Dparse-Plus-Loww-Ratt Quasi-Newton Method
solving entropiiereguularized OT is to study the dual prob- lem of (1), ie,, maximizing a function L(a,8), a E IR", 3 E R where
Hessian matrix is relatively dense. Rigorous convergence analysis is provided support the application of the pro- posed algorithm in large-scale OT problems. An efficient implementation of the method is included in the RegOT Python package'.
L(a, B) = TTa a + 3Tb
n m
exp(1 (pp(1-((;10 (a; + -3-- M,y)). (2)
11
Our main contribution compared to prior art is summarized as follows:
Contribution
Once an optimal solution (a", to problem (2) is sought, the primal optimal solution (1), denoted by T _=xm can be obtained as Tÿ- exp(" (af + +3 MM2))) Since ()), C(a + cl., 3 clm) for all E R. can remove the redundant degree of freedom by setting 3m = 0 globally. Then eventually, solving entropic-regularized OT is equivalent to solving
1. New theoretical results are developed to understand the mechanism of Hessian sparsification. Such a theoret- ical understanding also guarantees that broad range of sparsification schemes enjoy desirable properties.
2. A new quasi-Newton method is proposed to solve entropic-regularized OT that combines the advantages of Hessian sparsification and low-rank approximation, achieving fast convergence speed with low computa- tional cost.
min EE
(w)
(3)
where
3. We provide convergence guarantees for the proposed method, and conduct extensive numerical experiments to demonstrate its performance.
f(x) = L( =
expl11 (a: + -B)-
i=1 j=1
a T a-- BTb,
Notation Throughout this article we adopt the following notation. For n E N. denote [n] == {1.2...., ,n). a matrix A E Rmxm and a vector v E R we use Ā E Rnx(m-1) 1) to represent the first (m - 1) columns of A, and E R=-1 to represent the first (m 1) elements of u; inequality signs such as A > and v 0 are all elementwise. We use Amax and Xmia(-) represent the largest and smallest eigenvalues of real symmetric matrices, respectively. For matrix A, let A, be the vector of the i-th row of A, and A, be the vector of the j-th column of A.
and x = (o1...., On. .,...., Bm- Clearly, (3) a smooth and unconstrained convex optimization problem, which brings the possibility to use second-order methods for acceleration.
Given
However, the classical Newton method is not a realistic approach, since the Hessian matrix of of fx) is a dense matrix of the size (n + m - 1) x (n + m - 1). leading to an O(n + m) computational cost in computing the Newton direction. More recently, Tang et al. (2024) and Tang & Qiu (2024) tackle this problem based on the idea of Hessian sparsification, i.e., approximating the true Hessian matrix by sparse matrices, thus leading to efficient sparse linear systems to compute the search directions. This line of works show some promising results, but two major issues remain to be solved: first, there is limited theoretical understanding on the effect of sparsification; second, in cases where the true Hessian matrix is relatively dense, Hessian sparsification does not perform well.
2. Related Work
OT in machine learning Optimal transport has emerged as a powerful mathematical framework with diverse applica- tions in machine learning (Torres et al., 2021; Montesuma et al., 2024), is widely used for tasks such as domain adap- tation (Courty et al., 2017), generative modeling (Arjovsky et al., 2017; Genevay et al., 2018; Huynh et al., 2021), clus- tering (Laclau et al., 2017), and cross-domain alignment (Chen al., 2020), among many others.
In this paper, we mainly target on resolving the two issues above, and propose a new quasi-Newton method to solve entropic-regularized OT. First, we provide new theoreti- cal analyses on the sparsified Hessian matrices, and show that sparsification brings various benefits. The results also motivate a class of flexible sparsification schemes that sub- stantially generalize existing methods. Second, we propose a new model to better approximate the true Hessian ma- trix, which combines Hessian sparsification with low-rank approximation. We demonstrate that this method greatly en- hances the algorithm's performance in scenarios where the
Solving entropie-regularizee OT There is a rich collec- tion of algorithms developed to solve entropic-regularized OT. Following the seminal work Cuturi (2013) that popular- izes the Sinkhorn-Knopp algorithm (Yule, 1912; Sinkhorn, 1964), many extension methods are proposed (Altschuler et al., 2017: Dvurechensky et al., 2018; Guminov et al., 2021: Lin et al., 2022). In addition to first-order methods,
hhttps://githin com/yi http:/yixxanrrrottooeoo
2

---

The Plarss--luus-LowRaann Rank Quasi-Newton Method
Brauer et al. (2017); Tang et al. (2024); Tang & Qiu (2024) consider Newton-type second-order methods to solve the dual problem of entropic-regulariized OT. Another direction, as in the importance sparsification approach (Li et al., 2023), accelerates computation by constructing a sparse approxi- mation of the kernel matrix, thereby significantly reducing the cost of Sinkhorn iterations.
cost at the order of O + m))), which is too demanding for large-scale OT problems. Therefore, several Hessian sparsification methods have been proposed to approximate H by some sparse matrices, which lead to significantly reduced computational costs.
The SNS algorithm (Tang et al., 2024) sparsifies the Hessian matrix by a thresholding rule. Specifically, any entry in the Hessian matrix smaller than constant p is truncated to zero. However, although this procedure preserves the symmetry and diagonal dominance of the Hessian matrix, there is no guarantee on the positive definiteness of the sparsified Hessian. Specifically, 0 might be included in one of the Gershgorin discs, admitting À = 0 as a possible eigenvalue.
Quasi-Newton methods Quasi-Newton methods are a class of optimization techniques designed to efficiently solve large-scale smooth optimization problems (Nocedal & Wright, 2006). Unlike the classical Newton methods, they approximate the Hessian matrix by simpler structures in computing the search directions, thus reducing computa- tional complexity while retaining rapid convergence. These features make them highly effective for machine learning and optimization tasks involving high-dimensional data. For example, Cuturi & Peyré (2018) suggests using the limited- memory Broyden-FFltcher -oodaeee-eeee--oofarbb------ Shanno (L- BFGS) method (Liuu & Nocedal, 1989) to solve the dual problem of entropic-regularize OT.
The SSNS algorithm (Tang & Qiu, 2024) takes similar but different approach, which only sparsifies the off-diagonal elements of the Hessian matrix. Meanwhile, it controls the row-wise and column-wis wise approximation errors. This method provides guarantees on the positive definiteness, but it is not flexible to design the sparsity pattern. In particular, the density after sparsification is unknown in advance, thus making it difficult to control the computational cost. To this end, in Section 3.2 we carefully analyze the eigenvalue structure of the sparsified Hessian, which provides various new insights on the consequence of Hessian sparsification.
3. Understanding Hessian Sparsification
3.1. Background
Recall that our main objective is to solve the unconstrained convex optimization problem (3). is known that fx) has some good properties:
3.2. Eigenvalue structure of the sparsified Hessian
As can be seen from (4), the true Hessian matrix at (aT, )T) is determined by the matrix T = T(a, B3. There- fore, a natural way to sparsifying Hxx) is to sparsify the T matrix. Following the method proposed in Tang & Qiu (2024), we consider the off-diagonal sparsification frame- work as given in Definition 3.1. Definition 3.1 (Sparsification scheme). A sparsification scheme is defined by a set of coordinates Ω S l = {(i,)) i E [m].] E [m 1]]. In particular, the sparsified matrix Ta has elements
1. (xx is strictly convex, so if (3) has a solution, then the solution is unique.
2. The gradient of fxx), denoted by g(x) = V (z), has a closed-form expression:
T1 m
g(x) =
T TTa 3),
7 1,
b)
where the free variables are x = (a", B,,,, 3 = T(a, 3) is an n x m matrix with elements [Taa, F] = exp(1) + 3, - M,y)], and recall that T means removing the last column of T.
[Tg. (i.j) Ω 10. (1,5) $ Ω
(Ta)
and the sparsified Hessian matrix is given by
He = Ha(x) 1 [diag(T 1m)
T diag (T 1,))
3. fxx) is twice differentiable, and the Hessian matrix H(x) fx:) also has simple expression:
/iag(T7 (7'm) TT
T
Clearly, Ha depends on the current variable x, but we will omit it for brevity if no confusion is caused. It also worth noting that the diagonal elements of Ha remain unchanged compared to H Hxx), which are computed from the original T instead of the sparsified Ti. This structure is crucial for our theoretical analysis.
H(x)= = -1
(4)
diag (TT 1,))
The classical Newton method solves (3) by generating sequence of iterates } based on the update rule Ik+1 = In - a. 2 9k. where 9& = g(xx). H₂ = H (xx), and ag > 0 is the step size at iteration k. However, solving HE " for a dense Hessian matrix Hg has a computational
To gain insights on the effect of sparsification, we first con- sider the process of incremental sparsification: removing
3

---

The Charrs-Pluus-Low-Raaee Quasi-Newton Method
exactly two symmetric elements from the current approxi- mate Hessian Hn.- resulting in another matrix Ho,, where Go and Qi only differ by one element, and Q2 C Ro. Then the main theoretical result of this section, Theorem 3.3, claims that Ha, strictly decreases the condition number of Ha.
The results clearly demonstrate that during the incremental sparsification process, the minimum eigenvalue monoton- ically increases, while the maximum eigenvalue monoton- ically decreases, which is consistent with our theoretical predictions.
Hesslan Elgenvalues (n - 30, m 30)
Assumption 3.2. For sparsification scheme Q C 92 there exists positive integer p > 0 such that the p-th power of Ho has strictly positive entries, i..., (Hg)) > 0.
Hn Egenvalue
Haz Dignrwuue
25
Egenvalue
20
Given two sparsification schemes Ωo, St C G, suppose that 22 c Do and they only differ by ele- ment. If Assumption 3.2 holds for St then we have
Theorem 3.3.
15
10
Amax(Has) < Amax (Hsa): Amin (Han) > Amin (Hsg))
100
200
300
400
500
600
700
800
Sparsification Step
Hessian Condition Number (n = 30, 00.m= = 30)
which implies that Ha, strictly decreases the condition number of Hao
10 Condition Number
2.0
Logi0 Cundition Number
1.5
Suppose that we have chosen specific sparsification scheme Ω C li that satisfies Assumption 3.2, and then there must exist a sequence of sparsification schemes {Sp t E } that starts with 22 = 52 and ends with Qr = Ω satisfying P4+1 C Q4 E [T - 1]. and Q2 and $2 only differ by one element. Since Assumption 3.2 holds for Ω = Ωr Lemma B.4 implies that it also holds for all spar- sification schemes in the sequence, as Q2 S2 Vt [T]. Then by applying Theorem 3.3 repeatedly to Q4 and $2+11 for t = 1,..,, T - 1, we obtain the following corollary.
1.0
100
200
300
400
500
600
700
800
Sparsification Step
Figure 1. Plot of minimum and maximum cigenvalues of the spar- sified Hessian matrix during the sparsification process.
3.4. Summary
Corollary 3.4. For any sparsification scheme Ω C Ω satis- fying Assumption 3.2, the corresponding sparsified Hessian matrix Hq has the following properties:
We summarize this section by making few remarks on the significance of Theorem 3.3 and Corollary 3.4:
Amax HHa) S Amax (
Since the true Hessian matrix is positive definite, Corol- lary 3.4 implies that any valid sparsification scheme also maintains positive definiteness. This crucial for computing the search directions, as the approximate Hessian also needs to be inverted.
Amiu (Ha) 2 Amiu (H),
where 11 = Hgg = (x). The equalities hold if and only if Ω = Ω2
2.
The theorems indicate that the sparsified Hessian has a smaller condition number, which guarantees better nu- merical stability in solving linear systems. Even better, a smaller condition number makes iterative sparse lin- ear solvers (e.g. conjugate gradient method) converge faster.
3.3. Numerical Verification
We numerically verify Theorem 3.3 and Corollary 3.4 through a simple experiment that illustrates how eigen- values evolve during the sparsification process. Specifi- cally, we generate a random matrix T E TR with entries Try did Unif(O, 1), and construct the corresponding Hessian matrix H. We then iteratively set one nonzero element of T to zero each step, recalculating the minimum and maxi- mum eigenvalues of H at each iteration, until only the first row and first column of I' remain nonzero (see Theorem 4.1 for the rationale of this setting). We repeat this proce- dure five times, and illustrate the change of eigenvalues and condition number of the sparsified Hessian in Figure 1.
The controlled condition number is essential for the theoretical analysis of quasi-Newton methods; see for example Section 3.2 of Nocedal & Wright (2006).
4. Corollary 3.4 is valid for almost any sparsification scheme. This theoretical guarantee allows for highly flexible algorithmic designs, which greatly generalizes prior works on this direction.

---

The Chaare-Pllus-Low-Rae Quasi-Newton Method
5. The only requirement is Assumption 3.2, which is very weak: it can be satisfied with extremely sparse matrices. In Section 4.2, we provide a constructive method to let sparsification schemes easily satisfy this assumption.
Algorithm I The Sparse-plus-lo--rank quasi-Newton
method for entropic-regularized OT
Input: Initial point vo. maximum density Pinx E [0. 1,.
maximum shift Tmax > 0, stopping criterion Etol > 0
Output: Ik
4. The Spars--Plu-Low-Raall Method
1: Set Puain = 0.01 Pmax Po = 0.11 Pmax. Go = 52" (po)
2: Compute fo = f( 9o = g Ho = H
4.1. Overview
3: Sparsify Ho according to scheme Ro to obtain Hig
Despite the various benefits as explained in Section 3, the
4: Set To = min m [Tmax: ||
Hessian sparsification method has one intrinsic limitation:
Compute dg - (41 + -ood 'o0
it highly relies on the sparsity property of the true Hessian
6: Select the step size an and update 21 = zu ando
matrix. When the true Hessian is dense, no sparsified ver- sion would perform well, and we need better models to approximate the Hessian matrix. To this end, we introduce the paaee--plu-loow-rnn (SPLR) approach that improves the approximation by adding low-rank terms to the sparsi-
7: for k: = 1, 2, do
8: Compute fa = f(x)), gr = g(xx), H2 = (ax,)
9: if || < Exal then
10: return Zk
11: end if
fied Hessian. At a high level, we adopt the quasi-Newton framework solve (3) based the update rule
12:
Update pa according (8)
13:
Compute 24 52" (px) with Algorithm 2. and spar- sify H2 to obtain H
Ik+1 = Ik- By 9k+
14:
Compute Sk-1
11-1. Uk-1 = 9k- 9k-1 Compute a, b, v according to (7)
where B& is an approximation to the true Hessian matrix Ha. and ag > 0 is suitable step size that satisfies the Wolfe conditions given in (9). In our method, Bg consists of three parts:
15:
Let L (aauT + buot if yf > 10 1|yl2
16:
10,
otherwise
17:
Update T min T [|aall)
B₂ = Hq + (avn" T + boe")) + TI,
(5)
18:
Compute dy = - BB 9k. where By H +
19:
Select the step size Og and update xx+1 = I& + tayd
where Hq is the sparsified Hessian matrix according to some sparsification scheme Qt auu" + bve" is a rank-two approximation term, and > O is shift parameter. Note that all these terms may vary with the iteration number k.
20:
end for
Theorem 4.1. Define S2" - -i 1,
1}}, and let Har be the corresponding sparsi- fied Hessian matrix. Then Oo satisfies Assumption 3.2 with (Ha-) )d > 0.
[m].
The intuition behind (5) is that when H3 is truly close to a sparse matrix, Hn would be able to capture most of its information. And when this is not the case, the low-rank term can then compensate for the possibly dense difference Hgg The shift term 7/ is used to enhance the numerical stability when inverting the approximate Hessian matrix. Overall, BN is expected to perform as well as existing sparse Newton methods, and to show its advantage when H₃ is relatively dense.
Remark 4.2. The choice "i = 1 or j- 13" in the definition of Ω is not essential. One can use "i = i or j j" for any fixed values of i" and j Ω basically means keeping one row and one column of T and setting other elements to zero
Theorem 4.1 shows that it is reasonable to use Ω as our sparsification scheme. However, Ω may contain too few elements to provide good approximation to H. For better performance, we keep a fixed proportion P of the largest elements of 7, leading to a sparsification scheme M(p), where p E [0, 1] represents the density. To satisfy Assump- tion 3.2, we take the union of Ω* and Mp)), denoted by 2 (p) == ((p) U Qt. Then we can show that D2 (p) also satisfies Assumption 3.2 according to Lemma B.4.
In Algorithm 1, we first present our main SPLR algorithm to solve entropic-regularizedd OT, and then elaborate its de- tails in subsequent sections, such as the choice of Ω the specification of (a, b, u. v). etc.
4.2. Sparsification with density
One of the most important ingredients of the proposed SPLR algorithm is the choice of the sparsification scheme f2 in (5). First recall that an important condition for Corollary 3.4 to hold is that Q satisfies Assumption 3.2. Below we first show that there is a "minimal" scheme Ω² meeting this requirement, and then any scheme containing Ω² would also satisfy Assumption 3.2.
Formally, first define an operator se selectt.n large p). which takes an n x (m 1) matrix T and a density as in- puts, and outputs a set Ω consisting of the coordinates of the largest [pn(m - 1)J elements in T. Then the sparsification scheme with given density is obtained via Algorithm 2.
5

---

The Pharse-Plus-Loow-ankk Quasi-Newton Method
Algorithm 2 Sparsification scheme with a given density
4.4. Practical implementation
Input: Dual variable vector I (a" yT, density pa-
Shift parameter The shift parameter The in Algorithm I is not necessary for theoretical analysis, and in fact, one can globally set т = 0 without breaking the algorithm. However, in practical implementation, it brings various ben- efits, for example, stabilizing the linear system B Sk. and potentially accelerating iterative linear solvers such as the conjugate gradient method. This is because 1& + r has a smaller condition number than H
rameter E [0, 1]
Output: Sparsification scheme (p)
1: Compute T' = T(a, 3)
2: Compute (p)) = select
3: Set Q2" (p) = (pp) US
4.3. Low-rank terms
To avoid introducing large approximation error, we dy- namically set The to be the current gradient norm I19x]l, so that the The / term in (5) is negligible when is close to the optimum. An additional safeguard is to set a maximum shift Tiaxx in case |I9x|| is too large at the beginning. So overall, we take Th = min Fmaxx [| in each iteration.
To enhance the approximation quality when the true Hessian not well captured by sparsification alone, we incorporate low-rank correction term. Specifically, suppose that we are at the (A + 1)-th iteration of the Newton-type optimization procedure. We approximate H&+1 by a matrix B&+1 of the form:
H8+1 RJ Bx+1 U1ttt + +auu + bve + Tk+1d,
Adaptive density selection Thanks to the theoretical guar- antee presented in Corollary 3.4, we have a high degree of freedom to design the sparsification scheme in each itera- tion. In our implementation, the density parameter Pk varies according to [9x||- If ||9H|] decreases compared to the pre- vious iteration, it means that the previous B potentially provides a good approximation to H so we can try a more sparse H& in the current iteration, thus accelerating the search direction computation. Otherwise, we should increase the density to obtain more precise approximation to Hy. Based this idea, the update rule for Pk is:
where Hh+1 is the sparsified version of Hk+1 according to a scheme Ω, and a, b E I, u E Rn+m-1 are to specified later. Motivated by various quasi-Newton methods, especially the update rule, we determine a, b, u, v by the secant equation, namely the first-order approximation of y at 1k+11
(BFGS)
g(z) 9k+1 + Hk+1 (a zx+1) + O(1 xk+1ll²).
Replace with 2g. and we get
= 0*++ + (xF - Ta+1) + 011 Fx+1 119)
maxx 0.99p if || < 199--11 (8) min (Pmax: 1.1p
Pk
otherwise
Let s = Æk+1 Σy and yk 9k+1 - Ik- By ignoring the remainder term, we have y/& Rs H4+18 So it is reasonable to require the approximation of Hx+1, namely B to satisfy:
Line search method Finally, we use the
line search algorithm (Moré & Thuente, 1994) to compute the step size ag > 0 that satisfies the Wolfe conditions:
y/s = BB418h Sk = -H +1 + amns + boo + Th+1 D)kk- (6)
f(ax + ands) S f(xx) + (cian(oEda)) dx), [g(T + [Tdd))ddd de 2 (2 dx),
One solution to (6) is:
(9)
(H + Th+ T)sk:
yiks
1
where < < 1/2 and c 02 pre-specified constants.
a =
Vk 8k
b b=
of (u1 TK+1
(7)
Although such a modification makes the approximation dense again, its inverse can be computed very conveniently with sparse or vector-based arithmetics (see for example Section 6.1 of Nocedal & Wright, 2006):
5. Convergence Analysis
In this section, we prove that the proposed Algorithm 1 en- joys global convergence property, and the convergence rate is at least linear. Much of the theory has been developed in classical quasi-Newton literature such Byrd et al. (1987), but the key challenge here is to verify certain properties of the approximate Hessian matrix B
B TT (H
k+ + Tx+1l)- + ELSAST.
where EA- = 1/ 2x) and U = T- Since (Hb+1+ T is still a sparse and positive definite matrix, linear systems associated with it can be efficiently solved via either direct methods such as the sparse Cholesky decomposition, or iterative methods such as the conjugate gradient method.
In particular, the key to the convergence of quasi-Newtur methods is the condition number of Bx+1 in each itera- tion. Therefore, it crucial to bound both the smallest and
6

---

The Pharse-Plus-Loww-Rndd Quasi-Newton Method
largest eigenvalues of B8+1 for every k. We then have the following key findings.
6.1. Synthetic data
We first consider two synthetic datasets that have been ana- lyzed by existing sparse Newton methods (Tang et al., 2024; Tang & Qiu, 2024):
Theorem 5.1. Assume that there is closed set D such that
Amin (Hin(H(a))) (H(z)) 2 L, Amax S
Synthetic I: M = (M,j) has uniformly distributed entries, sid ~ (ni 1), and a = n 1,. b = m mm-1mm i.e., Mij
for some constants L., > 0 and all E D, and that (1 t) + txæ+1 D for all E [0,1]. 1]. Then
Synthetic II: This setting approximates OT between two continuous distributions: an exponential distribution with mean one, and normal mixture distribution 0.2-NV (1, 0.2)+ 0.8 . (3, 0.5). The vectors and b contain the discretized density function values of the two distributions, computed in the following way: let z; 5(i-1)/(n 1), E [ni]. and yj = 5 1)/(m - 1).j E [m]] be equally spaced points on [0, 5]. and let f1 and fz be the density functions of the two distributions, respectively. Then set a, = -fi fz;), by = f2())), a, a. (E kg). and by /// k= bk). The cost matrix is set to Mij- = (ar; )
Amin (Bx+1) 1 2 (2 + -
Xmax 1Bx+1)= S 2U + Tmax-
After bounding the eigenvalues of BE+1, we then obtain the global convergence of the proposed algorithm.
Corollary 5.2. Let To be an arbitrary initial value, and { be generated by Algorithm 1. Then
lim (1 = O. 10+000
Finally, we show that Algorithm 1 at least has a linear con- vergence rate.
We simulate different scales of the problems, 1000, 5000, and 10000, and do not run Newton or APDAGD for n 2 5000, as they are too time-consuming. The results are given in Figures 2 and 3.
Theorem 5.3. Let| f be the optimal value of I (x). Then for all k 22 there is constant 0 < r such that
As can be observed from the plots, the existing sparse New- ton method SSNS performs well in Synthetic II but is slow in Synthetic I in terms of the run time. This is because the Hessian matrix is relatively dense in I, and SSNS needs to keep large number of non-zero elements in the approxi- mate Hessian, leading to slow linear system solving. SPLR, on the other hand, performs well in both settings, showing its adaptivity to different OT settings.
f(æx+1) - S (r(f(z)) - f*.
Although the theory only gives linear convergence rate, our numerical experiments in Section 6 suggest that in many cases, the proposed SPLR algorithm achieves super- linear- like convergence speed.
6. Numerical Experiments
6.2. OT between a pair of vectorized images
In this section, we evaluate the performance of the proposed SPLR algorithm via series of numerical experiments, and compare SPLR with a number of widely-used algorithms for solving cetropic-regulariied OT: 1. the Sinkhorn algorithm (equivalent to block coordinate descent, BCD); 2. the adap- tive primal-dual accelerated gradient descent (APDAGD, Dvurechensky et al., 2018); 3. L-BFGS; 4. the Newton method: 5. the SSNS algorithm (Tang & Qiu, 2024).
We then study the OT problem between a pair of im- ages. Specifically, we randomly select two images from the MNIST (Lecun et al., 1998) or Fashion-MNIST (Xiao et al., 2017) dataset, and let the and b vectors be their flat- tened and normalized pixel values. The cost matrix holds the & fi-distances between individual pixels. These problems have a size of n = m = 784, with results shown in Figure 4 (for MNIST data) and Figure 5 (for Rishion-MNIST data), respectively.
We consider both synthetic and realistic OT settings, and fix the regularization parameter to be 1 0.001. To make 7) comparable for different problems, we normalize all cost matrices to have unit infinity norm, i.e., M <
The pattern in Figure and Figure 5 shows that the Sinkhorn algorithm (i.e., BCD) and the first-order method APDAGD have quite slow convergence speed, and the sparse New- ton method SSNS significantly accelerates the optimization. With this foundation, SPLR performs even better, as it com- bines the advantages of both sparse Newton methods and low-rank quasi- Newto- methods.
We use the gradient norm [[g()))] to measure the optimiza- tion error of the current iterate Zg. Additional test examples with different cost matrix settings and values are given in Appendix A.2 and Appendix A.3, respectively. The experi- ments in this section can be reproduced using the code our Github repository
httpss: htt///gtttuttt come/o//////oeeeeee Aoblex/ numeri ca 1 -exxeeereeeee
7

---

The Tparse-Plus-LLovy Quasi-Newton Method
Synthelic (nn = m= 10 0.001)
Synthetic (0== m= 5000, n000.n=0000 0.001)
Synthetic I In m= m 00000..=0.001
PCD
BCD
523
AFDAG
Duil
LaroS-Dua
Lop30 Gradient Norm
Log Crradient Norm
-3
Log10 Gradient Norm
SR.R
6IN1
5PU1
-4
-5
4
-5
20
40
60
80
100
6
10
3
6
Iteration Number
Iteration Number
Iteratoon Number
Synthetic (i1 - m- 1000, o 0.001)
Synthetic ((- in - 5000, n- 0.001)
Synthetic 10000. 0.001)
PCD
BCD Laros-oe
AFDAGO
LEEG Cual
Locto Gradient Norm
Logio Crndient Norm
U
Logi Gradient Norm
55NS
SPLA
-3
5TU1
5SNS
4
-4
-5
-6
3 Run Umelseconnds
10
12
14
10
15
20
25
30
Run time(seconds)
Run timetscconds)
Figure 2. Performance of different algorithms on synthetic data I. Top: Gradient norm vs. iteration number for different problem sizes. Bottom: Gradient norm vs. run time.
Synthetic 1000, 1y- 0.001)
Synthetic (nn 5000, n 0.001)
Synthatic 10000, p- 0.001)
18865-0 6065
(8905
Logio Gradient Nomn
18F00 Dual
Log10 Gradinnt Norm
Logi0 Grodient Norm
5065
Mortan
DSNS
3
6
200 Number
500
Iteration Number
Iteration Number
Synthetic
1000, g-
Synthetic (0n n- 0.001)
Synthetic 10000. g- 0.001)
0
18805-0.
10005
Logio Gradient Nomm
Wortan
(8F05 Dual
Leg10 Gradient Nomm
SP
5065
Log10 Gradient Norm
2
5P
5965
1.
-
b
4
6
5.
75
100
100
300 time/socandss)
700
timelseconds)
Figure 3. Performance of diffeerntt algorithms on synthetic data 11. Top: Gradient norm vs. iteration number for different problem sizes.
Bottom: Gradient norm vs. run time.
MNIST (ID1-2, ID2-54598, L.-orrm n- 0.0011
MNVIS (239, ID2-43062 Syrtorm, n- 0.0011
MNIST [ID1 1701-170 17300, ID2-49947, (y-torm n- 0.001)
500
ECD
BCD
AFDACO
AFDAGD
ACC
Logi0 Gradient Norm
2. 1.
Logi0 Gradient Norm
Log30 Gradient Norm
5516
5SN5
-2
3.
&
-4
1
4.
0.5
2.0
2.5
1.0
1.5
3.0
0.0
1.5
2.0
Run limelseconds)
Run timelserondds)
Run (imelsecoods)
Figure 4, Performance of different algorithms on the MNIST data.
8

---

The Phaarsc-luus-Loweee Quasi-Newton Method
FashionMNIST (D1-2, ID2-54698, Junorm, n- 0.001)
FashionMNIST ID1-239, ID2-43981. d-nomm. 0.002)
FashionMM (ID1-17390, ID2-49947, 4yttorm y 0.001)
APMMC
APDADD
Loçto Locto Gradient Norm
uero Dud
Logio Cradient Norm
L 1
usres Dut
Log310 Gratient form
1
4.
0.0
0.5
2.5
1.0
1.5 Run timesecoodds
3,0
0.0
2.5
3.0
timetseconds)
Run timelseconds)
Figure 5. Performance of different algorithms on the Fashion-MNIST data.
6.3. OT between two classes of images
ImageNet (Tench Cassette player. fy noror,, n= 0.01)
BCD
Finally, we reproduce the ImageNet experiment in Tang & Qiu (2024) that uses OT to characterize the difference between two data distributions. In particular, one class of images in the ImageNet dataset (Deng et al., 2009) is
APTACO
Log10 Gradient Norn
L CS-Dual
Newton
SSNS
treated as one distribution, and each image in this class is
an observation.
This experiment is interesting in that different values of the regularization parameter have a large impact the performance of algorithms, as shown in Figure 6. With 0.001, BCD and APDAGD again
-7
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
Run time (secann))
ImageNet (Tench vs Cassette player, Z--norm, ! 2 0. 0.001)
small regularization, have slow convergence speeds, and L-BFGS and SSNS perform quite well. The classic Newton method causes numerical issues,
APDAGD
Logi0 Gradlent Norm
-2
Newton SSNS
is not shown in the plot of f 0.001. When 7 is increased to 0.01, all methods converge faster except for SSNS. This is because a larger typically leads to more dense Hessian matrix, so SSNS cannot use sparse matrix to approximate the Hessian well. In contrast, SPLR deals with this situation well via its low-rank term, making SPLR perform consistently well on different v values.
SPLR
-4
-8
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
time (sccondds
7. Conclusion
Figure 6. Performance of different algorithms the ImageNet data. Top: 17 = 0.01. Bottom: 17 = 0.001.
In this paper, we propose the SPLR quasi-Newton method to solve large-scale entropic-regularized OT, as further extension of existing sparse Newton methods including the SNS (Tang et al., 2024) and SSNS (Tang & Qiu, 2024) meth- ods. The design of the SPLR algorithm is highly dependent on the deepened theoretical understanding of the Hessian sparsification technique, which may be of interest by itself. On the other hand, the low-rank term introduced in SPLR effectively overcomes the limitation of sparse Newton meth- ods in handling dense transport plans. In this sense, SPLR combines the best parts of purely low-rank-bart based methods methods
result for the Sinkhorn algorithm is given in Carlier (2022), which shows that the Sinkhorn algorithm has a linear con- vergence rate that only depends on /o/// and not the dimension (n, m). SPLR and other related solvers also have such properties.
of interest to understand whether
(e.g., L-BFGS) and purely sparsification-base (e.g.. SNS and SSNS). We anticipate that the technique developed in this paper would boost future exploration of highly efficient solvers for OT.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
One potential future research direction is to study dimension- independent convergence rates of OT solvers. One known

---

The Phaass--lussLoowRank. Quasi-Newton Method
Acknowledgements
gradient descent better than by Sinkhorn's is algorithm. In Proceedings of the 35th International Conference on Machine Learning, pp. 1367-1376, 2018.
Yixuan Qiu's work was supported in part by National Natu- ral Science Foundation of China (12101389), Shanghai Pu- jiang Program (21PJC056), MOE Project of Key Research Institute of Humanities and Social Sciences (22JJD1 10001), and Shanghai Research Center for Data Science and Deci- sion Technology.
Genevay, A., Peyre, G., and Cuturi, M. Learning generative models with Sinkhorn divergences. In Proceedings of the Twenty- Firenty-First International Conference on Artificial Intelligence and Statistics, pp. 1608-1617, 2018.
Guminov, S., Dvurechensky, P., Tupitsa, N., and Gasnikov, A. On combination of alternating minimization and Nes- terov's momentum. In Proceedings of the 38th Interna- tional Conference on Machine Learning, pp. 3886-3898 2021.
References
Altschuler, J., Wiles-Weed J., and Rigollet, P. Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration. In Advances in Neural Information Processing Systems, volume 30, 2017.
Huynh, V., Phung, D., and Zhao, H. Optimal transport for deep generative models: State of the art and research challenges. In Twenty-Ninth International Joint Confer- ence on Artificial Intelligence, volume 5, pp. 4450 4457, 2021.
Arjovsky, M., Chintala, S., and Bottou, I.. Wasserstein generative adversarial networks. In Proceedings of the 34th International Conference on Machine Learning, pp. 214 223,22222222222222 2017.
Brauer, C., Clason, C., Lorenz, D., and Wirth, B. A Sinkhorn-Newton method for entropic optimal transport. arXiv preprint 0rXiv:1710.066:: 2017.
Laclau, C., Redko, I., Matei, B., Bennani, Y., and Brault, V. Co-clustering through optimal transport. In Proceed- ings of the 34th International Conference on Machine Learning, pp. 1964, 2017.
Byrd, R. H., Nocedal, J., and Yuan, Y-X. Global con- vergence of cass of quasi-Newton methods convex problems. SIAM Journal on Numerical Analysis, 24(5): 1171-1190, 1987.
Lecun, Y., Bottou, L.,, Bengio, Y., and Haffner, Gradient- based learning applied document recognition. Proceed- ings of the IEFE, 86(1 (1:111)2227 11:2278-2322 1998.
Carlier, G. On the linear convergence of the multimarginal sinkhorn algorithm. SIAM Journal on Optimization, 32 (2):786-794, 2022.
Li, M. Yu, J., Li. T., and Meng, C. Importance sparsification for sinkhorn algorithm. Journal of Machine Learning Research, 24(247):1-44 2023.
Chen, L., Gan, Z., Cheng, Y., Li, L., Carin, L., and Liu, Graph optimal transport for cross-domain alignment. In Proceedings of the 37th International Conference on Machine Learning, pp. 1542-1553, 2020.
Lin, T., Ho, N., and Jordan, M. I. On the efficiency of en- tropic regularized algorithms for optimal transport. Jour- nal of Machine Learning Research, 23(137):1 42, 2022.
Liu, D. C. and Nocedal, J. On the limited memory BFGS method for large scale optimization. Mathematical Pro- gramming, 45(1):503-522, 1989.
Courty, N., Flamary, R., Tuia, D., and Rakotomamonjy. A. Optimal transport for domain adaptation. IEEE Transac- tions on Pottern Analysis and Machine Intelligence, (9):1853-1865, 2017.
Magnus, J. R. On differentiating eigenvalues and eigenvec-
Cuturi, M. Sinkhorn distances:
Lightspeed computation of optimal transport. In Advances in Neural Information Processing Systems, volume 26, 2013.
tors. Economelric theory, 1(2):199-191, 1985.
Meyer Jr, C. and Stadelmaier, M. Singular M- --matrices and inverse positivity. Linear Algebra and its Applications, 22:139-156, 1978.
Cuturi, M. and Peyré, G. Semidual regularized optimal transport. SIAM Review, 60(4):941 2018
Montesuma, F., Mboula, F. M. N., and Souloumiac, A. Recent advances in optimal transport for machine learn- ing. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2024.
Deng, J., Dong, W., Socher, R., Li, I., Li, K., and Li, F.FF. ImageNet: A large-scale hierarchical image database. In 2009 IEEE Conference on Computer Vision and Pattern Recognition, pp. 248-255. 2009.
Moré, J. J. and Thuente, J. Line search algorithms with guaranteed sufficient decrease. ACM Transactions on Mathematical Software (TOMS), 20(3):286-377, 1994.
Dvurechensky, P., Gasnikov, A., and Kroshnin, A. Com- putational optimal transport: Complexity by accelerated
10

---

The Sparse-Plus-Low Piae--llu--Lo-Raatt Quasi-Newton Method
Nocedal, and Wright, S. J. Numerical optimization. Springer, 2006.
Peyré, G., Cuturi, M., et al. Computational optimal trans- port: With applications to data science. Foundations and Trends@ in Machine Learning, 116--6)33556677, 2019.
Sinkhorn, R. A relationship between arbitrary positive ma- trices and doubly stochastic matrices. The Annals of Mathematical Statisties, 35(2):876-87,, 1964.
Tang. X., Shavlovsky, M., Rahmanian, H., Tardini, E., Thekumparampil, K. K., Xiao, T., and Ying, L. Acceler- ating Sinkhorn algorithm with sparse Newton iterations. In The Twelfth International Conference on Learning Representations, 2024.
Tang. Z. and Qiu, Y. Safe and sparse Newton method for entropic-regularized optimal transport. In Advances in Neural Information Processing Systems, volume 38, 2024.
Torres, L. C., Pereira, L. M., and Amini, M. H. A survey on optimal transport for machine learning: Theory and applications. arXiv preprint arXiv:2106.0199 0119963,, 2021.
Villani, C. et al. Optimal transport: old and new. Springer, 2009.
Xiao, H., Rasul, K., and Vollgraf, R. Fashion-MNIST: novel image dataset for benchmarking machine learning algorithms. arXiv preprint arXiv: 17X8.0/708.077: 2017.
Yule, G. U. On the methods of measuring association be- tween two attributes. Journal of the Royal Statistical Society, 75(6):579-644 1912.

---

The Llwwe---lu-- Quasi-Newton Method
A. Additional Experiment Details
A.1. Effect of low-rank terms
To verify that the low-rank term in the SPLR algorithm indeed plays a significant role in improving the sparse Newton method, we conduct an ablation study that compares SPLR with a purely sparsification-base method. We consider a special sparse Newton method that has the same sparsification scheme and hyperparameter setting as SPLR, and its only difference from SPLR is the lack of the low-rank term in the approximate Hessian. We examine the performance of the two methods on the synthetic data introduced in Section 6.1, with the results shown in Figures 7 and 8.
It is clear from the plots that in all the settings, SPLR converges faster than its counterpart, which highlights the benefits of the low-rank term.
Synthetic| (11 - m- 1000, .- 0.001)
Synthetic [n- m - 5000, n- 0.003)
Synthetic I in m- 10000,000000 r0000 0.001)
Spane Wodton
Sperss Now.or
Spuas Ntwton
51.2
SPLR
Logt0 Log Gradient Norm
Log Craditen Norm
-3
LogI0 Gradiint Norm
4
-4
-5
-5
5
-6
300
200
300
400
500
2
6
10
2
3
5
6
Iteration Number
Iteration Number
Iteration Number
Synthetic| (11 11-m m- 2000, y- 0.001)
Synthetic ((- m - 5000, np-
0.001)
Synthetic In 10000. 0.001)
Spane Wodoon
Speras Now.or
Spatsc Nonion SPLk
$P.2
Logi0 Log10 Gradient Nomn
Log 10010Gon Gradient Nom
Log10 Gradient Norm
-4
-5
-5
16
0
4
0.0
0.5
1.0
1.5
2.0
2,5
3.0
3.5
6
1.0
12
Run tmetseconds)
Run time(secands]
Run timesseconds)
Figure 7. Comparing SPLR and sparse Newton method synthetic data I. Top: Gradient norm iteration number for different problem sizes. Bottom: Gradient norm vs. run time.
Synthetic (n- m- 1000,00 . 0.001)
Synthetic BI (n - 5000, n 0.001)
Synthetic BI In mm- 10000. -- 0. 0.001)
Spame Mextton
Spame Nowcer
Sparn Nowton
Log10 Gradient Nom
Log10 Cradiont Norm
Log10 Gradient Norm
-2
-3
14
-4
6
300
200
400
100
200
300
400
500
100
200
300
400
500
Iteration Number
Iteration Number
Iteration Number
Synthetic (nx- m- 3000, . 0.001)
Synthetic (. m- 5000, n- 0.001)
Synthetic II In -m 10000. gm 0. 0.001)
Spame Mexttnn
Spame Neweer
1.
Spane Vevton
Logi0 Gradient Nom
Log10 Gradient Nomn
Log10 Gradient Norm
-2
-2
-3
-4
-4
0.00
0.25
0.50
0.75
1.00
1.25
1.50
1.75
10
20
40
50
60
70
50
100
150
200
250
Run timefseconds)
Run time(secoodss)
Run time]seccodds)
Figure& 8, Comparing SPLR and sparse Newton method on synthetic data II.
12

---

The Sparse-Plus-Low-RRae Rank Quasi-Newton Method
A.2. Effect of cost matrices
In this section, we include additional examples in which the cost matrices are formed by the Euclidean distance instead of the & distance. Figures 9 and 10 show the results on the (rishion->MNNIS data and ImageNet data, respectively. Clearly, in all settings, SPLR is among the fastest solvers.
MNIST (ID1-2, ID2-54698, Izinorm, q- 0.0011
MNVIST (D1-239. ID2-43961, fynorm, - 0.001]
MNIST IID1-17990 ID2-49947, fynorm, o) 0.001)
200
ECD
BCD
Lngio Gradiet Nomn
Log30 Gradent Norm
Neadee
Gradion
-2
55106
5516
12
50N5
&
3
0 1
1
0.0
0.5
1.0
1.5
2,0
2.5
3.0
0.0
0.5
1.0
1.5
2.0
2.5
3,0
0.0
0.5
3.0
1.5
2,0
2.5
3.0
Run timesseconds)
Run (me(seoonds)
Run timelssccons))
TashionNNIST (D1-2, (D2-54098, nonorm n= 0.001)
FashionMNIST IDI-239, ID2-43961. & nomm. 0.0011
Fas (ID1-17300, 1D2449947,
0.001)
KD
Logio Gradient Nom
Legio Coadient Norm
Log10 Gradient Norm
2016
1
-2
-3
-4
4
1.0
1.5
3,0
Run time(seconds)
Run time(secondds)
Run timelseccods))
Figure Performance of different algorithms on MNIST (top row) and Fashion-MNIST (bottom row) data using the En-norm to form the cost matrices.
ImageNet (Tench vs Church. & n-nomm... 0.001)
Imagetet et (Tench vs Garbage truck. n-norm, 0.001)
Imagelett (Tench vs Golf ball, 4 norm,, n 0.001)
BCB
ECO
BCD
AFDACO
Legio Gradient Nom
L
Legto Gradient Norm
Logi0 Gradient
3.
&
L
d
4.
3.5
3.0
4.0
1.0
2.0 timelseconds)
timelseconds)
Run timesseconds)
Figure 10. Performance different algorithms on the ImageNet data using the fg-norm form cost matrices.
A.3. Effect of regularization parameters
Finally, we further validate the performance of the SPLR algorithm under the regularization parameter setting r = 0.01. Figures 11, 12, and 13 show the experiment results on the synthetic data, (Fashion-)MNISS data, and ImageNet data, respectively.
In summary, all these experiments show similar patterns to those in Section 6, validating the desirable performance of SPLR.
A.4. Computing environment
All experiments in this article are conducted on a personal computer with Intel i9- 19-13900K CPU, 32 GB memory, and Ubuntu 25.04 operating system.
13

---

The SearrsePluu-LLo--Reee Rank Quasi-Newton Method
Synthetic| 1 (1 1000, 0.001
Synthetic In 5000, 0.02)
Synthetic (6n 10000. q-
BCD
AFDAGO
LSFGS-Ou
JancS Cual
Lo030 Gro.dlent Norm
L0T
Logi00 0 Croddent Norm
Log10 Graalent Norm
SSRG
SPIR
-5
-6
-7
0.0
0.2
0.4
0.6
0.8
1.0
1.2
0.0
0.2
0.4
0.6
0.8
1,0
1.2
1.4
Run tmelseconds)
Run time(seconds)
Run timelseconds)
Synthetic (n m-mm- 5500, a 0.013
Synthetic (p mn 10000. - 0.01]
Synthetic (n 1000, 0.01)
E00
AFDACE
Log Cradienn Norm
1
LSFGS-Dul
BC3
Logi0 Grodlent Nom
N
Lanos Dual
8:3
5215
Logio Gradient
-3
sPun
3
0.00
0.75
1.25
1.50
2.00
10
20
40
50
100
150
200
250
300
Run bmelsecondss
Run time/seconda)
Run timelzeconds)
Figure/ 11. Performance of different algorithms on synthetic data (top row) and synthetic data (bottom row) with 1) 0.01.
MMST (ID1-2, ID2-54638, N-nom, n- 0.01)
MNIST (ID--239 ID2-43961, I.-nom, 0.01)
MNIST 102007 (0244447 Aymoomm 0.011
BCD
BCD
BCD
Log10 Crassien Norm
Log30 cradiemt Norm
Logi0 Gradiant
SONG
b.
Un
0.00
0.25
0.50
0.75
1,00
1,25
1.50
1.75
2.00
0.00
0.25
0.50
0.75
1.00
1.25
1.50
1.75
2.00
0.00
0.25
0.50
0.75
1,00
1,25
1.50
1.75
2,00
Run tmelseconds)
Run (me(secones)
Run timelscconds)
FashionMNIST (ID1=2, ID2 152688 E-nom,, n 0.01)
FashionMNIST [ID1 -239, 102-43981. Z-rorm, qm 0.01)
Fastion/MNST 1731- I02 -09947 A-noom / 0.01)
E00
BCD
BCD
Loc30 Gradient Norm
Logi0 Gradirt Norm
Log30 crodent Norm
Nowten
Neaton
55M3
50W5
3
&.
b
0.00
0.25
0.50
0.75
1.00
1.25
1.50
1.75
2.00
0.00
0.25
0.50
0.75
1.00
1.25
1.50
1.75
2.00
0.00
0.25
0.50
0.75
1.00
1.25
1.50
1.75
2,00
Run dmelseconds)
Run Eimessccondd))
Run Eimelseconds))
Figure 12. Performance of different algorithms on the MNIST (top row) and Fashion-MNIST (bottom row) data with rj 0.01.
Imageltet (Tench vs Church, E-nomm 1 0.01)
ImageNet (Tench vs Garbage truck, (-norm, q- 0.01)
ImageNiet (Tench vs Golf ball, fy-norm, n- 0.01]
ECD
BCD
BCD
AFDAGO
APDAGO
Log10 Nomm
Log Cradient Nomn
Logi0 Gradient Norm
SSNS
5515
57U1
$7J1
-4
1,6
12
5
6
Run timalseconds)
Run time(seconds)
Run time//sccands)
Figure 13, Performance of different algorithms the ImageNet data with v 0.01.
14

---

The Sparse-Plus-Low-RRae Quasi-Newton Method
B. Proofs of Theorems
B.1. Technical Lemmas
Lemma B.1. Let M be a matrix of the form
D2 Re
R
M
D2
where D1 E Brrr and D2 E xR are diagonal matrices, and R E RR has strictly positive entries. Suppose that M is positive definite, and M-1 has the partition
Axxr BT
Brxs
M M- =
Caxa]
and then we have A > 0, C> 0, and B < 0, where the inequality signs apply elementwisely.
Proof. Since M is positive definite, all of its diagonal elements must be strictly positive, implying that both Dj and D2 are invertible. By the inversion formula of block matrices, we have
== (D1 - RD; RD) C= (D2 - Rd R R)-1
DD - R(D2
R DT = RC.
Clearly, J = Di - is the Schur complement of the block D2 of the matrix M. By the properties of the Schur complement, we know that J is positive definite. and hence it is nonsingular.
Moreover, since R has strictly positive entries and D2 has positive diagonal elements, we have that RD2 'R" has strictly positive entries. Therefore, J can be represented in the form J = 8/ L,, where I, > 0, 2 p L), and p(-) stands for the spectral radius. Clearly, J is irreducible, so by Theorem Aii) of Meyer Jr Stadelmaier (1978), A J-1 has strictly positive entries. The same argument can be used to prove that C >
Finally, recall that B DJ RC. Since R O, C O, and D, has positive diagonal elements, we conclude that B 0.
Lemma B.2. Let M be matrix of the form
Arxr - BT Csxs
- -Brx
M
where A> 0, B 0, C > O, and the inequality signs apply elementwisely. Then M has positive eigenvalue such that any other eigenvalue of M in absolute value is strictly smaller than r, and the eigenvector = (v) uz associated with where E R" and E R can be normalized such that 01 >> and U2 < O.
Proof. Define
I 0 Q =
-T
and then it is easy to show that Q-1 = Q, and
1r -
0
A
-B
4
0
A T
B
QMQ-'= = 0 11
M.
- B
C
0 I
Therefore, M and M are similar to each other, and hence they must share the same eigenvalues. Clearly, A is positive Perron-Frobeminn Frobenius theorem, it must have positive and simple eigenvalue such that any other eigenvalue of M in absolute value is strictly smaller than r. Moreover, M has an eigenvector = (Vf (oT, o7)"" such that on E R",
matrix, so by the Perron
22 R.a 0, and M rt.
15

---

The Sparse-Plus-Loww-aam Quasi-Newton Method
Now let Q and then
Mu Q- Q
Therefore, is an cigenvector of M. Note that
1
t = Q-Q-
0
02
022
So 01 01 0 and 12 -62 O, which implies the stated result.
B.2. Proof Theorem 4.1
Proof. According to the sparsification scheme, the first row and first column of Tex is strictly greater than zero:
(T ). > O. (Th- ).1 > 0.
Suppose that Ho. is of the form:
Anxn
B Bux (m
Ha-
((m--1)x
(m- 1) ()x(m-)))
so in Hn- we have:
A diag(A) 0, D diag(D)>0 B₁. > 0, B, >0 Ci. > 0,, C
To show H.- > 0, first condider elements of H2-
When n. S n. (Hg-)? (( (Hg-) 2 B,,1 >0. i. When i > n.j n. (Hr-) 13 ((Hn))i, (H ) 2 Ci-n,1 Cj-n.1 >
ii.. When i n. n + 1, (Hn- 313 ((Ha): ),, (HQ )z,) 2 C1, > 0. iv. When i n 1, S n, (Hn-) Y2 ((Ha)) )i, (Ha-) )}.) 2 Aj. 0. When i = 1.j > (Hl )/3 ((Ha)) )., (Hn-) )j-) > B B i-jn D j-n.j-n > 0.
vi. When i > n.j = 1, (H )/5 = ((Hg)) (Hg-) ) 2 B D Dinniin > 0.
So we can assume that H3. is of the form:
N/xnn
B' Bax(m-11)
H2g =
Cmm-1)xn
(m- (mm-1x((m-1)
where
A' 0, D' D' 00 By > 0, B > 0, C. >0 C2 0.
Then consider elements of H
When i n, S n. (H Dij = ((H))) ). (H)-);) 2 A',1 > O.
ii. When >n,jj n, (Ha-); )ig = ((H?-). ). (1)))))) 2 D' D1
1113--0 1.j-n > 0.
iii. When nn n, (Hit) = ( Ci, > 0.
iv. When n.j < n, (1 (HH.)) D'. > 0.
16

---

The Sparse-Plus-LowwRRand Quasi-Newton Method
Overall, have (Har) ) 0 for p 4, and hence Hg- satisfies Assumption 3.2.
Lemma B.3. Given a sparsification scheme Q that satisfies Assumption 3.2, the sparsified Hessian matrix Hy has the following properties:
1. The eigenvector ((Ha) corresponding to Xmax (Ha) can normalized to have strictly positive entries;
2. There exists neighborhood (HH)) of Hg such that Xmax (H), H E ((Haa is differentiable.
Proof. Since Hn satisfies Assumption 3.2, there exists integer l 0 such that H 0. According the Perron- Frobenius theorem, u( ()) can be normalized to have strictly positive entries and Amax(HHg is strictly greater than all other cigenvalues of Hţ. Since Hi and Hn have the same eigenvectors, we can conclude that u ((H)) = u (HH)) can be normalized to have strictly positive entries and the corresponding eigenvalue Amax Ho) is strictly greater than all other eigenvalues of He. This means that Amax (g)) a simple positive eigenvalue and therefore Amax differentiable at Hg according to Theorem
Magnus (1985).
Lemma B.4. Given two sparsification schemes Do and S2, suppose that 22 S Po. If satisfies Assumption 3.2, then Uo also satisfies Assumption 3.2.
O 1nxn Th
T Om - 1)x
Proof. Define Qg = Lo and A =
Then we have
Ha. = Ha, + A.
Since Q2 satisfies Assumption 3.2, there must exist an integer such that (Hg,) 0. Then for Hag, we have:
(Ha,)* 1 = (Ho, + +)* 2 (Hn,) )* + 2 (Ha,) )a 0,
where all inequality signs are elementwise. This means that Hn, also satisfies Assumption 3.2.
B.3. Proof of Theorem 3.3
Proof. First, we prove Amax (Ha)) Amax (Ha,).
Define = i, = n + f.B = (Hop) )pg. then the difference of Hessian is:
Ha. Hns = 3 (e + t = -BJ,
Define M(x) = Ho. - <J,[(x)= = Amax (aaxMM(M (x)); & E R, then we have:
Amax (HQ.) - Amax (Ho,) ((3) - 1(0),
Since Ho, satisfies Assumption 3.2, we can show that Amax is differentiable on (M)) E [0, []} according to Lemma B.3. Suppose the eigenvector associated with Amax MamxxMM (M))) is ug. then the derivative at M(s) is:
0 0xmax
uu
0M M=
thus ((()) = u,, -J). According to the Lagrange's mean value theorem, 3 E (0, 3) such that:
((3) - 1(0) (B()) 0) 0)
{n -JB8
=-3 (tr(uuuu + e
(10)
= =-3 [tr((u op ug) + tr(u c uz)
--28(uc),
17

---

The Pharse-Plus-LoowRRatt Quasi-Newton Method
According Lemma B.3, ug can be normalized to have strictly positive entries that (uz)p(Up) 0, which means 1(3) 1(0) Amax Ha, Amax (Hag) < 0. Then we prove Amin (H < Amin (Hn,).
According to Tang & Qiu (2024), we have M(s) > thus invertible, then according to Lemma B.1, (()) is of the form:
A B C
B
M(K) =
where A > 0,C 0, < 0. Then according to Lemma B.2, we can show that (max(M(()) is simple positive cigenvalue, the corresponding cigenvector Vx = =v}', v²'" where V1 E R",,, V2 E 1m-1 can be normalized such that Vi > 0, V2 0. Since V. also the eigenvector of (() corresponding A we can tell that (Vx)p(Vi))
where p S n, Q >
Define h(x) Amin (M (x)), then similar to (10) we have h(8) h(0) = -28((v)) where & E (0, 8). Recall that (Ve) < 0, so h(3) h(0) > O, which means Amia (Hn,) < Amin (.,,).
B.4. Proof of Theorem 5.1
By definition,
B&+1
Hatt + auns" T + bbv" + Tkt11,
where
yko (H2t Tk+1)sks a
Us 8k
of (H + Tk+1)sk
By the design the algorithm, we have uI Sk > O, so auu? and deu are rank-one matrices with a > and < Then we have
Amax (auu²) = au u, Amax (mox()
It is well known that Amaze (A + B) < Amax (A) + Amax (B) for symmetrie matrices A and B,
V yk y Ss
Amax ( (Be+1) Amax (H&+1)+
+
Tk+1-
Corollary 3.4 shows that Amax (Hh+)) I Amax(HH+in and then by the assumption of the theorem, we have
Amax (H Amax (Hx+1) U.
On the other hand, since f(x) is twice differentiable, the mean value theorem indicates that
ys g(xx+1) -(xx)) = G2 xn) GRSK.
where
G [ H(11 t)æy + t2x+1)dt.
Again by the assumption, (1 - f)xr + tx D for all t E [0, 1], so for any v E R=+m-1
sf GGGU== (:' of ( T) + T88+1)ed7 2 { Lu cdT Lof
and similarly, of G'G U v. This indicates that
0< L Auna (Gt) < Xmax )ma) S U.
As result,
yh yk
$ G
w Ggu
au T
1 5k
$ GLSE
T
w w
18

---

The Pharse-Plus-Low--Rndd Quasi-Newton Method
where (G²) 8//2 is well-defined since G2 is positive definite. Then by the properties of eigenvalues, get au u S (Gk) U. Overall, we have
Vk Uk Vi 8k
Amax (BR+I) S Amax (H (lttt +
+ Tk+1 S 2U + Tmax-
Now consider the inverse of Bx+1- Using the Sherman Moriannn--orosonn- Morrison-Wodduryy formula, we can obtain
B U* (H
g + +w+11)) +
1
vk vE 8k
where U - - (y -
- yusk Let H4441 1 + Th411, and then Hx+1 is positive definite H1 is positive definite. So for any vector E 1=+m-11
UU
Amax UU=) Amiu (Hx+1) -|10112
w - H Uv S Amax
Amii (H&+1)
Note that
yk yk (vk sx)²
UT = U (UkOR + tk +
Since U is positive semi-definite, we have Amax Amnx(UUU) = 110 where || Therefore,
represents the operator norm for matrices.
2 21 2||G&sH#f11
sw + skUiT)
5
2 sk
y Sk
$f Gk
G ||s
2||
$ sk
sf GaSSh
of GaSk
20
Similarly,
T yk
y yk
$ U
8888
T
(yl )2
V/ 8k
88 GLSE
I
So overall, ||TT|| & 1 3U/L. I,, and then
(U Amax(U
Amax Amin (H&+1)
1 (
30
Amax 10"H <
+
L
L
Finally, we obtain
3U
Tssk
1
3U
Amax B (Brt)) <
1
2
L
V 8k
L
Clearly, Bit positive semi-definite, so
1
3U -1
Amin (Bk+1)=
2
L.
Amax
Bx+1
L
B.5. Proof of Corollary 5.2
Consider the level set D (o: f(z) S f(zo)). Clearly, D closed convex set, and there exist constants L, U > 0 such that Amiu (H(x)) (x)) > L and Xinaxe (x)) U for all D.
Let 0g be the angle between Ik and the search direction Pk -B= 9k- Then clearly,
- Pk
P BuPE
Amin |
Amin (Bx)
cos cos 0g
(1I)
||
[BaPK]]
[1pa]]
Amaze BBx pk
Amax (Bx)
19

---

The Phpars-Pluus-Low-Rae Rank Quasi-Newton Method
Then by Theorem 5.1, have
Amin (Bk)
30
cos 0
2 =
2 2
0.
Amax (xx)
L
2U Tmax
Zoutendijk's 's theorem (see for example Theorem 3.2 of Nocedal & Wright, 2006) shows that
Σ coos 0z) 11)
00,
so we must have
2
Σ 11gx 1| < xc,
420
which implies that ||xx|| - 0 as At + 00.
B.6. Proof of Theorem 5.3
Similar to the proof of Corollary 5.2, consider the level set D = fx: f(x) S f(zo)}. and we have Amia (x)) 2 L and Amax H (x)) < for all r E D.
Lemma 2.1 of Byrd et al. (1987) shows that for any 2 1,
f(zx+1) - f [1 - Leje2 cos 0x] (f(xx) - f"),
where ट2 = (1 - e2)/U, and cos 0g is defined (11). Take
3U -1
L
r= 1 - c= 2
L
2U + Tmax
and then get the desired result.
20


Reviews:
[{"review_id": "n13VULPwx3", "reviewer": "Reviewer_2UKW", "summary": "The authors propose a quasi-Newton algorithm for solving entropic optimal transport (EOT) problems. The classical Sinkhorn algorithm enjoys linear convergence with a rate independent of the problem dimension but depends exponentially on the supremum norm of the cost function. There have been some recent works that aim to develop quasi-Newton methods for computing EOTs with super-linear convergence based on the idea of sparsely approximating the Hessian of the Kantorovich dual of EOT. The authors propose a flexible, theoretically justified sparsification scheme, and also combine the Hessian sparsification with low-rank Hessian approximation, similarly done as in L-BFGS. The authors deduce an asymptotic convergence guarantee and linear convergence rate for the proposed algorithm. The proposed algorithm seems to work very well in many synthetic and real-data experiments.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed review and thoughtful comments. We give our point-by-point responses below.\n\n---\n\n### Thm.3.3 and Cor.3.4\n\n1. We respectfully clarify that positive definiteness cannot be directly guaranteed via the Gershgorin circle theorem in this case. Specifically, for the sparsified Hessian $H_{\\Omega^*}$, the Gershgorin disc centered at $h_{n+1,n+1}$ has radius $\\sum_{j \\neq n+1} |h_{n+1,j}|$, admitting $\\lambda = 0$ as a possibility. Thus, the theorem alone does not preclude singularity. Moreover, you are correct that explicitly highlighting the improvement in the condition number would strengthen the impact of our results and we will emphasize this point more clearly.\n\n2. Thanks for this insightful suggestion. Indeed we have conducted numerical experiments on the condition number (this was actually done before we develop the theory), and we will add such experiments in the revised paper.\n\n3. There are two key reasons:\n\n   1. The assumption is analogous to the concept of *irreducibility* in Markov chains. Irreducibility is considered a mild condition because it only requires connectivity, not quantitative bounds on probabilities and it is *necessary but not sufficient* for stronger properties (e.g., ergodicity).  \n\n   2. Assum. 3.2 can be satisfied with *extremely sparse* matrices—_e.g._, retaining just one complete row and column (yielding a density lower bound of only $O(1/(nm))$).\n\n   We will add a brief remark in the revision to make this intuition clearer.\n\n### Thm. 5.1, Cor. 5.2, Thm. 5.3\n\nIndeed the current theoretical analysis still follows the classical framework of quasi-Newton methods, but we emphasize that one key ingredient of proving the linear convergence is to bound the condition number of the approximate Hessian matrix, and in our case it is highly related to Thm.3.3 and Cor.3.4, which are unique to OT problems. Your invaluable comments have encouraged us to actively explore novel technical tools to further optimize the linear rate constant.\n\nWe will revise the paper to acknowledge the significance of the theory in Carlier (2022).\n\n### Weakness 1\n\nSee our responses above.\n\n### Weakness 2\n\nYou are absolutely correct that the Sinkhorn algorithm has a linear convergence. Our phrasing was indeed misleading, and we will correct it in the revision. Our intended point was that, while the convergence is *theoretically* linear, the rate can *practically* resemble sub-linear behavior when the regularization parameter $\\eta$ is small ($1 - e^{-24\\\\|C\\\\|_{\\infty}/\\eta}$ approaches 1).\n\nAlso, Sinkhorn is indeed an alternating maximization procedure on the dual. We will revise for precision.\n\n### Weakness 3\n\nWe will update the manuscript accordingly.\n\n### Weakness 4\n\n$\\tilde{T}$ represents the first $(m-1)$ columns of matrix $T$, as defined in **Notation** (Line 68).\n\n### Weakness 5\n\nWhile thresholding the top x% of elements would indeed allow direct density control, we emphasize that existing Hessian sparsification methods do not employ this approach, nor do they provide theoretical guarantees for arbitrary density selection. Specifically:  \n\n1. SNS: Uses an elementwise threshold $\\rho$. Though we can set the Hessian to a desired density, there are no theoretical gurantees on invertibility or convergence rate.  \n\n2. SSNS: Sets a Hessian error threshold $\\delta_k = \\nu_0 \\\\|g_k\\\\|^{\\gamma}$ at each iteration, which depends on the gradient norm for theoretical guarantees. This scheme cannot be set to an arbitrary density a priori, as $\\delta_k$ varies with optimization progress.\n\nIn contrast, our proposed scheme explicitly enables any desired density level.\n\n### Weakness 6\n\nThe sparsification pattern $\\Omega^*$ consists of all coordinates in the first row and the first column.\nFor example, if $n=4$ and $m=4$, then $\\Omega^*=\\\\{(1, 1), (1, 2), (1, 3), (2, 1), (3, 1), (4, 1)\\\\}$.\n\nCould you please clarify in what specific way $\\Omega^*$ appears ill-defined? We would be happy to provide additional explanations.\n\n### Weakness 7\n\nWe sincerely apologize for this unintentional redundancy and we will carefully revise this section to eliminate repetition.  \n\n### Weakness 8\n\nThe discrepancy in our current plots likely stems from that each Sinkhorn \"iteration\" in our implementation includes both the $\\alpha$ and $\\beta$ optimization steps (two full matrix scaling operations).\n\n### Weakness 9\n\nWe thank you for this important technical observation. We will make the proof more rigorous in the revised manuscript. Below are the main ideas: we can first show that the optimal function value is finite, i.e., $f^*>-\\infty$. Then the set $D^*=\\\\{x:f(x)\\le f^*\\\\}=\\\\{x^*\\\\}$ is non-empty and bounded. Corollary 8.7.1 of [1] shows that $D_c=\\\\{x:f(x)\\le c\\\\}$ is bounded for every $c$, which implies that $D = \\\\{x: f(x) \\leq f(x_0)\\\\}$ is bounded.\n\n[1] R Tyrrell Rockafellar. Convex Analysis. Princeton University Press, 1970.", "claims_and_evidence": "Yes.", "methods_and_evaluation_criteria": "Yes.", "theoretical_claims": "**Thm. 3.3 and Cor. 3.4 on Hessian sparsification**\n\n1. The authors state that these results are significant because they guarantee positive definiteness of the Hessian after sparsification. But since they only sparsify off-diagonal entries and all Hessian entries are positive, diagonal dominance (or Gershgorin's circle theorem) trivially implies that the Hessian remains positive definite after any specification pattern. It does not, however, directly imply that the condition number improves. I think this latter point should be the main emphasis of these results. \n\n2. In the second and the third remark, the authors point out that their incremental sparsification decreases the Hessian's condition number, so it leads to numerical stability. While this is certainly true, there is no quantitative bound on the improvement of the condition number. If it is hard to justify theoretically, I recommend validating this claim by adding some experiments on plotting the condition number vs. iteration. \n\n3. It is not clear why Assumption 3.2 is \"very weak\" at this point. \n\n\n**Thm. 5.1, Cor. 5.2, Thm. 5.3 on the convergence guarantees**\n\n1. These results are derived from standard analysis of quasi-Newton methods, as the authors indicate in the appendix. Namely, the authors take the quotient space of the dual variables by setting $\\beta_{n}=0$ to get rid of the shift-invariance of the dual objective ($f(x)$ in their notation). While the dual objective is only convex and not strongly convex originally, after constraining on this subspace it becomes strictly convex. Further restricting this on the level set, one can obtain uniform bounds on the eigenvalues of the Hessian along the trajectory of the algorithm, provided that the objective values decrease monotonically. \n\nHowever, this standard analysis does not give linear convergence rate that is independent of the problem dimension. For instance, the rate $r$ in Thm. 5.3 may depend on $m$ and $n$, so the linear convergence stated there could effectively be sublinear for large-scale problems. The issue is that the condition number of the dual objective can be as large as order n. For instance, if at the current iterate the dual variables $\\alpha_k$ and $\\beta_k$ are the same (assuming $m=n$), then the $2n \\times 2n$ full Hessian has zero eigenvalue since the diagonal entries are $Cn$, which is the absolute row sums of the off-diagonal entries. Imposing $\\beta_{n}$ gets rid of one row from the (1x2) block matrix and the corresponding column of the (2 x 1) block, and the resulting off-diagonal sum is $C(n-1)$. The maximum eigenvalue is still of order n, but the minimum eigenvalue is of order 1. \n\nIn the state-of-the-art analysis of Sinkhorn (e.g., Carlier 2022 SIOPT), one addresses such issue by directly showing that the Sinkhorn iterates are uniformly bounded (by using closed-form expression of Sinkhorn iterates) and strong convexity of the entry-wise dual objective (in this case the exponential function) on a compact interval. Can the authors adapt a similar strategy and get a dimension-independent linear convergence rate of their algorithm?", "experimental_designs_or_analyses": "There are through synthetic and real-data experiments. Also, the additional experiment on the necessity of low-rank approximation provided in the appendix was very helpful.", "relation_to_broader_scientific_literature": "The method for computing EOT provided in this work could potentially be very useful for other researchers and practitioners.", "essential_references_not_discussed": "Guillaume Carlier, \"On the Linear Convergence of the Multimarginal Sinkhorn Algorithm\" SIAM Journal on OptimizationVol. 32, Iss. 2 (2022)10.1137/21M1410634", "other_strengths_and_weaknesses": "**strengths**\n\n1. The paper is relatively well-written and it was easy to follow. \n2. The proposed algorithm seems to be performing very well in various tasks, outperforming many existing methods. \n3. There are some theoretical guarantees on the Hessian sparsification and convergence guarantee. \n\n\n**weakness**\n\n1. Theoretical guarantees, especially for the ones for convergence, are weak. Please see the comments in the theory section. \n2. L043 \"Sinkhorn generality exhibits sub-linear convergence\": I don't think this is quire right. Do the authors have reference to back up this claim? Experimentally, in all plots Sinkhorn shows linear convergence. Theoretically, Sinkhorn enjoys linear convergence with a rate independent of the problem dimension but depends exponentially on the supremum norm of the cost function [Car 22]. Also, the authors write in the same paragraph that optimizing on the Kantorovich dual as a different approach than Sinkhorn, but Sinkhorn is just an alternating maximization on the dual and it looks like matrix scaling in the primal space. \n3. L125 and hereafter: Please just use $\\nabla f(x)$ instead of $g(x)$. \n4. L129: $\\tilde{T}$ is not defined. \n5. L163: The authors state that the existing Hessian sparsification by thresholding by values does not give a control on the density after thresholding. But one could easily threshold the top x % to get a desired density after thresholding, and this is essentially what the authors do. \n6. Thm. 4.1: the set $\\Omega^{*}$ is ill-defined. \n7. The first paragraph in Sec. 4.3 is a repetition from the previous page. \n8. Comparing the two plots in Fig. 1 in the second and the third column, it appears that BCD (Sinkhorn) is faster than SPLR in the iteration number. Since it is the cheapest per-iteration, it should be much faster than SPLR than in wall-clock time. The plots in the second row are misleading since SSNS is too slow and jams all other curves. \n9. L954: One needs to argue the sub-level set D is compact in order to obtain uniform bounds L,U on the eigenvalues of the Hessian over D. This is missing.", "other_comments_or_suggestions": "See above", "questions_for_authors": "See above", "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": "JPRLupLqeJ", "reviewer": "Reviewer_pXWg", "summary": "The paper concerns faster solver for the optimal transport (OT) problem, by proposing a new type of Hessian approximation within the quasi-Newton iterative solvers. Building on previous work that used sparse Hessian approximations, the authors introduce low-rank approximation added to the sparse format (hence sparse \"plus\" low rank in the title). The authors also provide convergence analysis, showing linear convergence.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed review and thoughtful comments. We give our point-by-point responses below.\n\n---\n\n### Weakness 1\n\nThanks for raising this question on motivation. Our approach is motivated by two key insights:\n\n1. **Problem-specific Sparsity**: The entropic-regularized optimal transport (OT) problem naturally exhibits (approximatly) sparse Hessian structure, particularly when $\\eta$ is small. This sparsity is unique to OT and related problems.\n2. **General Low-Rank Utility**: Low-rank approximations are successfully applied in quasi-Newton methods (e.g., BFGS, L-BFGS) because they efficiently capture dominant curvature information while maintaining computational tractability.\n   \nOur intuition is that the sparse-plus-low-rank (SPLR) structure may better preserve the problem's intrinsic geometry compared to pure sparse or pure low-rank approaches, leading to a faster convergence speed, as is shown in our empirical results. For example, if the Hessian \"residuals\" (i.e., the elements removed during sparsification) have some low-rank structure, then SPLR naturally captures this information.\n\n### Weakness 2\n\nWe thank you for highlighting this important consideration. Our primary limitation is that although superlinear-like convergence performance is observed, we currently lack formal theoretical justification for this observation. We will clarify such limitation in the revised manuscript.\n\n### Note 1\n\nWe thank you for this helpful suggestion to improve clarity. In the introduction, we will add a brief explanation of the cost matrix in OT problems:\n\n> The cost matrix $C \\in \\mathbb{R}^{n \\times m}$ encodes the pairwise transportation costs between source and target distributions, where $C_{ij}$ represents the cost of moving one unit of mass from source location $i$ to target location $j$. For example, when comparing two discrete probability distributions over spatial positions, $C_{ij}$ might be the Euclidean distance $\\| x_i - y_j \\|^2$ between points $x_i$ and $y_j$ in the source and target domains, respectively.\n\nWe agree that further clarification will help readers, particularly those less familiar with OT, and will incorporate it into the revised manuscript.\n\n### Note 2\n\nWe thank you for this observation, and appreciate the opportunity to clarify our notation system:\n\nIndeed, we already use iteration number as subscripts:\n\n- $H_k$: the true Hessian matrix at iteration $k$;\n- $B_k$: the approximated Hessian matrix with SPLR at iteration $k$.\n\nThe point is that we use $H_{\\Omega}^{k+1}$ to distinguish between the true Hessian and the sparsified ones. We acknowledge that better notations might be possible, but this is the clearest expression we have developed so far. When we need to express the power of $H$ (which only appears in **Assumption 3.2**), we have added a pair of parentheses around $H$. We will add a remark in the revised paper to avoid potential ambiguity.\n\n### Question 1\n\nWe thank you for this excellent question. Here are our comments:\n\n1. **Motivation**: See our response to **Weakness 1** above.\n2. **Connection to robust PCA**: The reviewer astutely observes the conceptual parallel to robust PCA. Indeed, both approaches rely on a _sparse + low-rank_ decomposition to provide a structured approximation of a given matrix. However, while robust PCA—formulated as a _Principal Component Pursuit (PCP)_ problem [1]—aims to _recover_ the low-rank and sparse components through optimization, SPLR adopts a more _constructive_ approach, providing these components via closed-form solutions. Given that robust PCA has been shown to effectively capture the essential structure of a matrix, we argue that SPLR’s approximation should likewise preserve the key information in the Hessian, ensuring a well-justified representation. Importantly, this intuition is further supported by our experimental results, where SPLR shows a superlinear-like convergence speed.\n\nYour comment has helped us recognize that this connection deserves more explicit treatment, and we appreciate the opportunity to elaborate on these important points.\n\n### Question 2\n\nWe thank you for raising this important point. As noted, the low-rank term (Eq. 7) is computed from the sparsified Hessian $H_{\\Omega}^{k+1}$​, meaning the random selection in Algorithm 2 directly influences the low-rank approximation.\n\nThe rationale behind combining Algorithm 2 with the low-rank term is to strike a balance between computational efficiency and information retention. While Algorithm 2 deliberately sparsifies the Hessian to reduce computational cost—inevitably losing some structural information—the low-rank term serves as a compensatory mechanism. This hybrid approach ensures that we maintain a meaningful approximation to the Hessian while remaining computationally tractable.\n\n[1] Candès, Emmanuel J., et al. \"Robust principal component analysis?.\" Journal of the ACM (JACM) 58.3 (2011): 1-37.", "claims_and_evidence": "Authors claim rigorous proof for linear convergence, and demonstrate super-linear convergence in examples.", "methods_and_evaluation_criteria": "Authors demonstrate the performance of the new sparse plus low rank Hessian approximated quasi-Newton iterations by applying to benchmark problems. These make sense, although there is some room for more challenging problems.", "theoretical_claims": "I only did a light check of the proofs, the claims seem to be reasonable and I did not spot any errors.", "experimental_designs_or_analyses": "The experiments are standard benchmarks (MNIST, Fashion-MNIST, ImageNet, etc), I find them suitable for analyzing convergence rates.", "relation_to_broader_scientific_literature": "The existing methods have used low-rank approximations to the quasi-Newton iterations (e.g. BFGS) or sparse approximations, but not both. This work uses the sum of both types of approximations in a systematic way to accelerate convergence.", "essential_references_not_discussed": "I do not have any suggestions.", "other_strengths_and_weaknesses": "Strengths\n\nThe proofs rely on clever combinations of elementary linear algebraic results. Presentation is clear, the proofs are well-written.\n\nWeaknesses\n\nThe overall motivations for seeking sparse plus low-rank perhaps can be explained for the benefit of the audience. Is there a more intuitive argument why such an approximation can converge faster?\n\nThe authors do not present a potential limitation to the method.", "other_comments_or_suggestions": "minor notes: \n- page 1: perhaps mention what a cost matrix is (or its form)?\n- The notation $H^{k+1}_\\Omega$ (e.g. page 5) is a bit confusing, since powers $(H_\\Omega)^{k}$ also appear in the text. Perhaps there is a better version that keeps the iteration number as subscripts?", "questions_for_authors": "What is the overall motivation for using sparse approximations? I understand the low-rank Hessian updates are well-known to be effective in general, but sparse approximations seem special. Adding a sparse matrix to a low-rank one seems to resemble robust PCA; is there any insightful relation between robust PCA and this choice of Hessian approximation?\n\nShould the addition of the low-rank part affect the choice of random selection criteria in sparse directions? i.e. Is there a reason to stick with Algorithm 2 even though now low-rank term was added in the Hessian?", "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": "xMKCmcmqDE", "reviewer": "Reviewer_Ybnz", "summary": "This paper proposes a Sparse-Plus-Low-Rank Quasi-Newton （SPLR）method for entropic regularized Optimal Transport (OT). The proposed algorithm improves the approximation of the Hessian matrix by adding a low-rank term, thus better solving the dense situation, effectively solving the entropic-regularized OT problem and additionally reducing the amount of computation. The theoretical analysis, experimental validation, and convergence guarantee in the paper are all  sufficient, demonstrating the superior performance of the SPLR algorithm.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed review and thoughtful comments. Please see our point-by-point responpses below.\n\n---\n\n### Weakness 1/Comment 1\n\nWe thank you for raising this important point about our complexity analysis. Let us clarify the theoretical aspects of our algorithm's computational efficiency:\n\n1. **Theoretical per-iteration cost**: As noted in Section 5 (Lines 267-270), each iteration's dominant cost is the inversion $(H_{\\Omega}^{k+1} + \\tau_{k+1} I)^{-1}$ in the computation of $B_{k+1}^{-1}$. According to Chapter 10 in (Shewchuk, J. R., 1994), when it is solved using conjugate gradient method, the time complexity is $O(\\rho n m \\sqrt{\\kappa})$, where $\\kappa$ is the condition number of $H_{\\Omega}^{k+1}$.\n2. **Convergence Guarantees**:\n   - Theorem 5.1 establishes linear convergence\n   - The convergence constant depends on $L$ and $U$\n\nYou are absolutely correct that the sparse pattern selection affects complexity, and we appreciate the opportunity to clarify this relationship. Our method provides explicit control over this through $\\rho$ while maintaining theoretical convergence guarantees.\n\n### Weakness 2/Comment 2\n\nWe sincerely thank you for identifying this important gap between our theoretical and empirical convergence results. The observed superlinear-like behavior indeed requires deeper investigation, and we offer the following possible causes:\n\n- *Approximation quality improvement*: As optimization progresses (especially near the solution), our SPLR approximation appears to capture increasingly accurate Hessian information. This is evidenced in Fig. 2 (columns 2-3) where the log gradient norm decreases sharply after certain iterations.\n\n- *Hybrid behavior*: The method may initially follow linear convergence (as proved) but transition to superlinear once the approximation error becomes sufficiently small.\n\nWe greatly appreciate your suggestion to pursue this theoretically, as it could lead to significant new insights about Hessian approximation methods. The gap between our current theory and observations points to exciting future research directions.\n\n### Comment 3\n\nWe thank you for this valuable suggestion. Below is a proposed comparison table we will include in the revised manuscript to clearly demonstrate the advantages of our SPLR method:\n\n| Method          | Convergence Rate                                 | Per-Iteration Cost                                                          | Memory Usage               | Hessian Approximation Type |\n| --------------- | ------------------------------------------------ | --------------------------------------------------------------------------- | -------------------------- | -------------------------- |\n| Newton's Method | Quadratic                                        | $O(n^3)$                                                                    | $O(n^2)$                   | Exact                      |\n| BFGS            | Superlinear                                      | $O(n^2)$                                                                    | $O(n^2)$                   | Dense low-rank update      |\n| L-BFGS          | Linear                                           | $O(nm)$ ($m$: memory)                                                       | $O(nm)$                    | Limited-memory low-rank    |\n| SNS<br>         | Superlinear (conjectured)                        | $O(\\rho(\\varepsilon) n^2 \\sqrt{\\kappa})$ ($\\varepsilon$: element threshold) | $O(\\rho(\\varepsilon) n^2)$ | Sparse only                |\n| SSNS            | Quadratic                                        | $O(\\rho(\\delta) n^2 \\sqrt{\\kappa})$ ($\\delta$: error threshold)        | $O(\\rho(\\delta)n^2)$       | Sparse only                |\n| SPLR            | Theoretically linear;<br>emperically superlinear | $O(\\rho n^2 \\sqrt{\\kappa})$                                                 | $O(\\rho n^2)$              | Sparse + Low-rank          |\n\nYour suggestion significantly improves our paper's ability to communicate the method's advantages, and we appreciate this constructive feedback.\n\n### Comment 4\n\nThanks for raising this question. The basic idea is as follows: since our method guarantees that the objective function value is nonincreasing, the iterates $\\{x_k\\}$ are restricted to the level set $D = \\{x: f(x) \\leq f(x_0)\\}$, which can be proved to be compact. Since $H(x)$ is continuous on $x$, it must be bounded on $D$. Of course, the bounds may depend on the initial value $x_0$.\n\n### Comment 5\n\nWe thank you for this insightful technical question. Our rank-2 approximation is motivated by the rank-2 modification scheme in the BFGS quasi-Newton method, and can indeed be naturally generalized to higher ranks with all theoretical gurantees. The question highlights an important implementation flexibility that deserves more explicit treatment and we will adjust our manualscript accordingly.\n\n### Comment 6\n\nWe thank you for pointing out this typo. It should be $auu^T + bvv^T$.", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "Yes", "theoretical_claims": "Yes", "experimental_designs_or_analyses": "Yes", "relation_to_broader_scientific_literature": "Quasi-Newton Method", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "**Strengths: **\n\n(1) The authors effectively enhanced the algorithm's efficiency by integrating a low-rank matrix into the sparse approximation framework for Hessian, with rigorous theoretical analysis demonstrating both convergence properties and computational complexity reduction, well-supported by theoretical analysis.\n\n(2) The empirical tests conducted on both synthetic and actual datasets, highlight the practical superiority of SPLR compared to conventional techniques such as Sinkhorn, L-BFGS, and SSNS. The ablation study, detailed in Appendix A.1, convincingly confirms the essential role of the low-rank component.\n\n**Weaknesses:**\n\n(1) It is mentioned that the algorithm has a lower time complexity, but it relies on the choice of sparse matrices during the computation process, and the paper does not provide specific theoretical analysis.\n\n(2) The paper mentions that during the experimental process, the algorithm can achieve super-linear-like convergence speed, but in the theoretical proof, the algorithm only reaches linear convergence speed. What is the cause of this phenomenon? Is it possible to further derive corresponding results in theory?", "other_comments_or_suggestions": "(1) It is mentioned that the algorithm has a lower time complexity, but it relies on the choice of sparse matrices during the computation process, and the paper does not provide specific theoretical analysis.\n\n(2) The paper mentions that during the experimental process, the algorithm can achieve super-linear-like convergence speed, but in the theoretical proof, the algorithm only reaches linear convergence speed. What is the cause of this phenomenon? Is it possible to further derive corresponding results in theory?\n\n(3) Although the paper provides a brief introduction to existing quasi-Newton algorithms, it does not compare their convergence rates with those of classical algorithms. It is recommended to add a comparison table to illustrate the superiority of the algorithm.\n\n(4) How is the boundedness of H(x) ensured during the iteration process?\n\n(5) In the article, a rank-2 approximation term was added to achieve a faster convergence rate. Can this approximation term be of a higher rank? \n\n(6) Is there an error in line 16 of Algorithm 1?", "questions_for_authors": "Please see the above comments and suggestions.", "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": "Kk178ixlgN", "reviewer": "Reviewer_aArm", "summary": "This paper is about solving entropic-regularized optimal transport problem via a quasi-Newton method. Entropic-regularized OT problem has applications in machine learning, but its solution is difficult to find. The paper first presents some theoretical analysis Hessian sparsification, which can be used to solve the OT problem in a more tractable way. The paper then presents a low-rank approximation to the Hessian matrix to facilitate computations. The paper later gives theoretical analysis of the method and validates its performance via experiments.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed review and thoughtful comments. Please see our point-by-point responpses below.\n\n---\n\n> Not always, as ... on machine learning.\n\n> Biggest weakness is ... of this conference.\n\nWe sincerely appreciate your feedback on our manuscript. Your observation regarding the limited discussion on machine learning is invaluable and has encouraged us to further explore this perspective to enhance our work.\n\nTo provide a more comprehensive overview, we have reviewed recent advancements in computational optimal transport (OT) within the machine learning community, particularly those presented at ICML. Our findings indicate a substantial body of literature [1-4] that focuses on this topic, underscoring its relevance within the field.\n\nFurthermore, the [ICML 2025 Call for Papers webpage](https://icml.cc/Conferences/2025/CallForPapers) includes \"Optimization (convex and non-convex optimization, matrix/tensor methods, stochastic, online, non-smooth, composite, etc.)\" as a topic of interest, highlighting the strong connection between our research and ongoing discussions in the machine learning community.\n\nOnce again, we sincerely appreciate your insights, and we would add more discussions on the importance of OT optimization in the revised paper.\n\n**References**\n\n[1] Dvurechensky, Pavel, Alexander Gasnikov, and Alexey Kroshnin. \"Computational optimal transport: Complexity by accelerated gradient descent is better than by Sinkhorn’s algorithm.\" ICML, 2018.\n\n[2] Lu, Haihao, Robert Freund, and Vahab Mirrokni. \"Accelerating greedy coordinate descent methods.\" ICML, 2018.\n\n[3] Lin, Tianyi, Nhat Ho, and Michael Jordan. \"On efficient optimal transport: An analysis of greedy and accelerated mirror descent algorithms.\" ICML, 2019.\n\n[4] Guminov, Sergey, et al. \"On a combination of alternating minimization and Nesterov’s momentum.\" ICML, 2021.\n\n> In Eq. (5), ... is it needed?\n\nThank you for your insightful question. We cannot guarantee that the low-rank approximation is always necessary. However, based on our numerical experiments, particularly in Figures 5 and 6, we have observed that adding the low-rank term generally leads to faster convergence of the algorithm. Specifically, when the curvature condition $(s^k)^T y^k > \\varepsilon \\|y^k\\|^2$ holds, we consistently find it beneficial to include the low-rank term.\n\n> In Eq. (6) ... behind the method.\n\nThank you for pointing this out! The correct notation should be $u = y_k$, and we will fix it in the main text.\n\nAs for the variables $(a, b, u, v)$ in Eq. (6) and (7), they are mostly notational: $u$ and $v$ represent two low-rank components, derived from the BFGS quasi-Newton method, and $a$ and $b$ are their respective scaling coefficients, determined by the secant equation (Eq. 6).\n\nThese terms collectively help ensure a more stable and efficient update, improving the convergence behavior of the method. We will expand the explanation in the main text to better highlight the motivation behind these choices.\n\n> The experimental results ... in the first place.\n\nThank you for your insightful observation. Indeed, our results confirm that the SPLR method is as good as alternatives in virtually all scenarios, and significantly outperforms others when the Hessian is sparse. This demonstrates the robustness of SPLR in worst cases and its high efficiency in favorable cases.\n\nThis behavior is due to the sparse-plus-low-rank design, which ensures stable and efficient performance across different Hessian structures.\n\n> Figure 4 seems to ... shows similar patterns.\n\nWe appreciate your observation about the impact of the regularization weight $\\eta$ on convergence speed. Indeed, as shown in Figures 4 and 11, the choice of $\\eta$ significantly affects all methods, including ours. While our paper focuses on optimizing the OT problem *given* a fixed $\\eta$, we acknowledge that selecting $\\eta$ is critical in practice—though this lies outside our current scope. Investigating systematic strategies for choosing $\\eta$ efficiently could be an interesting direction for future work.\n\n> Another concern is ... elaborate on this?\n\nThank you for your valuable observation. We appreciate your concern regarding the oscillation in the Log10 Gradient Norm observed in SPLR. However, we would like to clarify that such oscillations do not pose a risk of early or incorrect stopping. In fact, the gradient norm serves as a rigorous stopping criterion. Specifically, when the gradient norm is reduced below a predefined threshold, we can theoretically prove that the optimality gap is also bounded by a corresponding threshold. This ensures that the stopping condition is reliable and provides confidence in the convergence of the method.\n\nOn the other hand, our method guarantees that the dual objective function value is nonincreasing across iterations, so if the stopping criterion is on the function value, it will not have any oscillation.", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "Yes", "theoretical_claims": "Yes", "experimental_designs_or_analyses": "Not always, as the paper barely touches on machine learning.", "relation_to_broader_scientific_literature": "It is related to relevant literature in using OT in machine learning.", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "The strength includes a comprehensive analysis of using sparse Hessian and a low-rank approximation in solving entropic-regularized OT problem and a large amount of experiments. Another strength is to set the \\tau dynamically so when the solution is close to the minimum, \\tau is updated more gradually.\n\nBiggest weakness is the seemingly lack of relevance of this paper to machine learning. While OT is indeed widely used in machine learning, this paper barely covers any machine learning techniques or applications. This makes the paper somehow out of scope of this conference. \n\nIn Eq. (5), is the low-rank approximation always needed? If not, when is it needed? \n\nIn Eq. (6) and (7), how is u, actually u_k, determined? It seems it was never defined in the main text, though it was given in the supplementary material, and it seems u_k was a typo. And can the authors give some explanation about the meaning of a, b, u, v in Eq. (6) and (7)? That will help readers see the motivation behind the method. \n\n\nThe experimental results are not always clear to me. In Figure 1, especially second row, it seems the proposed SPLR method had similar run-time vs log10 Gradient norm curves as several existing methods. \n\nIn Figure 2, when the Hessian is not so dense, the advantage of SPLR method is more obvious. Considering Figure 1 and Figure 2 together, it seems to indicate the SPLR method works better if Hessian is not so dense in the first place. \n\nFigure 4 seems to indicate that the regularization weight is critical in determining each method's speed of convergence. While SPLR gave better performance than some methods, it converged slower than a few other methods. If so, that means a careful selection of the regularization weight \\eta is necessary, and this may increase the computational load of the proposed method (as well as the other methods). Figure 11 shows similar patterns. \n\nAnother concern is that it seems SPLR experienced large oscillation in reducing Log10 Gradient Norm, though some other methods had the similar behavior, but this may raise the risk of early or incorrect stopping of iterations. Can authors elaborate on this?", "other_comments_or_suggestions": "Please see strengths and weaknesses.", "questions_for_authors": "Please see strengths and weaknesses.", "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:

Meta-review: The paper introduces a novel quasi-Newton method that exploits sparse plus low-rank Hessian approximations for entropic OT, demonstrating superior empirical performance over existing methods (Sinkhorn, L-BFGS, etc.) in both synthetic and real-data experiments. While some theoretical gaps remain (e.g. dimension-dependent rates, superlinear convergence in practice), the method’s practical efficacy and solid analysis leads me to recommend that the paper should be accepted.
Decision: Accept (poster)