A Truncated Newton Method for Optimal Transport
Authors: Mete Kemertas
Amir-massoud Farahmand
Allan Douglas Jepson
Venue: ICLR 2025 Poster (2025)
Keywords: Computational optimal transport
numerical optimization
numerical linear algebra
Primary area: other topics in machine learning (i.e., none of the above)
TL;DR: A high-performing, high-precision, GPU-parallel optimal transport solver based on truncated Newton methods.

Abstract: Developing a contemporary optimal transport (OT) solver requires navigating trade-offs among several critical requirements: GPU parallelization, scalability to high-dimensional problems, theoretical convergence guarantees, empirical performance in terms of precision versus runtime, and numerical stability in practice. With these challenges in mind, we introduce a specialized truncated Newton algorithm for entropic-regularized OT. In addition to proving that locally quadratic convergence is possible without assuming a Lipschitz Hessian, we provide strategies to maximally exploit the high rate of local convergence in practice. Our GPU-parallel algorithm exhibits exceptionally favorable runtime performance, achieving high precision orders of magnitude faster than many existing alternatives. This is evidenced by wall-clock time experiments on 24 problem sets (12 datasets $\times$ 2 cost functions). The scalability of the algorithm is showcased on an extremely large OT problem with $n \approx 10^6$, solved approximately under weak entropic regularization.

Published as a conference paper at ICLR 2025
A TRUNCATED NEWTON METHOD FOR OPTIMAL TRANSPORT
Mete Kemertas
Amir-massoud Farahmand
Allan D. Jepson'
University of Toronto, Polytechnique Montréal, Mila Quebec Al Institute, Vector Institute
ABSTRACT
Developing a contemporary optimal transport (OT) solver requires navigating trade- offs among several critical requirements: GPU parallclization, scalability to high- dimensional problems, theoretical convergence guarantees, empirical performance in terms of precision versus runtime, and numerical stability in practice. With these challenges in mind, we introduce specialized truncated Newton algorithm for untropic-regularizee OT. In addition to proving that locally quadratic convergence is possible without assuming a Lipschitz Hessian, we provide strategies to maximally exploit the high rate of local convergence in practice. Our GPU -PU-parlllel algorithm exhibits exceptionally favorable runtime performance, achieving high precision orders of magnitude faster than many existing alternatives. This is evidenced by wall-clock time experiments on problem sets (12 datasets x cost functions). extremely large OT problem
The scalability of the algorithm is showcased with n se 10 solved approximately under weak entropic regularization.
INTRODUCTION
The optimal transportation problem has long been cornerstone of various disciplines, ranging from physics (Bokanowski & Grébert, 1996; Léonard, 2012; Levy et al., 2021) to machine learning and computer vision (Ferns et al., 2004; Pitic et al., 2005; Guirajani et al., 2017; Genevay et al., 2018). Traditional approaches (Pele & Werman, 2009; Lee & Sidford, 2014), while exact and theoretically robust, encounter significant computational hurdles in high-dimensional settings. The (re-) ree)iintroduttoon of entropic regularized OT (EOT), as pioneered by Cuturi (2013), has mitigated challenges in scalability by regularizing the classical problem, thereby enabling solutions via the GPU-friendly Sinkhorn- matrix scaling algorithm. This advancement has yielded substantial speed improvements, making it several orders of magnitude faster in high dimensions than traditional solvers. However, EOT methods necessitate a delicate balance between regularization strength and convergence speed, a trade-off that can compromise the precision of the solution.
Despite significant recent progress towards improving this tradc-off, many state-of-thhe-a solvers still struggle to outperform aggressively tuned Sinkhorn iterations in practice (Jambulapati al., 2019; Lin et al., 2019). While they offer superior theoretical guarantees, their practical performance often less compelling, particularly in terms of speed and scalability. Existing algorithms either suffer from high computational complexity or fail to leverage modern hardware capabilities, such as GPU parallelization, effectively. To bridge this gap, we develop new algorithm that remains numerically stable and converges rapidly even at extremely weak regularization levels, thereby enhancing precision in practice. By simultancously exploiting the inherent parallelism of GPUs and superlinear local convergence of truncated Newton algorithms, method scales effortlessly to high-dimensional problems, offering a pragmatic yet theoretically sound solution to the OT problem.
Our contributions are follows: (i) develop a specialized (linear) conjugate gradient algorithm for obtaining an approximation of the Newton direction for the EOT dual problem and analyze its convergence properties, (ii) we the approximate (truncated) Newton direction in conjunction with a helper routine to develop solver for the EOT dual problem and prove its superlinear local convergence, as well as per iteration computational cost, (iii) we develop an adaptive temperature annealing approach, based on the MDOT framework of Kemertas et al. (2025), to maximally
"Correspondence kemertas@@es.orootto.ccuuu Code available github.cco//mmeeee githan...ommmmeeekeeettttsmaooo

---

Published as a conference paper at ICLR 2025
exploit this fast local rate, and finally (iv) present compelling empirical results via wall-clock time benchmarking in a GPU setting against a large suite of alternative algorithms n = 4096 dimensions.
2 BACKGROUND AND RELATED WORK
Notation and Definitions. In this work, we are concerned with discrete OT. A. C RR denotes the (n-1) --iiple.,. The row sum of an n X n matrix P is given by r(P) = P1 and the column sum by c(P) == 11 Given target marginals r, c E A., the transportation polytope is written as u(r,c) c) = [P R220 rr(P) c(P) c}. Division, exp and log over vectors or matrices indicate element-wise operations. Vectors in TR are taken to be column vectors and concatenation of two column vectors 2: y is (z, y). Elementwise minimum and maximum of vector x: is written as min and maax Matrix and vector inner products alike are given by (:. An n x diagonal matrix with Φ E I= along the diagonal is written as (xx). We write x²(y]a) a+) for the x" -dvveregee given 2 lyy - x: For the square root of x x2 a)), we write x(yr(a)) with a slight abuse of notation. The Shannon entropy of
by (ar, (y/an (x)=²² )
Ei ry rg- We write DKL(aly = E; xz log(z;/u)) >; y/ - - 2; for the KL divergence between æ, E 150. We +),)).B,3
A, is denoted ((r)
denote LogSumExp reductions along the rows and columns of X by (X), = log (exp(X}1) and (X) = log (exp(X)) }1)
2.1 OPTIMAL TRANSPORT AND ENTROPIC REGULARIZATION
We study the discrete optimal transport problem, formulated as the following linear program:
minimize Peu((r
(P,C),
where we assume the cost matrix has entries Ci E |0, 1]. Cuturi (2013) re-popularized EOT, showing that entropic regularization can help quickly approximate the solution of (1) on GPUs:
-
1 H(P.,
minimize P E u(r,
(P,C)-
where the regularization weight > E Rz is called the temperature. It can be shown with ease that since the objective in (2) is strictly convex in P. problem (2) has a unique solution of the form
P(u, 0; exp(u1" +1v CC..
Using the form of the solution of (2), the following unconstrained dual problem can be solved instead:
minimize u, R
g(u, 0;
Σ P(u, u, r;
(u, (v. c),
where we keep the constant -1 convention. Solving (4) given initial u, v amounts Bregman projection onto U(r c) in the sense that P(ui" v²) arg min peure et al., 2025). Noting that Vng -r(P)= and V.g c(P) we write:
v)) (Kemertas
V2g D(r(P)) PT
P D(c(P))
(5)
2nx2n
where the Hessian is positive semi-definite (PSD) with one zero eigenvalue!'
Related Work. The SK algorithm has long been known to enjoy an exponential convergence rate for minimizing (4) (Franklin & Lorenz, 1989; Knight, 2008). However, at low temperatures, this fast rate does not predict non-asymptotic behavior well due to a large constant. Altschuler et al. (2017) provided a simple analysis, in which they proved rate =3n²==3) for the SK algorithm, where (P- C) S e. simple routine for rounding near-fea plans onto (r e) was introduced and now widely adopted. They also proposed new algorithm. Greenkhorn, and showed matching complexity bound. Unlike SK, Greenkhorn scales one greedily selected row/column at a time, which limits GPU utilization unless is extremely large. The complexity bounds for Sinkhom and
Since r(P) P1 D(r(P))1 and c(P) P ((c(())) have V (111111)
-1) 0.
2

---

Published as a conference paper at ICLR 2025
Greenkhorn were later improved to 3 (Dvurechensky et al., 2018; Lin et al., 2019). However, our experiments suggest Sinkhorn typically behaves like 1(n²==1). Moreover, Kemertas al. (2025) showed it can enjoy better performance at lower temperatures if tuned.
Dvurechensky et al. (2018) proposed an Adaptive Primal-Dual Accelerated Gradient Descent (APDAGD) algorithm for solving the dual EOT problem (4). Lin et al. (2019) provided refined rate of ( for APDAGD and proposed generalization APDAMD, which applied mirror descent to (4). The complexity for APDAMD was shown to be Vc/e,, where E (0, n) is a constant. Following Altschuler et al. (2017), Lin et al. (2019) measured speed in terms of the number of row/col updates in their experiments (rather than wall-clock time) and only considered a strong regularization (low precision) setting. Targeting higher precision, Jambulapati et al. (2019) proposed an algorithm for the OT problem that is not based on entropic regularization, with complexity O(n²==1). While this rate is theoretically state-of- of-theeart, Jambulapati et al. (2019) noted that Sinkhorn iteration, when aggressively tuned, outperforms all other methods empirically (including their own). Guminov et al. (2021) proposed an Accelerated Alternating Minimization (AAM) algorithm, combining Nesterov's momentum and Sinkhorn-type block coordinate descent, with complexity ((n//2=-1)
While APDAMD applied mirror descent to (4), MDOT of Kemertas et al. (2025) applied it to (1) and recovered connections to temperature annealing methods (see details in Sec. 2.2), such as those of Schmitzer (2019) and Feydy (2020). For instance, Alg. 3.5. of Feydy (2020) takes single Sinkhorn update every time the temperature is decayed; while this can compute rough approximations quickly at high temperatures, a single Sinkhorn update is insufficient for keeping the dual objective value in check, so that their approach hits a precision wall as we empirically show; see also Xie et al. (2020). Ballu & Berthet (2023) derive a similar algorithm, but they guarantee convergence by maintaining a running average of plans P' computed this way. While effective at low precision and easy to implement on a GPU, this algorithm exhibits 2 dependence to ours is the work of Kemertas et al. (2025), as we build on MDOT. In addition to the temperature anncaling framework, Kemertas et al. (2025) proposed an algorithm (PNCG) to minimize (4) at each new value of the temperature. based on a non-linear conjugate gradient method (Fletcher & Reeves, 1964). The approach introduced here has several benefits over PNCG, including added ease of theoretical analysis, faster runtime in practice and minimal line search overhead. While second order methods have been considered for OT (Mérigot, 2011 I; Blondel et al., 2018), they have not been implemented on GPUs with strong empirical performance in high dimensions to our knowledge. Indeed, Tang et al. (2024) also developed a 2nd order method recently, but their Hessian sparsification strategy is more amenable to a CPU setting, and as such was only tested on CPUs for n = 784.
error. Most closely related
2.2 TEMPERATURE ANNEALING AS MIRROR DESCENT
well-known strategy to deal with the difficulty of solving (4) under weak regularization is annealing the temperature y in (3) gradually towards zero. When viewed as mirror descent on (2), temperature annealing strategies (c.g., see Schmitzer (2019)) amount to a particular initialization of the dual variables (u, v) in successive instances of (4) given approximate solutions prior 7 (Kemertas et al., 2025). Each dual problem (4) a given (+++1) for t
2 1 is warm-started in some neighborhood of the solution given some near-optimal =()) (!)) E TR The MDOT framework of Kemertas et al. (2025) specifically initializes z (++1) via a Taylor approximation with respect to under backward finite differencing. Further, they proposed to use a more stringent tolerance O(>-P) for [[all,, (given some p 2 1), whereas prior work used 0(> (Altschuler et al., 2017; Lin et al., 2019). Their tuned choice 1.5 was shown to improve performance under weak regularization.
The pseudo-code for MDOT is shown in Alg. 1, where some routines are defined with "..." as a placeholder for extra parameters that may be required by specific implementations. Given (0) in each iteration: 2 1, the algorithm picks a tolerance Ea for the dual gradient norm (L4). Then, marginals r, are smoothed in L5 for numerical stability or improved convergence by mixing in the uniform distribution (with a combined weight of at most 2a/2) to stay away from the boundary of U((r c). L6 initializes the transport plan P to be the independence coupling FēT ii... the solution of (2) for - O over U(5 é). In 1.7, minimizing (4) to ea/2 tolerance for the smoothed marginals guarantees (V)(z;)))l S Ea by triangle inequality. Next, we add an AdjustSchedule routine in L8, whereas Kemertas et al. (2025) used a fixed decay rate, setting g (+++1) t- q(º) After the temperature decay (I.9), the dual variables are warm-started in L10 via an approximate 1st order expansion for the next
3

---

Published as a conference paper at ICLR 2025
iteration. Finally, the plan given by (3) is rounded onto U(, c) in L14 via Alg. 2 of Altschuler et al. (2017). Then, given that [[Vg(z: S "/r min(H(r), H (c)) the user is guaranteed error (P PP. .C) < 2 min(H(r) (c)) in the worst case (assuming r(P) = c(P) at loop termination). In Section 3. we develop and tune specific algorithms in tandem to carry out L5, 7 and 8 of Alg. 1 maximize efficiency, and maintain GPU parallelization and numerical stability.
2.3 TRUNCATED NEWTON METHOSS
For minimizing (4) in L7 of Alg. I, we will rely on truncated Newton methods, which are briefly reviewed here (see Ch. 7 of Nocedal & Wright (2006) and Nash (2000) for a survey). While Newton's 's method for non-linear optimization seeks exact solution dg of the linear sys-
Algorithm 1 rr,
1,4 1: t 1, done +- false, 7(1) + A g()) + 9,,,00) t
2: while not done do
3: done t- 7(1) = 7/f
4:
Ed t min
(7
5:
r, < Smoothh G, Ed/ ea/2: 2:...)
tem g d₂ yig. at each opti-
6:
if then x0)) < (log r. log c) 211) t 20)
mization step k (typically combined
7:
(() <- BregmanProject(= (t)
C,F, c, Ea
with line search or trust-region methods), truncated Newton meth-
g (++1) < AdjustSchedule(()
9:
7 ((+1) + (( V
ods find an approximate/inexact so-
lution by "truncating" an iterative 10:
x(1-1) + z(1) -(11 2(0) 2(0-2((-11) x(1-1))
solver of the linear system (L.S), e.g..
11:
1-t+1
linear conjugate gradient (linear
12:
end while
CG) algorithm. The particular ter- mination criteria for the LS solver dictates the order of convergence; more stringent criteria yields higher order of convergence (up to quadratic). but requires more iterations for the LS solver (inner-most loop).
13:
(u,v) + z (tt-1), + exp + +1, - -n
14: Output P t- (Pr..c
In particular, define the residual eg V gk de Vgk of the system in step le for the approximate Newton direction dy. Given 77 E 10, the LS solver is terminated when
| || S mh Vg |I
(6)
for some norm. If the Hessian is continuous in some neighborhood of the solution and we start sufficiently close to the solution, we have V || S (7ı + o(1) g/lll Then, choosing the
forcing sequence {m such that 7/s- = V convergence of order r assuming lime-+oe M) 0 (Nocedal & Wright, 2006).
for r E (1.2), one obtains Q-superlinear local
3 MIXING TRUNCATED NEWTON & SINKHORN FOR BREGMAN PROJECTIONS
This section is organized as follows. In Section 3.1, we first develop technique for obtaining the truncated Newton direction with convergence guarantees (Alg. 2) and introduce helper routine (Alg. 3) to improve the convergence rate of this algorithm. In Sec. 3.2, we integrate this approach with backtracking line search to arrive at Bregman projection algorithm (I.7 of Alg. 1). Its per-step cost and local convergence properties decay schedule (L8 of Alg. 1) is proposed to maximally exploit the high rate of local convergence and its use is empirically demonstrated. Lastly, in Sec. 3.4, we discuss precautionary measures for numerical stability of this technique in practice via marginal smoothing as in L5 of Alg. 1.
shown theoretically. Then, in Sec. 3.3 an adaptive temperature
3.1 THE DISCOUNTED HESSIAN AND THE BELLMAN EQUATIONS
Suppose we are interested in finding an approximate solution to the Newton system for the dual (Bregman projection) problem (4) with a linear CG solver. In this section, we propose particular positive-definite (PD) approximation of the PSD Hessian in (5) for two reasons. First, although linear CG should converge in theory despite the presence of a zero eigenvalue (Axelsson, 2003), infinitesimal numerical errors along this zero-eigendirection, as well as other unknown near zero eigendirections, may compound and lead numerical instabilities due to machine precision limits. Second, since our
4

---

Published as a conference paper at ICLR 2025
PD substitute for the Hessian allows for the of matrix inversion, this is convenient for obtaining theoretical guarantees. This approach also yields connections to reinforcement learning that may be of separate interest. Henceforth, define the p-discounted Hessian for p [0, 1):
(D((())
vpP
V99()) = VDPT
D(c(P))
By the Gerschgorin Circle Theorem, all eigenvalues (V 2 (1-vD) min(r( Plmin el m())mn)) are positive, so that (pg(p) is invertible. The block matrix inversion formula yields:
V ggg(p)~~~ (- (x(
-VDF,(D)1PP
(8)
-vPPeFr(())-
Fe(p)
where we used the following definitions:
Pr = D(r(P))- Pe = D[=[((())-1p
Pre = Pe Per = PeP (p(p) == D(r(P)) plrc), Fe(p) = D(c(P))(L pPcr).
(9)
Here, Pr, Pa: Pre and Per are irreducible row-stochastic matrices since all their entries are positive and we have P,I ((r(P)) P1 = D(r(P))- rr() = 1 (likewise for Pa). Now, multiplying on the left by (8) both sides of the discounted Newton system (pg(p) d = Vg and re-arranging, we arrive at the well- Bellman equations central in reinforcement learning:
du suv + pPrreuu de svu pPerdp
(10)
where Pre and Per serve as "transition matrices", and "reward vectors" are given by
8pu = - D(cc(P)) That is. the discounted Newton direction d = (dn, dp)) corresponds to the fixed point of the Bellman equations (or, "state-value function" of the finite Markov reward processes) in (10) and can be written as two vavariible PD linear systems rather than a single 2n-variable PSD linear system. Further, without loss of generality, assuming that Veg = 0 (or equivalently c(P) c, for example, following single Sinkhorn update) yields a more intuitive understanding and a practical advantage. Using the form of the inverse discounted Hessian in (8), it can shown that in this case (10) reduces to: du = (r/r(P)- - 1) + pPredu.
Suv (V
TPP/Vus) (11)
de = - -Vpaddu
(12)
where the second system is solved via single matrix-vector product (effectively for free), thus reduc- ing the problem further to a single n-variable PD linear system. Moreover, an intuitive interpretation of the reward vector emerges, as the system now assigns reward r:/r(P); 1 to row index (state) i. Next, we provide a theorem on the sufficiency of solving the discounted Newton system as proxy for the undiscounted system (see Appx. 1.1.3 for
technical version of the following).
Theorem 3.1 (Foreing sequence under discounting). Assuming e = c(P) and d, = Podu define residuals e., (p):=Pi (p)dd + V.ng (cf. (12)), and e == V g d + Vg (i... the Newton residual), For every 7j < // 2. 3po E [0, 1) such that Vp E [po, 1).
[|][(p)] I 7 Vgll, = |lell, >n V/ll.
(13)
In other words, discounting can be adopted while still enjoying the local convergence guarantees of truncated Newton methods (discussed Sec. 2.3) for the original problem.
To find a suitable discount factor, we propose Algorithm 2, which anneals 1 - p. approximately solving a sequence of linear systems satisfy the forcing inequality (6). Note that CG in L3 is terminated when the L1 norm of the residual is below ń /444. While the specific selection in L3 only for simplicity, the decay factor 4 in L4 minimizes a theoretical upper bound on the number of operations until convergence (see proof of Thm. 3.2 in Appx. 1.1.4). The following provides a rate on the overall cost of obtaining the truncated Newton direction with this approach.
17/4
Theorem 3.2 (Convergence of Algorithm 2). Suppose r r(P)min 2 ca/(4n) given Ea > 0 and each CG initialized with d -00 Alg. 2 terminates in -0.
step of Alg. 2 runs diagonally-
(1 ()))r -1-x2)n [[]([)) -| log= -1/2
δ
n²
(14)
operations, where X2 < is the 2nd largest eigenvalue of Pre and μ < 1 its smallest diagonal entry.
We refer the reader to Appx. A.1.3 for an intuitive description of the process represented by the stochastic matrices Pre and Par, and their technical properties.
5

---

Published as a conference paper at ICLR 2025
The proof of this theorem uses (i) the ((vr) convergence of linear CG (Shewchuk, 1994), where condition number
Algorithm NewtonSolve (V Pre: r(P).
O((1 p) p)
I: p t O, du t- - DD Vug
in our setting, and (ii) an observation from the proof of Thm. 3.1 that
2: while | (1)d + V u6 | > 77
V do
3:
du <
(F (p), V ug n/4)
ve
(I (1-
at termination
4:
P p=1- p)/4
in the worst-case. However, as we discuss in more detail and validate empirically in Appx. E, this worst-case dependence on (1 A2)
5: end while
6: Output du
can be overly pessimistic; (1 p) can be much better behaved than the spectral gap in practice. In other words, the Newton system can be discounted more aggressively than the worst-case analysis allows, thereby mit- igating possible large values of (1 x2)-1
Algorithm 3 ChiSinnhormn(u,,, u, ,,
I: while x (r(rr(P) > >x do u u + logr log r
2:
3:
v + log c (uE, CC
4: logr(P) + u + LSE, (In -00)
Furthermore, the convergence of CG can be much faster when preconditioned eigenval- ues are tightly clustered R (Nocedal & Wright, 2006). The rate in (14) captures this added efficiency only for the case when Pre ~ I so that its smallest diagonal entry μ = 1; for example. if plan P is approximately a one-to-one mapping as in the Monge discrete matching problem (Brezis, 2018).
5: end while
6: Output u, v, r(P)
Next, to control the possible dependence in (14) of the ratio propose run Sinkhorn iteration before Alg. 2 is called until (r(r(((P)) suitably bounded. Given dual variables (u, v) and current row sum r(P) where c = c(PP) by assumption, Alg. 3 performs this auxiliary task while maintaining e = c(P).
(|(p)))))///()---))r/)))))) on Ed. we
Lemma 3.3 (Convergence of Algorithm 3). Assuming that ||r(()/rllo < 0o and (ard|r//((())) 0o.
Algorithm 3 converges in O(n²/e, /ex) operations.
3.2 PROJECTING ONTO THE FEASIBLE POLYTOPE
Combining Algorithms 2 and 3 with backtracking line search, we arrive at the TruncatedNewton- Project algorithm shown in Alg. for solving (4). The choice of r) in L5 and the requirement that (P(r)(()) S Ed
2/5 by choosing Ex Ed in L4, to control the ratio in (14), together yield the following corollary Thmm. 3.2 and Lemma 3.3 (see proof in Appx. A.1.5).
2/5
Corollary 3.4 (Per-step Cost of Algorithm 4). If the backtracking line search in Alg. 4 converges in S iterations, then iteration of Alg. 4 costs Ónn + + 2 X2)-1 operations, where 22 < 1 the 2nd largest eigenvalue of Pre defined as in (9) and evaluated at u, v (cf. (3)).
In the next section, we outline an adaptive temperature annealing strategy to minimize the added cost of line search by initializing close the solution. First, we pause for the next theorem, showing that Alg. 4 enjoys local quadratic convergence if 1 = g(l)) as in L5. Theorem 3.5 (Per-step Improvement of Algorithm 4). Given descent direction d = (du, -Pdd such that ||ell = V2 g + V S gu ler a E (0, 1] be the step size found via backtrack- ing line search in the bth step of Alg. 4. Then, Vgk++ = Vg(u + adu, a duu) satisfies
1gk+1l1 < (1 a + an) Vgxll, || + avao O(
(15)
The result differs from typical quadratic convergence results in two important ways: (i) we did not assume Lipschitz Hessian, but instead leveraged the Armijo condition and the specific form of our descent direction, and (ii) we bounded the L1 norm of the gradient (rather than Lg), which is more commonly used in optimal tramsport algorithms as stopping criteria. Given (15), we select 1 in L5 of Alg. 4 with the local quadratic rate mind, but avoid over-solving the system when V/ is already close the target ea by taking 7) to be the maximum of / and 0.8ea/ ||
3.3 INITIALIZING NEAR THE SOLUTION VIA ADAPTIVE MIRROR DESCENT
In this section, we outline a practical strategy for initializing the dual problem (4) sufficiently near the solution that (i) the cost of line search is minimized (a = 1 is almost always admissible).

---

Published as a conference paper at ICLR 2025
Algorithm 4 TruneateedNewrooPPoojccttl C,r,
1: v + log LLE, -C)
> Ensure c(P)
2: logr(P) t- + (1, - -7C)
3: while VVll1 = - > sa do
4:
u, v, r()) - ChiSinkhornii( 0, 7, C,r,c,
P Choosing 8x
2/5
5:
T - Vgll v V.. Vgll,
P Expect || Vgg+illl |I
77 V
6:
du f NewtonSolve( V ug, Pe: (P).
>> See Algorithm
7:
d. f------- -odd
As in (12).
8:
a+1
Initial guess for step size.
9:
c(P) f- v ad,, LSE. (u ad,, , -7C)
10:
while ||c(P)]|
0.99a (300rr(Aoo V Vng, du,) do
b Armijo condition. See Appx. A.1.6.
11:
+- 0.5c
Backtracking line search.
12:
log c(P) + v + ad, + LSE adu) -CC)
end while
14:
u u u add, t- + adp
ot e log logc c(P)
D Ensure c(P)
16: logr(P) - LSE, (1, 7C)
17: end while
18: u logr logr(P)
Since rPP) is readily available, take Sinkhorn step.
19: Output u, v
and (ii) the last term in (15) is negligible. In Table I, we observe that with a sufficiently slow fixed temperature decay schedule (i.e., setting (+++1) < ((t) L8 of Alg. 1), the extra cost of line search (as well as ChiSinkhom) disappears almost entirely (given small enough (1)). However, with too slow temperature decay schedules (q too close to 1), the overhead due to relatively costly LSE reductions (in lines 1, 2. and 16 of Alg. 4) may slow down the overall algorithm (see 21/8 in Table 1). To eliminate the need for tuning q. develop an update rule adjust the schedule in L8 of Alg. I; we seek to minimize the number of mirror descent steps while staying in the superlinear convergence zone for consecutive instances of problem (4). To this end, we first define the following parameter of interest, which is the ratio of the actual reduction in gradient norm to the predicted reduction given by (15) for a (in the ideal case, dropping the last term) step of Alg.
V |1
V 19k+1l11
(16)
(1 )7) g
This formula for ò% is inspired by trust-region methods, which update the trust-region size based on the ratio of actual to predicted decrease under a model of the objective (rather than the gradient norm
here) (Nocedal & Wright, 2006). Let 5(t) be the smallest ōg among the iterates of Alg. at step / of Alg. 1. In L, of Alg. 1, we heuristically set
(++1) + 22 (1
if 8 (t) 0.95
((+1) +- (((t)
if (()
min < 0.8 otherwise.
(1+1) t (1)
q 2//8
2//4
2//2
21
22
Adaptive
Subroutine
NewtonSolve
1870
1725
1744
2223
3281
1831
Line Search
8
48
226
24
ChiSinkhom
8
8
8
96
8
Mirror Descent
105
53
27
14
23
Total (operations)
4623
3315
2760
3178
4369
2795
Total (seconds)
2.84
1.99
1.56
1.76
2.41
1.55
Table 1: Breakdown of number of O(n²) operations per subroutine and total wall-clock time for the upsampled MNIST dataset (n = 4096) with L1 distance cost. Alg. is called with 7 = 2-, = 218 and p 1.5. For the adaptive approach, initially. Results show median over 100 problems.
7

---

Published as a conference paper at ICLR 2025
Ratio of Actual-to- Gradient Norm Reduction Temperature Adaptive
Fixed
15
3.0
Decay rate
2
2
2
2
3.0
Y
Figure 1: Ratio 6g of actual to theoretically predicted reduction in V ll per step for fixed (left) and adaptive (right) temperature decay (initialized with q()) = q). Each 8x the median iteration of Alg. 1. Shaded areas show 80% confidence intervals around median over 100 random problems from the upsampled MNIST dataset (n 4096) with normalized L1 distance cost mmxx.j |Ca,yl 1).
In Fig. (left), we show that with sufficiently slow fixed schedule, óg typically remains near 1 throughout the execution of Alg. 1. Fig. (right) further confirms that the adaptive update rule arrives at such a schedule regardless of the initial value of g( thereby showing its urility in eliminating the need for tuning. Table verifies that the schedules found by the adaptive updating perform similarly to the best fixed schedule across OT problems. Section details the experimental setup here, and further adds benchmarking against a suite of alternative algorithms the literature.
3.4 ASYMMETRIC MARGINAL SMOOTHING FOR NUMERICAL STABILITY
Recall from (9) the form of the coefficient matrix (p) = D(r(P))(T pPre) of the
Wr
0.25
0.35
0.45
Total cost
n-variable PD linear system that are in-
terested in solving. Since the smallest diago- nal entry of (p) can be almost as small as
Median (ops.)
3431
2750
2795
90th %oile (ops.)
7622
3852
3869
(1 r)r (P) min in the worst case, diagonal preconditioning used Alg. 2 may cause nu- merical instabilities when solving the linear system (due to infinitesimal entries in r(P) for p s 1). To this end we require that, after smoothing, Fmin bounded away from zero, and find that running Alg. (F)) (i.e., the variance of F/r(P) 1) is sufficient for stable behavior in practice. Since log-domain Sinkhorn updates to ensure c P) = are numerically stable (Feydy, 2020), we allocate "smoothing budget" sa/2 mostly for the row-marginal Specifically, in L5 of Alg. 1,
Median (sec.)
1.96
1.56
1.55
90th %ile (sec.)
4.41
2.16
2.13
Table 2: Comparison of median and 90th percentile performance for varying smoothing weight wyr
advance to bound
(1 e. (wpSa/n)1,
è (1 ween (weaa/n)1,
where wy + W. 1/2 and w, > We. This is contrast the more standard symmetric smoothing Wp We (Dvurechensky et al., 2018; Lin et al., 2019). We repeat the experiments in Table 1, but this time ablating w, (using the adaptive schedule introduced in Sec. 3.3); Table 2 confirms empirically that asymmetric smoothing improves both stability and median performance of the overall algorithm.
4 EXPERIMENTS
In this section, we first provide the details of the MNIST experiments in Fig. and Tables 1-2. Then describe color transfer problem set used for wall-clock time benchmarking of the combined MDOTT-runcate/Newwtor algorithm. Our setup and baseline implementations follow Kemertas et al. (2025). All experiments
NVIDIA GeForce RTX 2080 Ti GPU with 64-bit precision. Appx. D provides benchmarking on 10 additional datasets from Schrieber et al. (2017), showing similar results with confidence intervals and including operation counts alongside wall-clock time.
4.1 EXPERIMENTAL SETUP
Upsampled MNIST. Each image is a probability distribution over pixels, with probabilities given by L normaaalzeed intensity values. The cost matrix C is fixed across OT problems and constructed from
8

---

Published as a conference paper at ICLR 2025
Subo gap wall-clock time (GPU)
upsampled MNIST. L cost
upsampled MNIST. L3 cost
10-1
10
10
10-2
10
10-1
NDOT'S
10-4
10
MDOTIPNOG
10 5
MBOT croeeckwee
10-5
Jimbrrsad es (2520)
10
Sikhan (pm
10
Sikhum (pm1
10
MSK
10
10
Peythy Map
APDAGD
10-4
0
10-9
color transfer, cost
color transfer, L3 cost
10
10-2
10-2
10-3
10
10
10a4
10
10-9
10-5
10-1
10a5
10-9
10
10-9
10-2
10-1
10 10a
10a
10-2
10-1
10
10a
10a
Wall-clock time (s)
Figure 2: Error vs. wall-clock time for various algorithms. Each marker shows the optimality gap and time taken (median across 18 problems) until termination at given hyperparameter setting, followed by rounding of the output onto (r,, c) via Alg. 2 of A Altschuler et al. (2017). Upsampled MNIST (top) and color transfer (hottom) problem sets (n 4096) using L1 (left) and 1.3 (right) distance costs. MDOT-TruncatedNewwe
outperforms others by orders of magnitude at high precision and exhibits much better practical dependence on error than best known theoretical rates Ö(n²==1).
pairwise L1 and L3 distances between pixel locations in 2D space. Scalar division by the max. entry ensures max., |Ciyl 1. MNIST images are upsampled to 64 x 64 resolution for benchmarking on higher dimensional problems (n = 4096). In contrast, Tang et al. (2024) benchmarked their CPU-based algorithm on original MNIST images (n = 784). Luo et al. (2023) ran their PDASMD algorithm on downsampled MNIST images for 100 at most; their code also runs on a CPU and includes an inner for loop of n iterations which limits parallelization.
Color Transfer. We define the color transfer problem with all marginals set to the uniform distri- bution over An. The cost matrix C' varies across problems and is constructed from pairwise L and I.3 distances between RGB values in 3D space. Scalar division by the max. entry ensures maxi Ciyl = 1. We use the 20 images provided by Kemertas et al. (2025), which were generated by prompting DALL-E 2 to produce vibrant, colorful images with intricate details. These are down- sampled to 64 x 64 resolution for benchmarking on = 4096 problems, except in the case of Fig 3, where the scalability of our algorithm is visually demonstrated the original 1024 x 1024 images (n ~ 10") on an individual sample problem.
4.2 WALL-CLOCK TIME BENCHMARKING
While theoretical analysis and computational complexity provide valuable insights, it is the practical performance, measured in wall-clock time, that often determines the viability and adoption of an algorithm. In this section, we present wall-clock time benchmarking of the proposed algorithm against a broad range of available alternatives. Here, we compare to Alg. 3.5 of Feydy (2020), the Mirror Sinkhorn (MSK) algorithm of Ballu & Berthet (2023), Sinkhorn iteration with typical and stringent tolerance settings (p (p=11 and p = 1.5 resp. for MDOT-Sinkhorn called with 7: = 7r), the APDAGD algorithm of Dvurechensky et al. (2018), the Mirror Prox Sherman Optimized algorithm of Jambulapati et al. (2019) and the AAM algorithm of Guminov et al. (2021). We omit comparison to APDAMD (Lin et al., 2019), APDRCD (Guo et al., 2020) and PDASMD (Luo et al., 2023), as they exhibit significantly longer convergence times in our high- -igeenimnsiinnl,, GPU-parallel setting. See Appx. E of Kemertas et al. (2025) for implementation details. In Appx. C, we also perform experiments with varying problem size n and show that the dependence is no worse than (nn)) for the problems considered.
9

---

Published as a conference paper at ICLR 2025
Image A
Image B
A B
B A
Figure 3: The Troe--TuunattedNeeonn algorithm applied to a large-scale color transfer problem on 1024 1024 images (n ) For this visualization, the cost matrix is given by the L distance in RGB color space, normalized so that max; |Ciy 1. Final temperature is 1/7r 2 2-10, Source images (top row) were generated with DALL-E 2. This figure best viewed digitally.
In Fig. 2, we observe that 2018 2-18-era GPU, MDOT--TruncatdNeeweer typically solves 4096 dimensional problems in 1-5 seconds to 6-decimal precision, demonstrating numerical stability up 9-decimal precision across a range of realistic OT problems. In the highest precision range, is more than 10x faster than the best alternative, MDOT-PNCG of Kemertas et al. (2025). We also implement a scalable (n) memory footprint version of TiDOT-TrnnateeNNeeto which computes via the PyKeOps package of Charlier et al. (2021) the cost matrix C' and the matrix P on-the-fly every-time they are used in Alg. 4. As evidenced by Fig. 3, this implementation can solve very high dimensional OT problems (n RS 10%) to high precision. It leaves a memory footprint of just se 600 MBs, but takes Re 10 hours. Regardless, we believe that this is an important step towards discrete OT in very high dimensions.
CONCLUSION
In this work, we set out to design a modular, practical algorithm to exploit the superlinear convergence of truncated Newton methods in weakly-regularizer EOT. To improve the conditioning of the dual Hessian, rather than amplifying its diagonal entries as in Tikhonov regularization, we dampened off-diagonal entries (discounting) with inspiration from reinforcement learning. Then, Alg. was presented for approximately solving the modified Newton system. This method of Hessian modification enabled superlinear local convergence rate in terms of the L1 norm of the gradient for the custom truncated Newton routine (Alg. 4) that used Alg. 2. We additionally introduced precautionary measures improve the numerical stability of Alg. 4, which is crucial for reaching high precision. Lastly, Alg. 4 was integrated into temperature annealing framework, MDOT (Kemertas et al., 2025), where adaptive temperature updates ensured superlinear convergence is maintained, hyperparameter (g in Alg. 1) is eliminated and line search overhead is minimized in practice. We implemented the resulting algorithm on a GPU and showed that it outperforms many recent algorithms by orders of magnitude in 4096 dimensions, exhibiting fast empirical rates ranging from e 1/6) to Dnn2-11/3) Furthermore, as visualized in Fig. 3, the algorithm holds potential to effortlessly scale to much larger problems (n 10%).
One avenue for future research is the development of variant of Alg. 4 with global convergence rate, which may yield an explicit rate in terms of the error for the overall algorithm. Another direction could be stochastic generalization of Alg. 4 that leaves O(nm) memory footprint, where m E [1, n] is a user-defined parameter given hardware constraints. This added flexibility may enable better of between runtime and memory footprint in very high dimensions (n >> 1000).
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
AMF acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) through the Discovery Grant program (2021-03701), MK acknowledges the support of NSERC through the CGS-D scholarship. Resources used in preparing this research were provided, in part, by the Province of Ontario, the Govermment of Canada through CIFAR, and companies sponsoring the Vector Institute.
REFERENCES
Jason Altschuler, Jonathan Wiles-Weed, and Philippe Rigollet. Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration. Advances in neural information processing systems, 30, 2017.
Owe Axelsson. Iteration number for the conjugate gradient method. Mathematics and Computers in Simulation, 61(3):421-433 2003. URL https: www. swwencceiiiiieeeeeoococ.eeeiieee rect, com/ science/ 7577744444000777777 6. MODELLING 2001 Second IMACS Conference Mathematical Modelling and Computational Methods in Mechanies, Physics, Biomechanics and Geodynamics.
Marin Ballu and Quentin Berthet. Mirror Sinkhorn: Fast online optimization on transport poly- topes. In Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learming Research, pp. 1595-1613. PMLR, 23-29 Jul 2023. URL
httpps: /ttps://prrceeeee ingeeediiig.... wwg..mm........................m ml p.reesppesss...... v202/ba 1u23a htmll.
Mathieu Blondel, Vivien Seguy, and Antoine Rolet. Smooth and sparse optimal transport, In Proceedings of the Twenty-Firsl International Conference on Artificial Intelligence and Statisticy, volume 84 of Proceedings of Machine Learning Research, pp. 880-889. PMLR, 09-11 Apr 2018. URL https: http///pooeeeeeeeeee/eeeeeee mlr. wwwpppressssseee http://4////ooo htt////ooeeeeeeee 1te.l11aa...t.....
O Bokanowski and B Grébert. Deformations of density functions in molecular quantum chemistry. Journal of Mathematical Physics, 37(4):1553-15733 1996.
Haim Brezis. Remarks on the Monge-Kantorovin problem in the discrete setting. Comptes Rendus. Mathématique, 356(2):207-213, 2018
Benjamin Charlier, Jean Feydy, Joan Alexis Glaunes, François-David Collin, and Ghislain Durif. Kemel operations on the GPU, with Autodiff, without memory overflows. The Journal of Machine Learning Research, (21):345773366 2021.
I. Csiszár and P.C. Shields. Information theory and statisties: tutorial. Foundations and Trends in Communications and Information Theory, 1(4):417-528, 2004. URL http: //x. doi. org/ 10.1561 01000000000000000000000000000000
Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems, 26, 2013.
Pavel Dvurechensky, Alexander Gasnikov, and Alexey Kroshnin. Computational optimal transport: Complexity by accelerated gradient descent better than by sinkhorn's 's algorithm. In International conference on machine learning, pp. 1367-1376. PMLR, 2018
Norm Ferns, Prakash Panangaden, and Doina Precup. Metries for finite Markov decision processes. In Proceedings of the 20th Conference on Uncertainty in Artificial Intelligence, UAI '04, pp. 169, Arlington, Virginia, USA. 2004. AUAI Press.
Jean Feydy. Geometric data analysis, beyond convolutions. PhD thesis, Université Paris-Saclay, 2020.
James Allen Fill. Eigenvalue bounds convergence to stationarity for nonreversible Markov chains, with an application to the exclusion process. The annals of applied probability, pp. 62-87, 1991.
Reeves Fletcher and Colin M Reeves. Function minimization by conjugate gradients. The computer journal, 7(2):149-15) 154, 1964.
11

---

Published as a conference paper at ICLR 2025
Joel Franklin and Jens Lorenz. On the scaling of multidimensional matrices. Linear Algebra and its applications, 114:717-735, 1989.
Aude Genevay, Gabriel Peyré, and Marco Cuturi. Learning generative models with Sinkhorn divergences. In International Conference on Artificial Intelligence and Statistics, pp. 1608-1617.
PMLR, 2018.
Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville. Improved training of Wasserstein GANs. Advances in neural information processing systems, 30,
2017.
Sergey Guminov, Pavel Dvurechensky, Nazarii Tupitsa, and Alexander Gasnikov. a combination of alternating minimization and Nesterov's momentum. In Proceedings of the 38th International Conference on Machine Learning. volume 139 of Proceedings of Machine Learning Research, pp. 3886-3898. PMLR, 18-24 Jul 2021. URL https: htttp:///erceeddngg..........s. http:nns.....rpes/s/sssss guminov2 21 html.
Wenshuo Guo, Nhat Ho, and Michael Jordan. Fast algorithms for computational optimal transport and Wasserstein barycenter. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pp. 2088-2097. PMLR, 26-28 Aug 2020. URL https: http:ttts///rocceeeeitiiie pttpp/eeeeenggn.oppeeepp.o ..s.s.m...... m www.r.presssvlec guo20a. html.
Arun Jambulapati, Aaron Sidford, and Kevin Tian. A direct ((1/e) iteration parallel algorithm for optimal transport. In Advances Neural Information Processing Systems, volume 32. Cur- ran Associates, Inc., 2019. URL https: proceedd ngs neurips cc/paper__ ffffifie///// paper/2019 9/ 9t/f///il//.. file/ 02 0222422699eenom 02444
Mete Kemertas, Allan D. Jepson, and Amir massoud Farahmand. Efficient and accurate optimal transport with mirror descent and conjugate gradients, 2025. URL https: /ttp:/:////xxvvo///. abs/2301 08507v3.
Philip A Knight. The Knoop algorithm: convergence and applications. SIAM Journal on Matrix Analysis and Applications, 30(1):261 275, 2008.
Yin Tat Lee and Aaron Sidford. Path finding methods for linear programming: Solving linear programs in (Vrank) iterations and faster algorithms for maximum flow. In 2014 IEEE 55th Annual Symposium on Foundations of Computer Science, pp. 424-433, 2014.
Christian Léonard. From the Schrödinger problem to the Manto-Kantooovic problem. Journal of Functional Analysis, (62(4):8879919) 2012.
Bruno Levy, Roya Mohayaee, and Sebastian Hausegger. A fast semidiscrete optimal transport algorithm for a unique reconstruction of the early Universe. Monthly Notices of the Royal Astronomical Society, 506(1):1155-1185, 06 2021. URL /tttps://doi..or/1....1. mnras/ sta61666...66
Tianyi Lin, Nhat Ho, and Michael Jordan. On efficient optimal transport: An analysis of greedy and accelerated mirror descent algorithms. In Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 3982- 3991. PMLR, 09-15 15 Jun 2019. URL https: http://prooeednin...........r.ee html.
mmr.
/ 9a.
Tianyi Lin, Nhat Ho, and Michael I. Jordan. On the efficiency of entropic regularized algorithms for optimal transport. Journal of Machine Learning Research, 23(137): 1-:37):::42,, 2022. URL http: // //mrr.. org/ pttpps/ppprs/333 20 277 html.
Yiling Luo, Yiling Xie, and Xiaoming Huo. Improved rate first order algorithms for entropic optimal transport. In Proceedings of The 26th International Conference Artificial Intelligence and Statistics, volume 206 Proceedings of Machine Learning Research, pp. 25-27 Apr 2023. URL ht / proceddi ngs ml www.rrrrrrss.....mmmmmmmmmmmmm.e presss v20 uo2. 3 html.
PMLR,
hnttpss:/proceeeiiii
12

---

Published as a conference paper at ICLR 2025
Quentin Mérigot. A multiscale approach optimal transport. Computer Graphics Forum, 30(5): 1583-1592 2011. URL https: /oniieeibbaary. wiley com/ doi/ abaa /00. j.1467-8659. 2011. 02032.X.
Stephen G. Nash. survey of truncated-Newton methods. Journal of Computational and Applied Mathematics, 124(1):45-59, 2000. URL https: /www. wwww.scci.. www.sciienccc.e rect com/ mennsciincce Numerical Analysis 2000. Vol. IV: Optimization and
article/pile httpp///////300000000/0//1 Nonlinear Equations.
Jorge Nocedal and Stephen J. Wright. Numerical Optimization. Springer, New York, NY, USA, 2e edition, 2006.
Marcel Nutz. Introduction to entropic optimal transport. Lecture notes, Columbia University, 2021.
Ofir Pele and Michael Werman, Fast and robust earth mover's distances. In 2009 IEEE 12th International Conference on Computer Vision, pp. 460-467, 2009.
Pitie, A.C. Kokaram, and R. Dahyot. N-dimensional probability density function transfer and its application to color transfer. In Tenth IEEE International Conference Computer Vision (ICCV'05) Volume volume 2, pp. 1434-1439 Vol. 2005.
Igal Sason and Sergio Verdú. f -divergence inequalities. IEEE Transactions on Information Theory, 62(11):5973 2016.
Bernhard Schmitzer. Stabilized sparse scaling algorithms for entropy regularized transport problems. SIAM Journal on Scientific Computing, 41(3):A143-AAI: 2019. URL hwtps::t//://tt //oi. org/
10. 10.1137/ 6M1106018.
Jörn Schrieber, Dominic Schuhmacher, and Carsten Gottschlich. Dotmark a benchmark for discrete optimal transport. IEEE Access, 5:271-282, 2017.
Jonathan R Shewchuk. An introduction the conjugate gradient method without the agonizing pain. Technical report, Carnegie Mellon University, USA, 1994.
Xun Tang, Michael Shavlovsky, Holakou Rahmanian, Elisa Tardini. Kiran Koshy Thekumparampil, Tesi Xiao, and Lexing Ying. Accelerating Sinkhorn algorithm with sparse Newton iterations. In The Twelfth International Conference on Learning Representations, 2024. URL https: / //oeerreview net. forum? 5dg
Yujia Xie, Xiangfeng Wang, Ruijia Wang, and Hongyuan Zha. fast proximal point method for computing exact Wasserstein distance. In Proceedings of The 35th Uncertainty in Artificial Intelligence Conference, volume 115 of Proceedings of Machine Learning Research, pp. 433- 453. PMLR, 22-25 Jul 2020. URL https: /ttp//eeeeenns...//ee/p/oo ml ptress/v11///iie html.
13

---

Published as a conference paper at ICLR 2025
A APPENDIX
A.1 PROOFS OF THEORETICAL RESULTS
A.1.1 ADDITIONAL NOTATION
Here, we describe additional notation used in the proofs that follow. Given a column vector a E IR" and an n X n square matrix A, || 2 = x
_ Ax. The operator norm of a matrix is denoted ||Al| suP=ER" ||A | To mean [|Allp.p we use the notation ||All, For an element-wise norm, we write ||vec[A)]] where vec(A) is the vectorization of the matrix A. vector formed by diagonal entries of a matrix P is denoted diag(P).
1.1 DERIVATION OF THE EOT DUAL (4)
In this section, we provide the derivations for (3-4), Recall the EOT primal problem given by (2):
(P.C) - _ H(P. 7
minimize P (
Observe that we can replace negative Shannon entropy (P, log P) in (2) with the KL divergence (written without assuming P E Anx)) to the uniform distribution U = Inxi
1 1) +
(KU(P (P, log P) (P,
Σ Pij
n
log P)
Prs
Pij + 2logn++ +1
That is. replacing the negative Shannon entropy term with the above in the primal objective only increases the objective by a constant 2 log n on the feasible set, since Eij P. = =1for all U(r, c). For convenience, we drop the 2 log n term and take the primal problem given by
minimize P u(r,c)
(P,C) +
(P.log P) - EP +1
of
To derive the dual problem, we write the (scaled) Lagrangian with dual variables u, E R corre- sponding to equality constraints r(P) = r and c(P) = C:
C(P, u. v) = x(P,)) (P, log P\ Σ Pa +1+ (u, r - r(P)) + (v, (17)
Taking the first derivative with respect to the isth entry P.j of P:
"Cg + 1 log Psj - 1 ut U/j
OP.
= -Cis + log Pij Fj.
Setting the partial to O:
0C aP,
0 + log P.j = u, + 03 - CCg
P. explu: + 03 70
where the last equation is the same (3). Now, plugging the above P into the Lagrangian (17):
-g(u, v) = (P, (P, + To -0) Σ Pij + 1+ (u, r(P)) + (0.c c(P))
= -u, + (0, Σ P.y (u, r(P)) - c(P))
ij
=1 Σ
P(u, U)1s + (u, r) + c).
17
14

---

Published as a conference paper at ICLR 2025
Maximizing -q(u, 0) with respect to the dual variables, obtain the dual problem (4):
minimize u U E IR"
g(u, v) - P(u, - (u,, r) (v, c)2
where we kept the constant -1 convention.
As an aside, note that if one assumes instead an L1 normalized form for P via the softmax function (given that the feasible set (r, c) is a subset of the simplex):
expf4 + U1 -C
P(u, v)s
Σx, exp(ux k + uz GGg)
one obtains an alternative form for the dual objective by plugging the above into the Lagrangian (17):
g(u, v) = log Σ explu: V/ >Ci) (u, r) (v,c). c).
Both the -f-oo-e--xponett and log-sum se--sum-off-xooentttt forms the dual appear in the literature; see for instance Altschuler et al. (2017); Dvurechensky et al. (2018); Lin et al. (2019) for the former and Lin et al. (2022) for the latter. Both objectives g and have the same value whenever P(u, v) as defined in (3) is on the simplex, since = log x = 0 for x = which is why the constant -1 in (4) was kept as a convention.
A.1.3 PROOF OF THM. 3.1
We start this section with intuitive example describing the role of the row-stochastic matrix Pre that was defined in (9). The discussion carries over to Per by symmetry. Next, we will list some mathematical properties that will be useful in the proofs that follow.
Suppose r and are disjointly supported on two sets of particles I1, Ln, and Un respectively, and let n2 tributions over and y particles. Recall the definition of Pre P. in (9) as product of P. D(r) x and P. = D(c)-IpT E IR Given some initial distribution E A., over --particlss, P. q E An, is distribution over y-particles after transportation accord- ing to P. Indeed, we can show casily that P q is the simplex:
n. That is, an n1 n2 transport plan P u(r, e) maps dis-
= P
D(r) - -1-r))---1-------
Similarly, P. PT - PT An, is again a distribution over x- x-particles after transporting back according P. This is the stochastic process represented by Pre. Given any initial distribution q. the process converges to the row-marginal r of if Pre. is applied repeatedly (as the next lemma shows). The second largest eigenvalue A2 of Pre determines how quickly this convergence occurs (see Lemma A.3 below).
Now, we list useful technical properties of Pre our setting, Analogous claims hold for the stochastic matrix Per and the column sum c(P) by symmetry, but are omitted for brevity.
Lemma A.1 (Properties of Prc). Given a matrix E Rxx with strictly positive finite entries define The following are true:
Pre P.P.
1. Pre is an irreducible row-stochastic matrix. Its second largest eigenvalue x2 is strictly less than one.
2. The stationary distribution of Pre is (PP).
3. Pre is reversible in the sense that D(r ((r(P))Pre Prr)), which implies all eigen- values of Pre are real.
4. PreD(r(P)) D(rr(P))-1
5. Given some p E [0, 1), we have (1 = -
Proof. We prove each claim in order.
15

---

Published as a conference paper at ICLR 2025
1. The vector of ones is right-ciggenvector of Pre with eigenvalue 1:
Pre1 -
Drr(P))=rPD(((P))L((( D(r
-rr((
= 1.
Since all entries of P are strictly positive, the same true of Pre. The claim follows.
2. The vector r(P) is a left-eigenvector of Pre with eigenvalue 1:
r(P) Tre -r(P) -r(()) DD((P))
=1 -1P ((c()) PT = =PP
=r(P)T.
3. The claim holds since
D(r(P))Pre
D PD(c(P))-\p'
PD(GP))-PPPD(())))-D(1>> Drr(()
4. The claim follows similarly the previous claim.
5. First, notice that Claims 3 and apply analogously to all powers of Pl for I 0. Indeed, for Claim 4 we have:
pl (r(P)) = ProDP))
PC ((r)D(r(P)) -ipT
(by Claim 4)
-(r-)),D(r((P)) T
(H=)
(again, by Claim 4)
-P(r(P))-
(by repeated application of Claim 4)
Then, from the Neumann series (I - (Pr
pl D(r()))
we obtain the claim by applying the above equality to each element of the sum.
Lemma A.2 (Properties of the coefficient matrix Fr(p)). Let E [0,1), Pre be row-stochastic matrix with strictly positive entries and min; r r(P): 2 for all [n]. Given an n x n matrix F,(p) D(r(())((T -Ppe)) pPre), the following hold true:
Fx(p) is symmetric, positive-definite matrix.
Amaze (F+(P) S (1 + p)rppr((pp)) max
Amin (Fr(p)) 2 (1 (1-p)r(P)) (P) min.
Proof. Observe that
Ft(p) = ((r(())([ pPre) = D(r(P)) pPPc D(r(P)) PPPD(((P))-pp
a sum of two symmetric matrices, which is also symmetric. The eigenvalue bounds follow from the Gerschgorin Circle Theorem. Since the smallest eigenvalue is positive, Fr(p) is dositive-definite
16

---

Published as a conference paper at ICLR 2025
Before move to the proof of Thm. 3.1, state the following Lemma, which follows immediately from Thm. 2.7 Fill (1991) given that Pre reversible by Claim of Lemma A.I. Lemma A.3 (Convergence the stationary distribution under Pre). Given some An.
(P ( 5 (rrr()))
(18)
where X2 is the second largest eigenvalue X2 of Pre.
Now, we provide more formal version of Thm. 3.1 presented in the main text followed by a proof.
Theorem 3.1 (continuing from p. 5). Assuming c= c(P) and dy = =Pddd define residuals eu (p) + u. and = V d + Vg. Further, let X2 E (0,1) be the 2nd largest eigen- value of Pre and Ç = Vuglll (x(rl(())) I 1. For any 3 E (0, 1), suppose:
0, - (1- K2)
max
s
(19)
A2(1 K)
3
leu(p)l1 - 2 Vgll1
(20)
where = 5B8n < 1. Then.
|ell |l,, (1)11 S Vgll1
(21)
Proof. To establish necessary conditions for bounding ||e|| first, we write it out explicitly:
V'gd + Vg
D(((P)) P D(c(P)))
dhi d.)
Vuy Veg
pJ
+
D(r(P))
du
V
Vug
PT
D(c(P)))
PPdd
+
(by construction)
D(r(P)) - PPodu
Vug
(since D(c(P)) PT by definition.)
0
(I Pre)du
Vu9
(22)
0
0
ea (1)
which proves the equality on the LHS of (21), where the last equality holds given the definitions (u())==ff(ppd, + Vug and (y(p) = =D(r(P))(1 pPre). Hence, bounding [[)]())]|, suffices. Now, given the definition of e (p) and the invertibility of Fr(p) for p < 1 by Lemma A.2, have:
d, = Fp(p) ((((()))(e()))) -
Plugging into the top half of (22), observe that:
en (1) F,(1)PF(0) - Vng) + Vny -(p)eu(() + (T- (())
where we defined ?(p)= (.(1)F Then, have
|(oo|l ||en(p)ll + v())
[ea <
(23)
[(o)
First, we prove that the operator norm
2.
Î(p) = D(r(P))(( Pre)(1 pPre) -P D((((P))
-
= (r1p)))((PP))PPpp) Prc(T -ppT--pPree-tt -Dr
I (1-11)))pp pPTT)=
(24)
17

---

Published as a conference paper at ICLR 2025
where the last equality follows from claims and 5 of Lemma A.I. Recalling that ||A||] is the maximum absolute column sum of matrix A:
T(p) +1+ p) Prell ||(T ppre) -11
=1+(1-p)
(1 pPre) -1
=1+(1-p) p)
Σ p (PT))
$1+ P) Σ
(PTx)))
1=0
2.
(since P is stochastic matrix for all / / 0)
Hence, (23) simplifies to
(eell) 2 + 111 Vo)
(25)
Next, we turn to the second term on the RHS:
(T - Vug = (1 p) p) - -pPre) pPpo)-1 -'V
(From (24))
(1 p) Σ
(P 11 ( (+(() -()
(1- Σ
(((P) (r(P)- (pTo)
1=0
where the last equality is due to the fact that r (P) is the stationary distribution of Pre by Claim 2 of
Lemma A.I. Then,
||-- ((- us/ 5(1-p)
Σ p
r(P)) (1P)-(PPT) -(PT))+1))
1=0
< (1 - p)
Σ
(((P) - (1tT))t17))
1=0
S - p) Σ
x N x(rlr(P))
(By Lemma A.3)
1=0
(1 p)>₂
1 - pb2 x(rlr(P))
Plugging this bound back into (25) and continuing with the main conditions given in the theorem:
|ee(11)] 21(eem(p)1l + (1 p)>r
x(r(r(P))
$ (1 B)n |gll + (1 -1- p)12 x(rlr(P))
(Since (20) holds by construction)
- pdz
$ (1 B)m Vgll, + 1|Vll
(Since (19) holds by construction)
which concludes the proof. Above, the parameter E (0, controls a trade-off between how precisely the discounted system is solved and how aggressively the original system is discounted. In Algorithm 2 of the main text, 3 was fixed at 1/2 for simplicity. Also, for intuition on the effect Xg. observe that the second term vanishes X - 0 (if Pre mixes quickly) so that the Hessian can be discounted more aggressively with smaller As we see next, this improves our guarantees on the condition number of the linear system.
18

---

Published as a conference paper at ICLR 2025
A.1.4 PROOF OF THM. 3.2
Lemma A.4 (Spectrum after preconditioning). Let (p) preconditioned coefficient matrix, where M ((d)ia(FF()))), Further, let μ -diag(P Pra) be the
(- 1/2 be the diagonally
diagonal entries of the stochastic matrix Pre E 1R Then, eigenvalues (p(p) of (p( satisfy
p(1 Mminn
<x((p) 11+ p(1 Mmmin)
Vi [n]
(26)
PPumin
P/mmiu
k(p) Amax (p)
Mmiu
2 Ummi +
I 1- P
(27)
Amia (p)
1
Proof. The proof of (26) follows straightforwardly from the Gerschgorin Circle Theorem. Due to diagonal similarity, the spectrum of (p) coincides with that of Fr(p) F, Clearly, diagonal entries of (p) all equal to 1, so that all Gerschgorin disks are centered around unity. Then, all eigenvalues must be inside the biggest disk, which contains all of the smaller disks. First, consider row i (p) = D( pPre):
Σ - r(P):p(11 -1). iti
Since Ma Fr(p)u = r(P).(1 p)4), we then have:
(r(c)) p(1 m4)
Σ AZ
(1- pHL)
where the biggest Gerschgorin disk corresponds to Pmin. so that (26) holds for all E [n)]. Then,
2
<(p)= S Amin (p)
I
1
p(1- Minin
(I Opmnm
Mmin
Pmnin
1 p
1 Mmmia
0
as p 1
- P
Lemma A.5 (Equivalence of norms). Suppose d', satisfies (p)dd = -Vg. Let Fx(p) = M (/) 1/2 and M = ((diag(F.(p)) as in Lemma A.4. Define the reparametrization d - d, given some d E R", and the residual eu = (p) + Vug. We have,
leull
V (1 - p)r(P) S
$ (11 (1+1 pmn.
(28)
du
dig
P-(0)
19

---

Published as a conference paper at ICLR 2025
Proof. For this proof, drop p from Fx(p) and (p) for convenience. First, observe that du d = Fr Cu. Then,
du dzl 2 P.
-|| F eu P.
=
-e M 1/2 2 eu
e e.lt F.FT Eu
eu eu
Amin (Fr-
| l
| 1
Amax Fr)
| 2
(by Lemma A.2)
(1 + p)r
)maa
l 2
(1 + p)
||eu ll
(since ||>|| S vn |/|12|12 Vx E R".)
n21 + p)
which is equivalent to the upper bound desired result. The lower bound follows similarly the reverse direction, using Amin (F.) and ||a:|| I ||=]]1-
Lemma A.6 (Convergence of CG). Suppose diagonally-preconniiione conjugate gradient method
initialized with ! 0 for the linear system F, (p)d -V where p 10, Let μ be the largest diagonal entry of Pre- Assuming that r(P) min 2 ea/(4n) for all E [n] given some
constant ed O, the residual satisfies [[F+(p) Vugl, = [()l())]l S V after at most ceil( (k) steps, where
k 5
1 up
log (6n(1 - pp)-1/21 -/
_222 = )/2)==0 (0(1
T- ==
(1T=
as p 1.
(29)
1
1-p)
Proof. Once again, we drop p from (p) and (p) for convenience. Using the same definitions as Lemma A.5, recall the equivalence of the aaggonally-precondiii linear system:
F,d, -
M 1-1/2 --M-1/22 n
M 1/ 1/2 1/ 1/2 M 1/2 du = - - ug
F, -
To guarantee [|eull n V Lgll1, is sufficient to have (given the upper bound in Lemma A.5)
d., a. R < n1/Vugll 22. Vn(1 + p)
(30)
It is well-known the that linear conjugate gradient method ensures l dull p.00) most ck) steps (Shewchuk, 1994), where
< 2s after
1 k: == 2 Vm((, ((F,)) log
-l
(31)
2
20

---

Published as a conference paper at ICLR 2025
Then, have
k 3 2 V x((F))) log
1
d e-
Ar
1 V m(FF) log
Eu
E
$ 2
(by the lower bound in Lemma A.5)
V(1 p)r(P min
1 m(FF)) log
Vuglli E-1
(since d = - eu = | (u9ll
2
p)r(P))
1
(1+ p)n
V ((Fr) log 2n² V (1 (p)rr())mn
(using from (30))
2
1
$ V x(FF)) log 3n 112 1)-1/211 1r(mmnn
(simplifying constants)
$ 1
k x(Fr) log 6n p)-
1122 -- 1/2
(since r(P)min 2 ea/ by assumption)
up
log 6m - p) p) 1/2 -1/2 1 3 11/2
$ 2(1 p)
(by Lemma A.4)
which concludes the proof.
Theorem 3.2 (Convergence of Algorithm 2). Suppose r(P) =a/(4n) given ea > 0 and each =0. Alg. 2 terminates in step Alg. 2 runs dlagonally-prrcondiitonn CG initialized with "
(1 ()((r(((P)) δ n V (1 X2) - rll
log ed
1/2
(14)
operations, where d < is the 2nd largest eigenvalue of Pre and μ < 1 its smallest diagonal entry.
Proof. An easy way to why Thm. 3.2 holds is by noticing that the lower bound condition on p given in (19) will hold after finite number of iterations depending logarithmically on remaining problem parameters. Then, since by Lemma A.6 each iteration requires O(11 p)-1/2) steps of CG (with each step costing O(n))), and - P will be at most O >22), the result follows. Now, we provide more detailed analysis.
First, define ( = 1 p(0) for the 4th iteration of Alg. 2. Since p(0) = 0 - f(o) 1, we have T = 4-1 given the update rule in L5. Then, by (19) of Thm. 3.1, Alg. I performs the final iteration (in the worst case) when
34--4's 1->2)K (1 A2) A2(1 K)
That is, in the worst case we terminate after / steps, where
log(1 - X2) log(4 4)
= ceil
(32)
In the worst case, we are guaranteed by the final step final that,
( = 4' (n- (1-A2)) 22
4 ,(1 -
(33)
Thus, for each integer I' [4], have
-4" (1
-44-11
(34)
21

---

Published as a conference paper at ICLR 2025
Then, the condition number at given step l' 2 0 satisfies:
1 1
S μr +
(From Lemma A.4)
1 p(0) 1-p(3
S 1z + 4 42+1-11 1 1
(From (34))
= =4 44' -IRos
where
1 H
Ko 0
as 22
(35)
(1 X2)
By Lemma A.6, given step 1' < of Alg. 2 takes at most
k((1) = 21+ -1 Vvro log 3n -0 6d 1/2
< 2 /ro log (nn/kon 1 sd
-1/2
(since 1' 1)
=2 -1((-))((VRO)))
steps of the conjugate gradient algorithm. Taking sum over l':
Σ k(() - ((v/5)))))))) )2 2
P=0
P=0
(2 -
=
(36)
2-1
= - - 2
S <Ō(VKo) V(1 (1-
(due to (32))
Ō(VKo).
(37)
Note that the choice of decay factor 4 is not arbitrary; if one chooses instead factor L 1 and repeats the above steps, one arrives at (cf. (32)):
log(1 - 12)
= ceil
log(L)
and the following (cf. (36))
VI i+1
Σ = =(V=))VV 1-1
VL-1
P=0
(1 1 k
s <ō(vKu)
VT. -1
< <((vm)) VL-1'
and = the global minimizer of h(x)) = x/(Vx - 1). The last inequality is strong when
convergence is bottlenecked by X2 R 1, so that the choice L 4 becomes near-optimal when
optimality most needed.
Continuing from (37), the overall complexity is:
Ō(Vio)) = (V
-
(1->2) (1
22

---

Published as a conference paper at ICLR 2025
The result follows by explicitly writing and simplifying K
k- A2 - K)
K
212(1 - -mm/2)
(Since we choose 3 = 0.5, where = (,31)
5n
2
$ {n
((r)r(P))
2
Vuugl
where the last equality holds by definition of Ç.
A.1.5 PROOFS OF LEMMA 3.3 AND COROLLARY 3.4
Lemma 3.3 (Convergence of Algorithm 3). Assuming that ||r(P)/rll < 00 and ||/rr() 00, Algorithm 3 converges O(n² /Ex) operations.
Proof. First, recall that by Eq. (169) of Sason & Verdú (2016), the ratio (r Pl) is bounded both x (P(())
above and below by constants depending on (r(P)/rllen < and (r/r(P)]l < oc. Notably, the
ratio converges to 1/2 as r(P) r as shown in Thm. 4.1 of (Csiszár & Shields, 2004).
Since we know that after ls steps Sinkhorn iteration,
(rl(((((P)) S A/k
for some constant A >>> (see Corollary 6.12 of Nutz (2021)) and x²(r|r(P)) ((Dri(()r(P))) we conclude that after & steps, x2 (2(r((())) = O(k-1). The result follows as each Sinkhorn iteration costs Onn²) and Alg. 2 terminates when x² (rlr(P)) Ex
Corollary 3.4 (Per-step Cost of Algorithm 4). If the backtracking line search in Alg. 4 converges in S iterations, then iteration of Alg. 4 costs (n + e 22/ x22 -1- operations, where X2 < is the 2nd largest eigenvalue of Pre defined as in (9) and evaluated at u, (ef. (3)).
Proof. First, from Lemma 3.3, the ChiSinkhorn routine in L6 of Alg. 4 has complexity 0 2//5 as we chose Ex = // Lines L9 and L18 each cost ( which leaves the cost of line search between lines 1.12-14 and the NewtonSolve routine L8. Since the former assumed to take steps, each costing (nn²), it remains to show the cost of NewtonSolve (or Alg. 2).
First, consider the case when 1 0.4ca/ |/ wgll as per L7. From Thm. 3.2, after dropping the logarithmic terms and the linear n² term, we have the following linear term:
(1 - (1-))x(r((PP))
V (1-M)x(r((PP)) (1 (1- Xe)Ea
(1
1-m 1/5
V (PP)
= 1 (1-A2)
-
2/5 222
since (2(r|r(P) < < 2/5 Next, consider the case when Vuoll 0.48a/ V 11 i.e., ||Vulll > 0.4Ea. so that L7 assigns 17 = | instead. Inserting n into the denominator on the LHS above and applying ? >> 0.4ca yields the same complexity.
A. PROOF OF THM. 3.5
First, we present simple lemma on the line search condition used in L12 of Algorithm 4. Lemma A.7 (Sufficient Decrease Condition). Assuming c c= given a descent direction d (du, -Pdd the Armijo condition for line search over objective g is satisfied for a step size a E (0, 1] and constant parameter e E (0, 1) if and only if:
vecc (P(u + ad a aPPd)))ll - 1 < c))a(- -Vug,
(38)
23

---

Published as a conference paper at ICLR 2025
Proof. Recall that the Armijo condition requires (see Ch. 3.1 of Nocedal & Wright (2006)):
g(u + adu; all - g(u, v) + cia d).
(39)
For brevity, let P(a) P(u adu; a.Pdu). We show the equivalence of the above statement to (38) step by step:
glu + adu, u a Pcdu) S g(u, v) +
l|vec(P(a))ll
(u + od r) (v aPedu, c) S (u, (0, + \vee(P(a))] 1 - a((dur r) (Podu, cga(V
|vve(P(a)) - 1 v a((du r) (du; Po c)) + (caa
|vee(P(a))]li - 1 I a(du. r) (d r (P))) + Vaa d)
(Since P.T c = (occr(p)) r(P))
||vec(((a))ll - 1 (d Vug) + tqa du), (Since c(P) = c by assumption)
which is equivalent to (38).
Theorem 3.5 (Per-step Improvement of Algorithm 4). Given descent direction d = (du: - Podu) such that ||ell, = V g + Vg S gullp ler a E (0, 1] be the step size found via backtrack- ing line search the hth step of Alg. 4. Then, Vgk+1 = Vg((u + adu aaPodw) satisfies
W < (1 + an)) 1 + av O(II ya 111)
(15)
Proof. For convenience, let = (u, v) given some (u, v) E R and (29(t)= = V2g(z + tad).
Recall that from Taylor's Theorem, we have:
Vgl=+ = Vg(z) + / (t))))
Vg(2)) + +V 1" [v'g(1) vig(0) d dl
- (1-a)Vg(z) + ae 1. [V'g(t) 4Z v²g(0)] d dt.
Now, we define the following:
((t) = (/' - -V)(()))))(0)) d.
Then.
Vg( + ad)l (1 (1--a) [gg / L/ n(E)at
$ (1 a + an) [g(2))1 + 1/// n(t)at]
(By construction)
L [lb(t)]l2 dt
< (1-a+on (V
(40)
We will bound [| in terms of tE [0, 1] and evaluate the integral. Again, define P(t) = P(z tad) for convenience, and let P = P(t) - P, where P = P(O). Then,
h(t)= D(r(AP(()) AP( D(c(AP(t))) ) (- du -Poduu (( (AP((t) AP(() dm,
AP(()
(
D(cc(AP(t))
(1((t)
(hs(t)
(2()))
(t((t)
24

---

Published as a conference paper at ICLR 2025
where we have [|t(t)]] 2 Ei=1 [11((t)l2 Consider (t),
| (1th = E (r(AP()):(d)))
- IIdulloo E ((AP(())
= |dull> |r(A((t)]]
- |dullo vcc(AP(t))ll
Since we have || loc - ||dull bo. the same bound holds for (122()) by symmetry. Next,
||t((t]lı ||AP(t) dull2
Σ (AP(t),; du>l
(A,g denotes the j column of A.)
- ||dullo 2 [1AP(t);ll
||du Do |vve(())(t))l)
Similarly, the same bound holds for [1aa(t) by symmetry. Then,
(()]) S 4 || ||vec(AP(()))>>
(41)
From Pinsker's inequality, we have:
2 vec(AP(())1 < (A(P))(()
(Bregman divergence under negative entropy)
(ivee((())))))) 1 + (P/PP(())P(() (Given c c= c((), E An..)
(vee((()))l) 1 to ((r( du) (c. (c.-Po
1 - to ((r(P), du) (Po c, duu))
(vec(P(t))ll: -1
(Since P. c= = PD(c) PP=
0.996o(-V dl,,;
where the last inequality is due to (38) as we assumed satisfies the Armijo condition, and if step size a feasible, so any step size to E [0, a] given that the objective is convex. Plugging the above into (41):
[[((t)l1 S 4 44|duuloo v2 0.99ta( (ee9tta(VVug-ae du)))
< 44 | 3/2 |V 1/2 V7.
Hence,
3/2 | L [ln(t)] di < 4v2a
/uolll/2 S( Vidt
3/2 2a
|d
|u
Val
[|d. 8//2
| 1/2
Vug)ll |V //nglll/2
S va12 (p)( 3/2 |uull?
(Since Alg. 2 ensures ||eu(p)]l 5
V.gll
VaO((
Plugging the above into (40) yields the desired result (15).
25

---

Published as a conference paper at ICLR 2025
Cost
Wall-clock time (s)
Adaptive
Po
Median
2.09
6.26
11 dist.
90th %ile
3.12
12.07
10th %ile
1.48
4.70
Median
5,49
25.74
L dist.
90th %ile
10.18
48.27
10th %ile
1,40
4.56
Table 3: Comparison of median, 90th and 10th percentile performance for adaptively and naively initialized p over 30 problems from the upsampled MNIST dataset with L1 and L3 costs 2 7 = 218 ,p
(n = 24,
pp= = 1.5, and g(1) = initially).
B ADAPTIVE INITIALIZATION OF THE DISCOUNT FACTOR
Recall that in Alg. 2, we initialize the discount factor at 0 and anneal (1 -P-) p) by taking:
- 11
in 14 of the algorithm until the forcing inequality (6) is satisfied. Here, we describe a simple, practical strategy reduce the overhead associated with this annealing procedure and the solving of sequence of linear systems (see also the proof of Thm. 3.2 in Appx. A.1.4).
In particular, we initialize NewtonSolve (Alg. 2) with initial guess po (rather than 0) in practice. Each call to NewtonSolve returns the final discount factor found by the algorithm in addition to the descent direction du. Then, the next time NewtonSolve is called, we call it with
Po
max (0, (1 (1-pold) 4),
(42)
where pood is the discount factor returned by the previous NewtonSolve call. That is, the annealing starts from the second last annealing step of the previous call. As the linear system has changed since the previous call, this allows for a smaller discount factor to potentially replace the previously feasible one, if appropriate. We find that this simple change in the implementation improves performance empirically, shown in Table 3. This version of the algorithm is used in the main experiments presented in Sec.
C PROBLEM SIZE EXPERIMENTS
Here, we conduct experiments with varying problem size n to empirically study the dependence of MiOOr-TruncateeNeewtor on n. Fig. 4 shows the behavior over MNIST and color transfer problems with L1 and L3 cost functions, and problem size adjusted by down- or up-sampling the images. In all experiments here, fix regularization weight at "r 212 In addition empirical behavior of MroT--TunnttedNeetto.. we include a polynomial (n) = an passing through the empirical curve at the largest n:; the curve explains the behavior of the algorithm well for large performs worse than (n²) empirically for the problems considered in Fig. 2.
D ADDITIONAL BENCHMARKING ON DOTMARK
In this section, we extend the study in Fig. 2 with 10 more datasets from the DOTmark benchmark introduced by Schrieber al. (2017) for benchmarking of discrete OT solvers. Schrieber al. (2017) proposed 10 different image sets, "to represent wide range of theoretically different structures, while incorporating typical images that are used in praxis and/or have been used for previous performance tests in the literature". Example image sets include various kinds of randomly generated images, classical test images and real data from microscopy; each dataset consists of 10 grayscale images, yielding total of 45 discrete OT problems, where the marginals r, are formed based on pixel values (Schrieber et al., 2017). The cost matrix is constructed similarly to the MNIST dataset from distances in 2D pixel locations. While Schrieber et al. (2017) proposed only evaluate on both L, and L costs functions for consistency with Fig. and for the sake of broader evaluation. Once again, for consistency with Fig. 2, we take 64 x 64 images, which yield 4096.
the L3 cost function,
26

---

Published as a conference paper at ICLR 2025
Wal-clock time vs problem sizw (upsampled MNIST, cost)
Wbll clock time VS. problem Sic (upsomp led MNIST. LJ cost ciCc
The Seoons
10-1
Wallelock time problem sizw (culor transfer, La cosl)
clack problem (color time
Wolkclok tine
30
Tino 5000E5
309
200
10
32-
207
307
Pribi 500
Predeen
Figure 4: Log-log plot of wall-clock time for MDOO--TrucateeNewtt vs. problem size Each marker shows the median over 60 random problems from the MNIST (top) and color transfer (bottom) problem sets with normalized L1 (left) and L3 (right) distance costs. Error bars show 10th and 90th percentiles, For all problems, = 2 Dashed lines show a polynomial f(n) = an², where is selected so that an equals the median time taken the largest considered. Above, the algorithm behaves no worse than (n²²).
For each of 20 problem sets (corresponding to class images and a cost function), we sample 20 random problems out of the 45 possible problems. Figs. 5-14 show the median time converge for each algorithm at a given hyperparameter setting, and the error (P P C) after rounding the output of the algorithm onto U(r c) with the exception of Alg. 3.5 of Feydy (2020); see Appx. of Kemertas et al. (2025). The wall-clock time plots for the respective cost functions (L1 and L2) follow the same trends seen in the two datasets considered Fig. 2. Following Kemertas et al. (2025), we include 75% confidence intervals along both axes here, and also show that MDOT- TruncatedNewton generally robust even at high precision, where maintaining numerical stability can be more challenging. Our conclusions based on Fig. 2 regarding the wall-clock time convergence behavior of MDOT-Trunce HIDOT-TruncatedNeettr and how compares baselines remain unchanged.
Operation Counts. In addition to wall-clock time, count here the number of primitive operations costing (n²) for each algorithm. Examples of such primitive operations involving n X n matrices include row/column sums of matrices, matrix-vector products, element-wise multiplication of matri- http:///wwwuuuctttoootttppppiiia
ces, element- wise exponentiation/logaritt between matrix and a scalar, over all entries of matrix, summation over all entries of matrix, etc. We count the number of primitive operations rather than a higher level function call such as the number of gradient evaluations due to inherent differences in the design of the various baseline algorithms; c.g., some methods require costly line search or inner loops between gradient evaluations. For all 20 problem sets displayed in Figs. 5-14, we find that the total number of O(n²) operations predict wall-clock time very well (high correlation), especially when the algorithms are run for long enough, as seen visually upon comparing top and bottom rows of the same column. All algorithms follow the same trend seen in Fig. 2, so that our conclusions once again remain the same.
matrices,
27

---

Published as a conference paper at ICLR 2025
Densityensity
L1
10
cost
10
VIOT Fnnttns
10-4
10
Sidbher
Optimality
10-1
10-2
L. cost
10
10
10
107
167
13
13a
10"1
20-1
221
107
O(nn) operations
Wall-clock time (s)
Figure 5: CauchyDensity problem with L1 (top) and L3 (bottom) costs, showing optimality gap vs. number of (nn²) operations (left) and wall-clock time (right).
Classiclmages
16-1
10-3
Li
10:4
cost
acal (2019)
16
Srknam (p=1.0)
Seknom 1p-5..5
16-9
MOK
APDAGO
Fodd Al.. 355
Optimality (P-
16-4
10-1
12-4
cost
10-
10-2
12.
1.7
157
3.a
207
Number of O(n") operations
Wall-clock time (s)
Figure 6: Classic Image problem with L, (top) and L3 (bottom) costs, showing optimality gap vs. number of O(n²) operations (left) and wall-clock time (right).
28

---

Published as a conference paper at ICLR 2025
GRFsmooth
12-
15a
NBOT-mnkeern
Li
NBOTONOS
cost
164
jantu apot (2019)
Srkhhm (p-1..)
Srknom ip-1.55
10-1
MSK
Feydly No.
Optimality (P
10:0
16-
cost
12-
16-
15-9
10a
10
10"
12
10
Number of O(n?) operations
Wall-clock time
Figure 7: CRFSmooth problem with L1 (top) and L3 (bottom) costs, showing optimality gap vs. number (n?)) operations (left) and wall-clock time (right).
GRFmoderate
164
12-4
L
NDSTSminer
NBST-NOG
cost
124
Inburatal et (2019)
164
Srkhom (pm).))
Seknem (p-1.4)
10-1
MSK
10-4
reddy Al.
Optimality
15-9
10-1
10-4
cost
10-4
12
10*
10
10-1
10"
Number of O(n?) operations
Wall-clock time
Figure 8: GRFModerate problem with L1 (top) and L3 (bottom) costs, showing optimality gap vs. number O(n?) operations (left) and wall-clock time (right).
29

---

Published as a conference paper at ICLR 2025
GRFrough
15
123
Li
10-4
NBOTOWOS
cost
juntu (2019)
Srkhhem 1p-1..
124
3rkhom (p-1.55
10-
MSK
Feydly No.
Optimality
12
12-
10-4
cost
10-
16-3
10-4
207
30a
10
10
30
30'
20'
Number of O(n?) operations
Wall-clock time
Figure 9: GRFRough problem with L1 (top) and L3 (bottom) costs, showing optimality gap vs. number (n?)) operations (left) and wall-clock time (right).
LogGRF
10-
L
NDSTSmiter
NDSTPNOG
cost
1
Inburaatte er (2018)
Seknom (p-1.4)
Srkhom (pm).))
NSK
10-4
reddy Nl.
Optimality
10a9
12
12
cost
16-4
12-1
124
10
15a
10
10
10-7
102
Number of O(n²) operations
Wall-clock time
Figure 10: LogGRF problem with L1 (top) and L3 (bottom) costs, showing optimality gap vs. number of O(?²) operations (left) and wall-clock time (right).
30

---

Published as a conference paper at ICLR 2025
LogitGRF
13-
Li
NBOTONOS
cost
164
bunta (2019)
10-
Srkhhm 1p-1..0
Srhnom ip-1.5)
MSK
10-M
Feydly No.
Optimality
10-
10-1
12-4
cost
16-
16-1
10
10
201
10'
Number of O(n?) operations
Wall-clock time
Figure 11: LogitGRF problem with L1 (top) and L3 (bottom) costs, showing optimality gap vs. number of (n²) operations (left) and wall-clock time (right).
Microscopylmages
T-11
16-
L
NDSTSmiier
NDST-PNG
cost
1
Inburaatle er (2018) (201))
Seknom (p-1.40
Srkhom (p-).))
MSK
1214
reddy Al.
Optimality (P
10-
12-1
10
cost
12
12
12a
18*
10
10-7
107
107
107
Number of O(n²) operations
Wall-clock time
Figure 12: MicroscopyIIagge problem with L1 (top) and L3 (bottom) costs, showing optimality gap vs. number of O(n²) operations (left) and wall-clock time (right).
31

---

Published as a conference paper at ICLR 2025
Shapes
12
16a
MBOT-nnkeer
Li
NBOTOWOS
cost
164
juntu aped (2019)
Srkhhm (p-1..0
3rkhom ip-1.55
10-1
MSK
10-
Feytly No.
Optimality (P
10-0
12-
cost
10-
16-3
16-
20
127
20
10
13a1
107
10
Number of O(??) operations
Wall-clock time
Figure 13: Shape problem with L1 (top) and L3 (bottom) costs, showing optimality gap vs. number of (n²) operations (left) and wall-clock time (right).
WhiteNoise
10-2
10-
L
124
NDSTSminer
NDST-NOO
cost
pnburaatien er (2018)
Seknom 1.-1.40
164
Srkhom (p-).))
10-1
MSK
reddy Al.
Optimality
10-9
12
12:4
cost
16-
12-1
107
100
10'
10
10-1
16
Number of O(n?) operations
Wall-clock time
Figure 14: Wh WhiteNoise problem with L, (top) and L3 (bottom) costs, showing optimality gap vs. number of (nn?) operations (left) and wall-clock time (right).
32

---

Published as a conference paper at ICLR 2025
E ON THE SPECTRAL GAP 1 - 22 AND DISCOUNTING
In Section 3, we demonstrated that the convergence of Algorithm 2, as established in Theorem 3.2 (and subsequently in Corollary 3.4), exhibits worst-case dependence of O((1 11/)) on the spectral gap Xg. In this section, we first explain why this dependence may be overly pessimistic. We then illustrate this observation through the experiments presented Fig. 16.
Recall from Lemma A.6 that the convergence of CG for solving p-discounted Newton/Bellman
system du = Suv + dy is in fact O((1 -) In Thm. 3.2 we expressed the overall complexity of Alg. 2 in terms of X2 using the fact that the largest p found by Alg. 2 via annealing satisfies (1 p)
O(1 2²) - in the worst-case; see (19) and the proof of Thm. 3.2. However, 1 p tends to be much better behaved in practice and effectively mitigates this worst-case dependence on 1 A2 We introduce the rationale for this behavior via an example.
Suppose, given n = the stochastic matrix Pre E TR. 20
has the form:
00775
0.25
0
0.25
0.75
0
0
Pre =
Q =
0
0
0.75
0.25
0
0
0.25
0.75
where A R444 is appropriately selected, tiny perturbation matrix with infinitesimal entries, and the Markov process given by Pge can be (approximately) decoupled into two separate processes corresponding to the top-left and bottom-right blocks; the spectral gap of Pre is tiny.
Consider also, for convenience, a stationary distribution r(P) = [0.25 0.25 0.25 0.25] T and
an initial distribution amount of total mass as the stationary distribution r(P). For such a pair, r(P) r would be
[0.5 0 O 0.5] which already assigns each partition the correct
orthogonal to the first eigenvectors V1 [1 1 0 0] and 12 0 1 1 of Q. We can then easily show that the convergence of r to r(P) via repeated application of Pre would be governed by Az)= = 2 rather than (1 A2) O(IIA| Indeed, we can write:
Q vnnj + 12 + 0.5 + 0.5m4u,
which implies (with simple calculations):
[|r(P) -r - Pt [[(((P) - rr p ((+(P) + +II [[(((P) Il) = =0.5'x(rlv(P)) + O(IIA] [|r(P) h(P))l)
The idea here applies more generally to multiple eigenvalues near one, i.e., when the Markov process may be approximately separated into a larger number of groups. Since we call Alg. 2: NewtonSolve near the solution, where V [|r(P) - rll, is small, the gradient (i.e., the RHS of the Newton system) tends to be nearly orthogonal to the subspace corresponding to near-zero eigenvalues of Pre- Consequently, the Newton/Bellman system can be discounted more aggressively than the worst-case suggested by the condition in (19), ultimately yielding (1 - p) < (1 A2)
In Fig. 16, we demonstrate this effect via ex- periments over datasets from the DOTmark benchmark (Schrieber et al., 2017). where dis- play (1 (1-p)) p) 1) (1 x2)-1/2 11/2 and the number of CG iterations until convergence in the same plot. Since P tends to get sparser with increasing generally observe de-
(Tang al., 2024), creasing spectral gap of Pre with higher y (across all problem sets) as communication between par- ticle groups is increasingly bottlenecked. In all cases, (1-p) -1/2 predicts the total number of CG iterations better than (1 (1-32 -1/2 Espe- cially for the Shapes problem, the marginals r, contain lots of near-zeros and often yield non-communicating particle groups. Even though
Figure 15: Sample images from Shapes.
(1 - >2) 1/2 explodes quickly in this case, the method remains robust and (1 1--p) pp- 11/2 continues to be predictive of good performance.
33

---

Published as a conference paper at ICLR 2025
Convergenco Alg and spectral gap (CouchyDersity.. 1,
Convergence of Alg. spectral gap
(CauchyDensttyy. 13
X-
NewtonSolve
Number of iterationss Number its taraioons
nerations
Number of Borations
200
spectral
Convergence of 2 and specital (Clasicmaagss Lj cost)
(Carstdmage, cost)
-G ittrations in NowvonSohre
CS terations NowtonSolve
Numbeer of Berations
1.2
X3
204
12*
15a
12
10a
10a
10a
12a
12'
V
Convergencs 2nd spectall gop (Githough cos))
Convergente
IGAFrough Li cost)
107
CB Iterstions Newtonn
Numb of Itarattons
CG Btrations NewtonSolive
10
100
Cenvergance of and [Shapes, L, cost
Converponce of Aip and spectral g2o
[Shapes E com)
1.7
Derationg Solvon
200
Rerations NowtonSolve
Number ileralionss
10
101
Number
302
167
16'
10
107
22
Figure 16: Over 4 tasks from DOTmark (n 1024) with L2 (left) and L3 (right) costs, we run MDoT--runcatedNeettn with a fixed temperature decay rate of g = 21/16 and compute the spectral gap 1 A2 explicitly via eigenvalue decomposition every-time Alg. 2: NewtonSolve is called (at fixed intervals of on a log-scale). On the LHS of the y-axis, the (1 A2) 12)-11/ 1/2 term seen in Thm. 3.2 and Cor. 3.4 and (1 -1/2 are displayed, where is the discount factor found by Alg. 2. On the RHS of the y-axis the total number of CG iterations taken by NewtonSolve shown. While 11/2 term shown in Sec. 3 may be adversely affected by extremely large values (see
the (1- -1/2 bottom row), the practical rate dictated by 1 - is much better and more accurately describes the convergence of the algorithm. Lines show the median over all 45 problems and shaded areas show 90% confidence intervals. See Appx. E for further discussion.
34


Reviews:
[{"review_id": "1efJxg9mbZ", "reviewer": "Reviewer_9Fpa", "summary": "This paper introduces a truncated Newton method, a second-order algorithm, to solve the Entropic-regularized Optimal Transport problem.", "questions": "Various notations used in the manuscript is not clear or hard to find the definition. For example, what is $\\mu$? Some intuitive explanation of $P_{rc}$ would be nice also. I recommend adding more to the notation section and explain these better.\n\nIn the experiment section, can the authors compare the method with other second-order methods, if any?\n\nAlso, can you give an intuitive explanations of how you do the proofs (main ideas or a high-level analysis like an ODE analysis)? I'm very interested in the results, and maybe in the later iterations of the paper you can incorporate these proof sketches to help the readers also.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The paper has interesting theoretical results, such as the connection to the Bellman equation (though I'm not sure if this is the first work to establish that connection) and the convergence in Chi-square distance of Sinkhorn algorithm (again, I'm not sure if this result is previously established and to the best of my knowledge this is the first time I see this result).\n\nTheorem 3.2 is interesting since it's a log(1/eps) convergence result, though I skeptical of the dependence on other terms could harm the convergence. The authors should discuss this result, as to the best of my knowledge, the only other equivalent result is in this work [1] (UOT can be viewed as a generalization/approximation of OT, I think this work should be cited in discussion to this work as well).\n\nThe paper has strong experiment settings that show superior result on wall-clock time and the color transfer result is very nice. The authors also have a detailed description for an efficient implementation of proposed algorithm.\n\n[1] \"On Unbalanced Optimal Transport: Gradient Methods, Sparsity, and Approximation Error\".\nQuang Minh Nguyen, Hoang Huy Nguyen, Lam Minh Nguyen, Yi Zhou\nJournal of Machine Learning Research, (JMLR), 2023", "weaknesses": "The complexity in Theorem 3.2 depends on the spectral gap of P_{rc}, but there is no control over this spectral gap. I suspect that this could lead to very poor convergence should the spectral gap is very bad. I recommend having some analysis on this term to make the results stronger. Additionally, it is possible for this term to be $n$ dependent too, so the $n$ dependence of the complexity might more than just $O(n^2)$. If that this is the case, the method is a tradeoff between $n$ and $\\epsilon$.\n\nThe paper is a bit hard to read and it is hard to follow the intuition of the paper. The proof should have been explained better, such as the chi-square convergence. Furthermore, the notations are not clearly explained or presented (see below).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qCcEBuINHg", "reviewer": "Reviewer_A2Da", "summary": "The paper introduces a truncated newton method for solving entropic optimal transport problem. The algorithm is GPU-parallel and is scalable for large OT problem.", "questions": "1. Could the author make the theorems self-contained? For example, what is $\\lambda_2$ in Corollary 3.4?\n2. What does helper Sinkhorn mean in Line 196?\n3. The definition for LSE is not clear in Line 066-067.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The proposed algorithm has good practical performance under GPU settings. \n2. The author has done a good job in discussing the related work.", "weaknesses": "1. Please clarify the definitions of $\\epsilon$ and $\\epsilon^\\prime$. Current presentations are confusing.\n2. The proposed algorithm lacks a theoretical convergence analysis and a comprehensive assessment of its computational complexity. While popular OT solvers like Sinkhorn and APDAGD offer promising guarantees on total computational complexity, the theoretical complexity of this algorithm remains unclear to me.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "S7ZptoWmOP", "reviewer": "Reviewer_4av9", "summary": "The paper presents a new Newton-based optimization method for entropic regularised optimal transport. The authors present a new three-stage method that combines truncated Newton method, conjugate gradient and mirror descent. Finally, the authors present experiments with fast convergence of the proposed method.", "questions": "see weaknesses", "limitations": "", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "In the proposed framework, the authors propose a lot of engineering improvements to make the method more stable and effiicient from the practical point of view.", "weaknesses": "The Weaknesses and Questions are together.  \n\n1. The section 2.1. is concerning and confusing.  First of all, lines 70-92 are almost exactly copied from  Kemertas et al. (2023), which is concerning. Secondly, the problem 2 is wrongly defined. It is not convex. It should be $+H(P)$. Thirdly, The dual problem is also unusual for me. Normally, the dual problem has a form of LogSumExp. In the paper, it is the sum of exponents. Please, provide how do you get formulas 3 and 4 in Appendix. \n\n2. The major results from my perspective are more focused on practice and involve a lot of engineering improvements. Hence, the major proof of success should be experiments. However, from my perspective, the experiments are not good enough for a practical paper. a) First of all, there is only one figure with 2 examples, which is not enough. b) Secondly, in Fig.2, the proper loss should be a duality gap, not a primal function value. As the primal function value excludes constraints' violation and can be misleading. Also, the figure with feasibility may also help in Appendix, as well as for the dual loss. c) The presented results in Fig.2 are presented without statistical measures. It means that it was only one run, which can be quite stochastic, especially with respect to wall-clock time. d) All methods are starting from different points, which is very concerning. e) The results are only presented with respect to wall-clock time. I would recommend adding another measure that is not based on code performance, for example number of calculated gradients or something like that. Otherwise, MDOT code can be optimized for performance and compete with other methods that are not optimized. f) Finally, I didn’t find any attached code to verify the performance.\n\n3. Do I understand correctly, that at the beginning, in practice, there is a part of pure adaptive mirror descent to find a near solution as a starting point, and then the algorithm performs Algorithm 4?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)