A Memory-Efficient Hierarchical Algorithm for Large-scale Optimal Transport Problems
Authors: Wenzhou Xia
Ya-Nan Zhu
Jingwei Liang
Xiaoqun Zhang
Venue: ICLR 2026 Poster (2026)
Keywords: optimal transport
linear programming
multiscale framework
first-order methods
Primary area: optimization

Abstract: We propose HALO, a memory-efficient hierarchical algorithm for solving large-scale optimal transport (OT) problems with squared Euclidean cost, particularly effective in moderate-dimensional settings. 
The core of \ours lies in combining a hierarchical representation of the OT problem with parallel-friendly linear programming solvers, within which an active pruning technique is integrated to further reduce memory usage and computational cost.
Theoretically, we establish a scale-independent iteration-complexity upper bound for the refinement phase, which is consistent with our numerical observations. 
Numerically, experiments on the image dataset \dataset and the 3D point cloud dataset \datasetnongrid demonstrate that \ours effectively alleviates the memory and scalability bottlenecks of existing solvers.
Our method demonstrates significant advantages compared to state-of-the-art baselines: for images with $n=1024^2$ pixels, it achieves an $8.9\times$ speedup and $70.5$% reduction in memory usage under comparable accuracy; for 3D point clouds at scale $n=2^{18}$, it achieves a $1.84\times$ speedup and an $83.2$% reduction in memory usage with $24.9$% lower transport cost.

Published as a conference paper at ICLR 2026
A MEMORY-EFFICIENT HIERARCHICAL ALGORITHM FOR LARGE-SCALE OPTIMAL TRANSPORT PROBLEMS
Wenzhou Xia', Ya-Nan Zhu Jingwei Liang 3.
1..t, Xiaoqun Zhang!.3.4;
ABSTRACT
We propose HALO, a memory-efficient hierarchical algorithm for solving large- scale optimal transport (OT) problems with squared Euclidean cost, particularly effective in ing a hierarchical representation of the OT problem with parallel-friendly linear programming solvers, within which active pruning technique is integrated further reduce memory usage and computational cost. Theoretically, we establish a scale-independent
settings. The core of HALO lies in combin-
complexxty upper bound for the refinement phase. which is consistent with our numerical observations. Numerically, experiments on the image dataset DOTmark and the 3D point cloud dataset ModelNet10 demon- strate that HALO effectively alleviates the memory and scalability bottlenecks of existing solvers. Our method demonstrates significant advantages compared state-of-the-an baselines: for images with n 1024 pixels, it achieves an 8.9x speedup and 70.5% reduction in memory usage under comparable accuracy; for 3D point clouds at scale 2'8, it achieves 1.84x speedup and an 83.2% reduction in memory usage with 24.9% lower transport cost.
INTRODUCTION
The Wasserstein distance, defined via the Kantorovich formulation of Optimal Transport (OT) prob-
lems, a powerful metric for measuring the similarity between two probability distributions.
has been widely adopted in various fields such as generative modeling Arjovsky et al. (2017); Ko-
milov et al. (2024): Tong et al. (2023); Hui et al. (2025), color transfer Solomon et al. (2015); Pitié &
Kokaram 2007), texture synthesis and mixing Dominitz & Tannenbaum (2009) Rabin et al. [201]]
registration and deformation Haker et al. (2004); Rehman et al. (2007). image restoration He et al.
(02T), , domain adaptation Montesuma & Mboula (2021): He al. (2024) and transportation-bant bbasdd morphology metrics Basu et al. (2014). However, modern large-scale applications push OT solvers regimes where classical formulations become computationally prohibitive and memory intensive.
Given probability measures on domains S and D with ground cost e : S x D + R, the Kan- torovich problem defined by
inf ( c(s, d) ds(ssd),
=Ell(u,)) LexD
(1)
where [[p, v/) denotes the set of couplings on S with marginals μ and on and ID, respectively.
In the discrete case, the supports are S = (s Sm } and ID = [di...., dx with corresponding probability vectors u T and e R² defined by u, p(s,) and Uy v(d,). A coupling represented by a nonnegative matrix X mx whose row and column sums equal u and v. The ground cost encoded in the cost matrix C = (c(s;, dz)) E Rm=x Consequently. this setup leads to the following discrete OT problem in standard Linear Programming (LP) form:
min e x;, S.I. = 220
(2)
School of Mathematical Sciences, Shanghai Jiao Tong University, Shanghai, China; School Math- ematics, Harbin Institute of Technology, Harbin, China:" Institute Natural Sciences, Shanghai Jiao Tong Institute of Artificial Intelligence, Shanghai, China. Emails: xw2. hhh@s jtu jtu.jeuuu. stu.edu.cc. en, yananzhu6@gma com, jinwwei. 1 i ang@s jtu. edu. en, xqzhangesjtu. edu. en. Corresponding authors.
University, Shanghai, China:

---

Published as a conference paper at ICLR 2026
1 8 1m In 1
where x =ee(X), vec(C)) q = - ( and A
E IR (
This LP problem has mn variables and m + constraints. Considering OT problems between two images of size n = X r, the variable count in X scales as n², reaching approximately 10 when 1024. Such an enormous scale results in prohibitive computational cost and memory consumption. Moreover, traditional solvers such as the network simplex (aboo & Tarjan (1991) and interior-point Pele & Werman (2009) methods further fail to exploit the modern GPU architec- tures, leading to growing gap between classical OT algorithms and the scalability demanded in contemporary applications.
Related work Existing algorithms for large-scale OT fall into three broad categories:
The first category consists of approximation methods, which trade accuracy for efficiency. Represen-
tative examples include entropy-regularized algorithms such as Sinkhorn Cuturi (2013); 2 D
sky et al. () Schmitzerr (2019); [ CCen et al. (2022), low-rank algorithms Scetbon et al. (2021);
Halmos al. (2025) as well as approaches based on approximate metrics such as the approximated earth mover's distance Shirdhonkar & Jacobs (2008) and the sliced Wasserstein distance Kolouri
et al. (2019); Nadjahi et 1202 These methods achieve good scalability but at the cost of re-
duced accuracy.
The second category comprises LP-based algorithms, which aim to accelerate OT by improving LP solvers tailored to this structured problem. Representative methods include randomized block coor-
dinate descent methods Xie al. (2024). semi-smooth Newton-typo approaches Li al. (2020). and
first-order algorithms such as Douglas Rachfos-Rachford splitting Mai et al. (202T), dual extrapolation Jam-
bulapati al. (2019),
methods Zhang et al. (2022) Chen al. (004)).
Halpern iteration Zhang et al. (2025), and primal-dual hybrid gradient (PDHG) Esser et al.
2010
Chambolle & Pock (2011). Among them, factorization -free first-order LP solvers, such as PDHG
and its variants Applegate et al. (2021); Du[L et al. (2023 2025), are particularly attractive for large- scale OT, since their computations are dominated by matrix-vector multiplications and are therefore
well suited for GPU parallelism. Nevertheless, despite these advances, LP-based solvers continue to suffer from severe memory bottlenecks in uatra-large-scale settings Lu & Yang (2024).
The third category consists of multiscale algorithms that exploit the pyramid structure of Mérigot (201T); ; Gerber Maggioni (2017): Leclaire & Rabin (2019), with representative methods such as ShortCut [Schmitzer] (2016) and the multiscale semi-smooth Newton (MSSN) method Liu cet al. (2022). Although these methods effectively reduce computational cost in practice, the lack of iteration-compleiit bounds leaves their worst-case efficiency unresolved, and the dependency on CPU-based LP solvers limits their scalability on modern parallel hardware.
Overall, developing OT solvers that simultaneously meet the demands of memory efficiency, large- scale parallel computation, and high accuracy remains an open challenge.
Main contributions We propose HALO, a memory-efficient and parallel-friendly Hierarchical Algorithm for Large-scale Optimal Transport problems with squared Euclidean cost, which is par- ticularly effective in moderate-dimensional settings. The key idea is to exploit the hierarchical struc- ture of OT by solving a hierarchy of OT problems across resolutions: coarse-level solutions warm- start finer levels, and each level alternates between updating active support and solving a restricted LP the active support with GPU UPU-friendly solvers.
The significance of HALO lies in the following aspects:
Memory efficiency: HALO requires only ((n) memory, matching (to the best of our knowledge) the lowest space complexity among existing GPU-based OT solvers.
Parallelism: HALO uses LP solvers based on factorization- first-order methods, whose domi- nant operations are matrix-vector multiplications, enabling efficient parallel implementation.
Theory: Under certain assumptions, HALO enjoys a -ndeeinndeeeddent iteration-complextty bound per level, which also validated empirically.
Empirical performance: On DOTmark at n = 10242, HALO achieves an 8.9x speedup and 70.5% GPU-memory reduction; on ModelNet10 with n = 218 it achieves a 1.84x speedup, an 83.2% 2% GPU-memory reduction, and 24.9% lower transport cost.
2

---

Published as a conference paper at ICLR 2026
The rest of the paper is organized follows: In Section introduce the pro- posed algorithm, where Section 2.3 provides the ccale-independentt lterationccompeext bound. In Section we present the experimental results, and Section concludes the paper.
Organization
Notation S, D denote the source and target spaces in OT problem. In discrete setting, given {(ste dz.) 1 < & < } C S DD As denotes the sub-matrix of the constraint matrix A restricted the columns indexed by N, and ay as the corresponding sub-vector the vector R Let Topx(C) denote the operator that selects K pairs from the set C C $ x ID, where these pairs have the top associated values.
2 HALO: A HIERARCHICAL ALGORITHM FOR LARGE-SCALE OT
2.1 OUTLINE THE PROPOSED FRAMEWORK
We first outline the overall workflow of HALO. The key reason behind the efficiency of HALO lies in two aspects: (i) an inherent multiscale structure that enables a coarse-to-fine hierarchical represen- tation, and (ii) sparsity of the transport plan which allows to incorporate active-support detection to the solver. These two aspects jointly contribute to the significant reduction of the computational complexity and memory demand.
Scale
solveRestricted
PDHG Salver
updateActive
prolongation
0(1)
Scale L-1
updote
N
solve
updote/salvee holle
solve
coarse-to Porse-to-fine fo-fine prolongation
X O(1)
Scale
N
update
solve
/anive/
salve
Output
Figure 1: Architecture of HALO.
Figure Sesureeprovee provides a high-level overview of our HALO method:
Hierarchical Structure (Outer Loop over Levels): We construct an (L + 1)-level hierarchy of OT problems, ranging from the coarsest L-th level to the finest 0-th level. For each level & E 10.,, L). let (u(") v(() c(e) denote the marginals and costs, and (f((), g(()) denote the dual potentials. Based on the I.P formulation (2), the equivalent min-max formulation of the OT problem at the E-th level reads:
min ER meee g( (O"R"" ( eRa
max
(1) p(0) + (0 (
(0)) + (c(9,x10) (0) (f(9) g (1))
(3)
where ( is the OT constraint matrix level f, and my, ny denote the marginal scales. We leverage this hierarchy via a coarse-to-fine strategy: For f = L, solve problem (3) and denote the
solution as (x (z(), f(t) .g" (T));; for E (0.... L 1}. the solution (z!++1), f f((+)), (((F)))g(((())) ((++))) from the coarser (E + 1)-th level is prolongated to initialize (x()), {(), (f)) for the finer L-th level: see Section 2.2 for details.
3

---

Published as a conference paper at ICLR 2026
Scale
Sodlc
Figure 2: Left & Middle: in updateActive green region indicates possible unshielded d, and arrows denote pairs (s, d) added by the shielding-based component.
Hierarchy construction of images and point clouds. Right: Shielding
Refinement Iteration (Inner Loop at Fixed Level): At each specific level E fo.... 1}, we refine the initialized solution to solve the OT problem (3). To fully exploit the sparsity of the transport plan, we alternate between updating the active support and solving the corresponding LP problem using GPU-based solvers (instantiated by PDHG-based algorithm Lu et al. (2025)): see Section 2.3 and Algorithm for details.
Summarizing the above discussions leads to the proposed HALO method.
Algorithm HALO: Hierarchical Algorithm for Large-scale Optimal Transport
1: Input: marginal distributions u and ground cost number of levels L.
2: Output: optimal coupling (")
3: Build hierarchical OT problems {uu(), (((), Ho
4: Solve OT with (u(2), v(L), c(l1) on the coarsest level L to obtain ( ( g(4).
5: for 8= T. , do
6: Initialize (x) ((), (()) by prolongating (I++1) f(++1) g (++1)).
7: repeat
8:
1). Updating active support, detailed in Algorithm 3
N(0) + sopdateActive(tt ( ((),g(0 (0) N()
9:
10:
2). Solving restricted OT, detailed in Section 2.3
11:
(x(*), ((t), (g) + so vollveeestrii retsstriiteiii. (2)
(0) 10, u00
v (00 )())
12: until (0) meets the termination criteria.
13: end for
14: return
For the rest of this section, we discuss the details of the key steps in Algorithm including the construction of the hierarchy, the active-support updating strategy and its convergence property.
2.2 HIERARCHY: COARSE-TO-FINE MULTISCALE STRUCTURE
A high-quality initialization greatly accelerates the inner loop in Algorithm | To this end, we exploit the geometry OT and compute such an initialization via a coarse-to-fint multiscale scheme.
Consider discrete OT problem with supports $ {*1,-.. " ID (di.....dd, and marginals μ = Σn u,02, 2, // = Σ v. In our hierarchy, the supports at the é-th level are formed by representative points of corresponding neighbor groups from the (e - 1)-th level. Accordingly, the marginals (e) and (f) are obtained by aggregating the masses of these constituent groups, and the cost ( is defined by the distance between representative points.
As illustrated in Figure 2 in the image setting where $ and D correspond to pixels on regular grid, the hierarchy is built by recursively merging x pixels, using the barycenter of each block its representative. In non-grid settings (e.g.. point clouds). we instead emplov spatial partitioning structures, such as 2d -r-tr or kd-trees Finkel Bentley [974]; Bentley (1975), to construct the hierarchy, detailed in Appendix C.I
4

---

Published as a conference paper at ICLR 2026
2.3 SPARSITY: SOLVE RESTRICTED OT ON ACTIVE SUPPORT
Classical LP theory [Bertsimas & Tsitsiklis (1997) states that the optimal solution x" E R of nonzero entries, which sparse when m, n are large. If the support supp(æ") were known a priori, we could simply solve the problem restricted to this set. Therefore, letting active support N x D be an estimation of supp(x"), ) we define the restricted OT problem on N as follows:
problem 2 has most
Definition (Restricted OT on active support N)
min CN Asaz = 9. EN 0 Eye 0,
(4)
where Ay is the submatrix of A restricted to columns indexed by N, CN is the corresponding sub- vector of c, vy and Fya are the corresponding subvectors on N and NO.
Remark Since the OT problems at each level share the same formulation and differ only in scale, omil the level superscript (f) and simply use the form (2) for illustration.
the active support N exactly corresponds to supp(æ"), solving the restricted problem provides the global optimizer of the original OT problem. However, finding supp(æ") ) is challenging as solving the problem. To address this difficulty, we consider Algorithm 2 a refinement procedure that alternates between updating the active support based on the current solution and solving the restricted OT on the newly updated active support.
Algorithm A refinement at each level
INPUT: initial coupling x! and dual potentials (f") go);; u, and e. OUTPUT: global optimizer
+ 0, No + 0
repeat
N++1 +- updateActive(e"" if g", f 1 colveRestrictediaa f*, N++1; u, & ++
until x* meets the termination criteria
return N**
The efficiency Algorithm hinges on the design updateActive. previous work, Short [Schmitzer (2016) constructs active support by augmenting the support of the previous coupling with shielding-based component. This component relies on the shielding con- dition, defined as follows:
Definition 2 (Shielding condition) Given S, Di and cost e S ID + R. Let ES, ED. We say (s', d') shields s from d if
c(s, d) c(s', d') > c(8, + c(s',d)
To obtain provable (1)) per-level iteration bound and improve robustness in practice, we intro- duce two modifications: (i) expanding the active support to include the entire previous ones, rather than only the previous coupling's support, and (ii) performing dual-violation correction, adding pairs suggested by significant violations in the dual potentials. The complete update procedure is summarized in Algorithm.
Although the first improvement makes the estimated active support larger, provides a scale- independent iteration-ommplexity bound, detailed in Theorem [ ensuring that the growth of the active support does not affect its sparsity. Additionally, factorization-free and farenlel-frienll LP solvers, which is less sensitive the problem scale com- pared to traditional LP solvers. As a result, our algorithm achieves excellent practical performance while benefiting from theoretical complitn guarantee.
the implementation, we use the
The second improvement involves a dual-violation correction mechanism that accelerates conver- gence by adding pairs suggested by significant violations in the dual potentials. This approach

---

Published as a conference paper at ICLR 2026
ensures that the active support includes supp(æ") much as possible, while maintaining sparsity. As a result, the algorithm refines the solution more efficiently, particularly in some challenging in- stances. Unlike MSSN Liu al. (2022) which relies on threshold parameter, employ a Topx operator to select pairs with the largest dual violations, ensuring the active support remains sparse throughout.
Algorithm 3 updateActive
INPUT: current coupling æ vee (X), dual potentials (f,g), cost E Rum with e: S D - R, active support N, dual-violation hyperparameter /3 OUTPUT: updated active support N'
1). start from current active support
N' + N
2). shielding-basee based
construct map: S + IB, map(s:) = arg maxd,er Xij for all 8; ES
Define R(s) C as the local neighborhood of (8-neighbors for images and KNN for point clouds, Appendix C.2).
for s E S do
(ss) + map(s') E R{{))
D(s) + (dell D: is not shielded from & by (8', map(s')), Aal E R(s))
if D(s) = ø then
D(s) arg min DEDDDED c(map(s),d)
end if
(s)) + (ss) U (ss) ((s)
N' - N' U ((s,d) de D(s)
end for
3). dual-violation correction Pricing problem, C +- ((8,, dz) EEX 1: h. > ag) K + [S) CK +- (C)
N' + CKCK
return N'
Theorem (Scale-independent iteration-complexity bound) Let (S, (S, p.v) be discrete OT in- stance with squared Euclidean ground cost. Consider Algorithm with the update rule from Algo- rithm 3 which produces solutions {(x* f* ))20 and active supports (N* Jk20 via the mono- tone update
N++1 < updat
Assume there exist constants q (0, oo, p< oc, L < ac, Ro < such that:
1. (Directional coverage) For vs E S and Vo E R 3s E S(s) C S with (u, s' - -) 2
||||| [|s' >|1
2. (Bounded radius of S(s)) ||8/ - ا|l < for all s and all s' E S(s).
3. (Bounded density of D) DD n Bz S (Bn) for all z E R. and R > 0.
4. (Uniform Lipschitz regularity) ||map, - mapx (l'll S L|[s Il, Vs, 8' E ES 0.
5. (Coupling stability) (s) map, [s) < Ro for all k and all E S.
Then there exists a constant > such that ay is a global optimizer for all k 2 C, and N S C|S for all k 0.
hold in our setting; by contrast, Assumptions 45 are stronger yet natural under the multiscale warm-start. Since each level refines a warm start prolongation from the previous coarser level, the per-level updates mild, making uniform Lipschitz condition and a bounded drift stability assumption intuitively justified. The proof of Theorem is deferred to Appendix B.l
Remark 2 Assumptions

---

Published as a conference paper at ICLR 2026
At the end of this subsection, we note that in Algorithm 2 the restricted OT the active support is solved by solveRestricter
which dominates computational cost and memory usage, making the LP solver choice critical. To achieve both memory efficiency and parallelism, adopt a factorization-free first-order method, and we highlight the following property of the restricted OT under Pock Chambolle rescaling Pock & Chambolle (201T)
Proposition 1 Let B be the constraint matrix of restricted OT problem, and define the rescaled matrix B by the Sack-Chambolle rescaling:
D, = diag( (vri..... v/ Voi, Von,, D B= BDD,'DD,1
with Tiy Cj the row/column degrees of B. Then B 12
This result justifies the of a constant stepsize in PDHG-based algorithm, eliminating the need for norm estimation and reducing overhead in HALO. The proof is provided in Appendix B.2
3 EXPERIMENTS
To evaluate the scalability and efficiency of HALO, we conduct experiments on two representative datasets: the 2D image dataset DOTmark (Schiieerr et al. (2016) and the 3D point cloud dataset ModelNet10 Wu et al. (2155). For image experiments, we compare HALO against three state-of-the art solvers: HCT Zhang et al. (2025). ShortCut Schmitzer 2016), and M3S Chen et al. (2022). For non-grid experiments, we compare against the state-of-the-art solver HiRef Halmos al. (2025)
and the standard Sinkhorn Cuturi ) Flamary et al. (202I)) Detailed dataset settings are provided in the Appendix E
All experiments run on dual Intel Xcon Gold 6330 CPUs (2.0GHz), 503GB RAM, and an NVIDIA RTX 4090D (24GB). Host memory is limited to 100GB and GPU memory to 24GB (violations reported as OOM). Each instance has a 3600s wall-clock limit (violations reported as TO).
For HALO, the flexibility of the LP solver choices are provided in Appendix H For the other baselines, we use their open-source implementations with default parameter settings, with more details provided in Appendix D
cuPDLPx Lu al. (2055) by default, while additional results demonstrating
Since different methods may use different internal stopping rules, we report solution quality using unified metrics:
gap |/c, - (c, a.5)
and feas max |mii 0)ll |Ax al
Re. an)) +1
+ [r
|1 }
where 26 is a high-accuracy reference obtained by solving the reduced OT model with Gurobi (Barrier with crossover) at a tolerance 10 following Zhang et al. (2025).
3.1 THE IMAGE DATASET DOTMARK
For image dataset DOTmark, we use denote the image resolution; the problem size thus scales as n m = 2 Table summarizes four metrics: runtime (s), peak memory (GB), relative objec- tive gap, and feasibility error. While providing comparable solutions in terms of gap and feas, HALO is memory-efficient and delivers short wall-clock time, clearly outperforming strong baselines, es- pecially at large scales. Against HOT, HALO leads on runtime, memory, and gap. Specifically, at = 512 it is 7.02 x faster with 89.2% less memory, achieving 2-3 orders of magnitude tighter gaps across scales. Despite ShortCut's is CPU memory efficiency, HALO is much faster at a comparable gap. yielding 37. x speedup at 512. Finally, while the entropy-regularized M3S also solves the r 1024 instances, HALO shows simultancous advantages runtime, memory, and accuracy. At 1024 is 8.92x faster, uses 70.5% less memory, and attains 1-2 orders tighter gaps.
Figure reports runtime and memory scalability on log-log axes. HALO shows an approximately straight runtime curve with slope 1, indicating near-linear growth and outperforming other base- lines. Its memory curve has a terminal slope x 2, matching M35, yet HALO maintains the lowest memory usage among GPU-based methods across resolutions.
7

---

Published as a conference paper at ICLR 2026
Table 1: The numerical results on mOTmark. "GPU/CPU memory" denotes GPU VRAM for GPU- based methods and CPU RAM for CPU methods (shown in gray). Time is reported in seconds (s) and memory in gigabytes (GB). gap at 512, 1024 are unavailable because solving the reduced model with Gurobi runs out of memory.
Metric
Resolution
64
128
256
512
1024
HATO
1.50
2.20
4.31
11.17
27.73
HOT
1.56
2.12
14.32
78.43
OOM
time
ShortCut
0.25
2,41
25.74
438.14
TO
M3S
1.95
3.44
8.51
39.32
247.22
HALO
0.38
0.48
0.76
2.07
6,25
GPU/CPU
HOT
0.84
1.09
3.10
19.25
OOM
memory
ShortCut
0.01
0.04
0.14
0.56
TO
M35
0.71
0,95
1.92
5,78
21.21
HALO
1,23E 6
1.51E-5
1.41E-5
HOT
6.77R-4
6.03E-3
-
gap
ShortCut
1.56E-6
4.04E-6
2.27B-5
M3S
1.89B-4
1.90E-4
3.31E-4
HALO
3.66E-8
2.45B-7
1.28E-7
1.06E-7
6.98E-8
HOT
5.51E
7.91E-7
7.61E
3.12E
OOM
feas
ShortCut
1.65E
1.65E-18
18
9.55E 19
3.67E 19
TO
M3S
2.40E-6
1.20E-6
6.13E-7
3.03E-7
1AME-7
Time Resolution
Memory Resolution
10a
(a=0.95)
104
HOT
HOT
M35
M35
ShoriCut
Memory (GB)
ShoriCut
10
109
10-1
10
10-2
24
25
2 Resolution
2
29
210
2
2
2
2 Resolution
2
2
210
(a) Runtime scaling.
(b) Memory scaling.
Figure 3: Scalability DOTmark across all evaluated methods. Left: Runtime scaling; Right: Memory scaling. In the left legend, a is the fitted slope for HALO. ShortCut uses CPU memory. shown in gray in the right panel.
To validate Theorem 1 Table lists the average number of inner iterations per scale Algorithm 2 across resolutions. The averages never exceed 2 and even tend to decrease at finer scales, corrobo- rating the (11) bound and illustrating the practical efficiency of updateActives
Table 2: Number of inner iterations per scale.
Scale
16
32
128
256
512
1024
scale0
1.00
1.91
1.82
1.22
1.12
1.05
1.00
scalel
1.00
1.80
1.65
1.15
1.09
1.05
scale2
1.00
1.75
1.63
1.17
1.09
scale3
1.00
1.78
1.66
1.17
scale4
1.00
1.77
1.66
scale5
1.00
1.78
scale6
1.00
8

---

Published as a conference paper at ICLR 2026
To understand the impact of data geometry solver efficiency, analyzed the runtime of HALO across different image classes in DOTmark (see Figure ! for examples). As shown in Table S the runtime is highly correlated with pixel intensity sparsiry, defined as the percentage of pixels with strictly zero mass. Classes with low pixel intensity sparsity, such as Classiclmages, converge fast, while those with high pixel intensity sparsity, such as Shapes and Microscopy, require significantly more time. This phenomenon can be theoretically explained by Theorem [ high pixel intensity sparsity often corresponds non-convex supports with singularities Luto al. (2022), which implies relatively larger Lipschitz constant L in Assumption 4 resulting in a larger constant in the iteration bound and thus longer runtime.
Table 3: Performance breakdown by image class in DOTmark at resolution 1024 x 1024. Metric sparsity denotes the pixel intensity sparsity, defined as the percentage of pixels with strictly zero mass. Time is reported in seconds (s)
Metric
WhiteNoise
GRF / Log (Avg)
Classiclmages
Cauchy
Shapes
Micro.
sparsity 0.00%
0.00%
0.01%
0.00%
45.3%
42.0%
time
18.29
21.76
23.78
25.82
44.39
56.22
Table presents an ablation that isolates the effects the multiscale framework and the GPU- based LP solver cuPDLPx. When cuPDLPx is disabled, use Gurobi 's barrier with crossover, as uwdateeeeee. teAt i ve relies on the sparsity of solutions. Disahling CuPDLPx in HALO results in 36.9x increase in runtime at
256. Removing the multiscale framework from HALO also causes an 85.6x slowdown at r = 64 and leads to OÓM at higher resolutions, Taken together, the multiscale framework and cuPDLPx are both indispensable, yielding short wall-clock time and low memory across all tested resolutions.
Table Ablation on HALO. An 'X the PDHG-based column indicates that cuPDLPx replaced by Gurobi's 's barrier method with crossover.
Multiscale
PDHG-based
Resolution
32
64
128
256
time
0.88
1.50
2.19
4.31
GPU memory
0.38
0.39
0.48
0.76
time
0.91
2.56
27.68
159.06
CPU memory
0.07
0.07
0.29
1.18
time
2.87
128.4
OOM
OOM
X
GPU memory
0.54
3,00
OOM
OOM
time
6.88
126.94
OOM
OOM
x
CPU memory
0.78
12.69
OOM
OOM
Table further further presents an ablation of the dual-violation correction in updat ve. We report the maximum and mean runtime over all DOTmark instances. The augmentation markedly improves robustness on difficult cases; at r = 1024, for instance, the maximum runtime falls to 24.3% of that without this component.
Table 6 reports the speedup of HALO with constant stepsize over the power-iteration choice in CuPDLPx. We show results for resolutions from 256 1024; at r 1024, the constant step- size yields 1.97 speedup, indicating that eliminating per-iteration norm estimation substantially reduces computational cost.
3.2 THE POINT CLOUDS DATASET MODELNETIO
To validate generalization, we evaluate HALO on the 3D point cloud dataset ModelNet10 (see App. C.3 for 2D results). As shown Table 7 HALO demonstrates superior scalability and effi- ciency non-grid domains. While Sinkhorn runs out of memory
2 and HiRef fails at = 219, HALO successfully scales to n = 29 consuming only 2.99 GB memory. Specifically at 218 HALO achieves 1.84x speedup and 83.2% reduction in memory usage compared to HiRef. Crucially, HALO attains significantly lower transport cost, improving upon HiRe f by approximately 24.9%, demonstrating the superior precision of HALOO

---

Published as a conference paper at ICLR 2026
Table 5: Ablation on updat erct ive: dual-violation augmentation (V) improves robustness
Table Speedup of constant step-size over power iteration.
dual-violation
Resolution
256
512
1024
Resolution
256
512
1024
Max time
11.99
54.36
154.20
Speedup 1.80x
1.65x
1.97x
Avg time
4.31
11.17
27.73
Max time
29.07
151.98
633.34
Avg time
5.47
15.23
52.85
Table 7: Performance Non-Grid 3D Data (ModelNet1O). gap denotes the relative objective differ- ence: for n = 214, the reference is the exact solution computed by the standard EMD solver Flamary al. 2021 ; for 2 where EMD solver is intractable, the reference is the solution of HALO.
Metric
Method
2
2
2
2
2
2
HALO
15.54
26.25
47.42
88.51
229.7
444.3
time
Hikef
23.60
46.40
94.40
189.9
422.7
OOM
Sinkhorn
29.50
OOM
OOM
OOM
OOM
OOM
HALO
0.54
0.64
0.77
1.12
1.83
2.99
memory
EiRef
0.91
1.03
1.67
3.60
10.92
OOM
Sinkhorn
10.60
OOM
OOM
OOM
OOM
OOM
EALO
+5.33E
gap
EiRef
+4.19E 1 +3.51E-1
+3.13E-1 OOM
+2.77E-1 OOM
+2.49E-1 OOM
OOM
Sinkhorn
+3.20E-2
OOM
OOM
DISCUSSION
In this paper. we presented HALO, scalable and memory-effiiint solver for large- optimal transport problems with squared Euclidean cost. Our work offers key insight into large-scale OT: by synergizing a hierarchy framework with GPU-based LP solvers, is feasible simultane- ously achieve memory efficiency, computational speed, and high precision, without resorting to regularization or approximation methods. Specifically, HALO integrates a hierarchical frame- work with rigorous active-support update rule and a parallel-friendly LP solver. Furthermore, the duol-violtiion correction enhances robustness, a proven sale-independent iteration bound explains the fast convergence, and the Pock-Chambolle rescaling justifies the rationality of constant stepsize.
Regarding the extension to high dimensions, although HALO faces challenges scaling to high dimen- sions due the growing shielding based component, feasible direction is prioritize and further optimize the dimension-independden dual-violation correction. Our preliminary results suggest that the algorithm holds great potential for scaling high dimensions in the order of thousands.
Finally, for general transport costs (e.g., Euclidean cost), a current limitation is the reduced sparsity of solutions returned by first-order solvers compared the squared Euclidean case, which dimin- ishes the efficiency of the refinement phase. To address this, potential solutions include integrating future GPU-based crossover algorithms to recover sparsity designing more flexible active-support detection strategies.
10

---

Published as a conference paper at ICLR 2026
ETHICS STATEMENT
We confirm adherence to the ICLR Code of Ethics. This paper develops a hierarchical algorithm for large-scale optimal transport and does not involve human subjects, clinical data, personally identifiable information. All datasets used are publicly available under their respective licenses; we applied only standard preprocessing and did not attempt re-identification or attribute inference. The method is general-purposs. so any faimesss concerns stem from downstream data and deployment contexts; practitioners should audit subgroup performance their applications. To reduce misuse risks, we restrict experiments to benign public datasets and will document intended use and limita- tions. Experiments were run on local institutional hardware without external data services; compute configuration and runtimes are reported in the appendix to encourage reuse and minimize redundant computation. To the best of our knowledge, this work complies with applicable laws and institu- tional policies; IRB approval was required. The authors declare no conflicts of interest and no third-party sponsorship that could unduly influence the results.
REPRODUCIBILITY STATEMENT
Complete proofs are provided in the appendix: Theorem i proved in Appendix B. and Proposi- ii proved in Appendix B.2] Implementation details for all baselines are given in Appendix Dataset-related details are provided in Appendix
ACKNOWLEDGEMENTS
This work supported by the NSFC (grant no. 12090024), the Natural Science Foundation of Chongqing, China (CSTBR023NSCQQL1XX0055
CODE AVAILABILITY
Our implementation of HALO is available at
REFERENCES
David Applegate, Mateo Díaz, Oliver Hinder, Haihao Lu, Miles Lubin, Brendan O' O'Donoghue, and Warren Schudy. Practical large-scale linear programming using primal-dual hybrid gradient. Advances in Neural Information Processing Systems, 34:20243-20257, 2021.
Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In International conference on machine learning, pp. 214 214-223. PMLR, 2017.
Saurav Basu, Soheil Kolouri, and Gustavo K Rohde. Detecting and visualizing cell phenotype differences from microscopy images using transport-based morphometry. Proceedings of the National Academy of Sciences, (11(9)344883455 2014.
Jon Louis Bentley. Multidimensional binary search trees used for associative searching.
Communications of the ACM, 18(9)500-5517, 1975.
Dimitris Bertsimas and John N Tsitsiklis. Introduction linear optimization, volume 6. Athena scientific Belmont, MA, 1997.
Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications imaging. Journal of mathematical imaging and vision, 40(1):120-145, 2011.
Kaihuang Chen, Defeng Sun, Yancheng Yuan, Guojun Zhang, and Xinyuan Zhao. Hpr-lp: An im- plementation of an hpr method for solving linear programming. arXiv preprint arXiv:2408.121 12179, 2024.
Yidong Chen, Chen Li, and Zhonghua Lu. Computing wasserstein-p distance between images with linear cost. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 519-528, 2022.

---

Published as a conference paper at ICLR 2026
Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems, 26, 2013.
Ayelet Dominitz and Allen Tannenbaum. Texture mapping via optimal mass transport. IEEE transactions on visualization and computer graphics, 16(3):419-433, 2009.
Pavel Dvurechensky, Alexander Gasnikov, and Alexey Kroshnin. Computational optimal transport: Complexity by accelerated gradient descent is better than by sinkhorn's 's algorithm. In International conference on machine learning, pp. 1367-1376. PMLR, 2018.
Ermie Esser, Xiaoqun Zhang, and Tony F Chan. A general framework for a class of first order primal- dual algorithms for convex optimization in imaging science. SIAM Journal on Imaging Sciences, 3(4):1015-1046, 2010.
Raphael A Finkel and Jon Louis Bentley. Quad trees a data structure for retrieval on composite keys. Acta informatica, 4(1):1-9, 1974.
Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z Alaya, Aurélie Boisbunon, Stanis- las Chambon, Lactitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, et al. Pot: Python optimal transport. Journal of Machine Learning Research, 22(78):1-8, 2021.
Harold N Gabow and Robert E Tarjan. Faster scaling algorithms for general graph matching prob- lems. Journal of the ACM (JACM), 38(4):815-833, 1991.
Samuel Gerber and Mauro Maggioni. Multiscale strategies for computing optimal transport. Journal of Machine Learning Research, 18(72):1-32, 2017.
Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2024. URL https: /wwww. gurobi com
Steven Haker, Lei Zhu, Allen Tannenbaum, and Sigurd Angenent. Optimal mass transport for reg- istration and warping. International Journal of computer vision, 60(3):225-22 2004.
Peter Halmos, Julian Gold, Xinhao Liu, and Benjamin Raphael. Hierarchical refinement: Optimal transport to infinity and beyond. In International conference machine learning, 2025.
Bishi He, Yuanjiao Chen, Darong Zhu, and Zhe Xu, Domain adaptation via wasserstein distance and discrepancy metric for chest x-ray image classification. Scientific Reports, 14(1):2690, 2024.
Ruiqiang He, Xiangchu Feng, Xiaolong Zhu, Hua Huang, and Bingzhe Wei. Rwrm: Residual wasserstein regularization model for image restoration. Inverse Problems & Imaging, 15(6), 2021.
Ka-Hei Hui, Chao Liu, Xiaohui Zeng, Chi-Wing Fu, and Arash Vahdat. Not-so-optimal transport flows for 3d point cloud generation. arXiv preprint arXiv:2502.12456 2025.
Arun Jambulapati, Aaron Sidford, and Kevin Tian. A direct tilde (O)(//esiion) iteration parallel algorithm for optimal transport. Advances in Neural Information Processing Systems, 32, 2019.
Tero Karras, Samuli Laine, and Timo Aila. A style-based generator architecture for generative adversarial networks. Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pp. 4401-4410, 2019.
Soheil Kolouri, Kimia Nadjahi, Umut Simsekli, Roland Badeau, and Gustavo Rohde, Generalized sliced wasserstein distances. Advances in neural information processing systems, 32, 2019.
Nikita Kornilov, Petr Mokrov, Alexander Gasnikov, and Aleksandr Korotin. Optimal flow match- ing: Learning straight trajectories in just one step. Advances in Neural Information Processing Systems, 37:104180-104204 2024.
Arthur Leclaire and Julien Rabin. A fast multi-layer approximation to semi-discrete optimal trans- port. In International Conference on Scale Space and Variational Methods in Computer Vision, pp. 341-353. Springer, 2019.
12

---

Published as a conference paper at ICLR 2026
Xudong Li, Defeng Sun, and Kim- Chuan Toh. An asymptotically superlinearly convergent semismooth newton augmented lagrangian method for linear programming. SIAM Journal on Optimization, 30(3):2410-244 14400 2020.
Yiyang Liu, Zaiwen Wen, and Wotao Yin. A multiscale semi-smooth newton method for optimal transport. Journal of Scientific Computing, 91(2):39, 2022.
Haihao Lu and Jinwen Yang. Pdot: A practical primal-dual algorithm and a gpu-based solver for optimal transport. arXiv preprint 1xXi:22407.1968) 2024.
Haihao Lu, Jinwen Yang, Haodong Hu, Qi Huangfu, Jinsong Liu, Tianhao Liu, Yinyu Ye, Chuwen Zhang, and Dongdong Ge. cupdlp-c: A strengthened implementation of cupdlp for linear pro- gramming by language. arXiv preprint MrXiv:2232.14882 2023.
Haihao Lu, Zedong Peng, and Jinwen Yang. cupdlpx: A further enhanced gpu-based first-order solver for linear programming. arXiv preprint arXiv:2507. 14051, 2025.
Zhongxuan Luo, Wei Chen, Na Lei, Yang Guo, Tong Zhao, Jiakun Liu, and Xianfeng Gu. The singularity set of optimal transportation maps. Computational Mathematics and Mathematical Physics, (2)3))1313-13333 2022.
Vien V Mai, Jacob Lindbäck, and Mikael Johansson. A fast and accurate splitting method for optimal transport: Analysis and implementation. arXiv preprint arXiv:2110.11738, 2021.
Quentin Mérigot. A multiscale approach to optimal transport. In Computer graphics forum, vol- ume 30, pp. 1583-1592. Wiley Online Library, 2011.
Eduardo Fernandes Montesuma and Fred Maurice Ngole Mboula. Wasserstein barycenter for multi- source domain adaptation. Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 16785-16793, 2021.
Kimia Nadjahi, Alain Durmus, Pierre E Jacob, Roland Badeau, and Umut Simsekli. Fast approxi- mation of the sliced wasserstein distance using concentration of random projections. Advances in Neural Information Processing Systems, 14:12411-12422 2021.
Ofir Pele and Michael Werman. Fast and robust earth mover's distances. In 2009 IEEE 12th international conference computer vision, pp. 460-467. IEEE, 2009.
François Pitié and Anil Kokaram. The linear konte-kantorovitch linear colour mapping for example-based colour transfer. In 4th European conference on visual media production, pp. 1-9. IET, 2007.
Thomas Pock and Antonin Chambolle. Diagonal preconditioning for first order primal-dual al- gorithms in convex optimization. In 2011 International Conference Computer Vision, pp. 1762-1769, IEEE, 2011.
Julien Rabin, Gabriel Peyré, Julie Delon, and Marc Bernot. Wasserstein barycenter and its appli- cation texture mixing. In International conference scale space and variational methods in computer vision, pp. 435-446. Springer, 2011.
Tauseef Rehman, G Pryor, and A Tannenbaum. Fast multigrid optimal mass transport for image registration and morphing. British Machine Vision Conference, 2007.
Meyer Scetbon, Marco Cuturi, and Gabriel Peyré. Low-rank sinkhorn factorization. In International Conference on Machine Learning, pp. 9344 9354. PMLR, 2021.
Bernhard Schmitzer. A sparse multiscale algorithm for dense optimal transport. Journal of Mathematical Imaging and Vision, 56(2):238-259, 2016.
Bernhard Schmitzer. Stabilized sparse scaling algorithms for entropy regularized transport probemsss SIAM Journal on Scientific Computing, 413))A144-A11481 2019.
Jörn Schrieber, Dominic Schuhmacher, and Carsten Gottschlich. Dotmark a benchmark for discrete optimal transport. IEEE Access, 5:271-282, 2016.
13

---

Published as a conference paper at ICLR 2026
Sameer Shirdhonkar and David W Jacobs. Approximate earth mover's distance in linear time. In 2008 IEEE Conference on Computer Vision and Pattern Recognition, pp. 1-8. IEEE, 2008.
Justin Solomon, Fernando De Goes, Gabriel Peyré, Marco Cuturi, Adrian Butscher, Andy Nguyen, Tao Du, and Leonidas Guibas. Convolutional wasserstein distances: Efficient optimal transporta- tion on geometric domains. ACM Transactions on Graphics (ToG), 34(4):1-11, 2015.
Wei Ren Tan, Chee Seng Chan, Hernan Aguirre, and Kiyoshi Tanaka. Improved artgan for condi- tional synthesis of natural image and artwork. IEEE Transactions on Image Processing, 28(1); 394-409, 2019. doi: 10.1109/T1P.2018 2866698. 2806761..21188866669888888888888
2018 2866688
Alexander Tong, Kilian Fatras, Nikolay Malkin, Guillaume Huguet, Yanlei Zhang, Jarrid Rector- Brooks, Guy Wolf, and Yoshua Bengio. Improving and generalizing flow-based generative models with minibatch optimal transport. arXiv preprint arXiv:2302.00482 2023.
Zhirong Wu, Shuran Song, Aditya Khosla, Fisher Yu, Linguang Zhang, Xiaoou Tang, and Jianxiong Xiao. 3d shapenets: A deep representation for volumetric shapes. Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 1912-1920, 2015.
Yue Xie, Zhongjian Wang, and Zhiwen Zhang. Randomized methods for computing optimal trans- port without regularization and their convergence analysis. Journal of Scientific Computing, 100 (2):37, 2024.
Guojun Zhang. Yancheng Yuan, and Defeng Sun. An efficient hpr algorithm for the wasserstein barycenter problem with o(dim(p)/s) computational complexity. arXiv preprint arXiv:2211.1488. 2022.
Guojun Zhang, Zhexuan Gu, Yancheng Yuan, and Defeng Sun. Hot: efficient halpern accel- crating algorithm for optimal transport problems. IEEE Transactions Pattern Analysis and Machine Intelligence, 2025.
A THE USE LARGE LANGUAGE MODELS (LLMs)
We used large language model solely to polish writing and to assist with small code snippets. It did not generate ideas, proofs, experimental designs, or results. All LLM-assisted content was authored, reviewed, and validated by the authors.
B PROOFS OF THEORETICAL RESULTS
B. PROOF OF THEOREM 1
Proof 1 (Proof of Theorem 1 Fix s E S and let Dx+1(s) C D denote the set of targets d that are
shielded from by the pairs ((s', ma)x(s')) s' E (()) Hence every (s, d) with d E Ds+1(s) must be explicitly added to Ne+1-
We first establish a geometric localization bound. Take d E Dx+1(s) and set = d - (s). By assumption 7 there exists s' E <s(s such that
(v, 5' - || - s|l
Insert and subtract mapy and apply assumption 4 together with Cauchy- to obtain
(v. s -) 8) = (d map + map (s) - map -)
2 I|d mapx (s)ll ||s - s|l ||map (s) (map -|s' -sl
> (q ||d mapp (s)]l - L||s' - =l)) I|s' =sll-
By assumption 2 we have I|s - =1l < D, hence
(d - (s), s'' x) 2 Ildd map (5)]l L.D) ||s =sll-
14

---

Published as a conference paper at ICLR 2026
LD/q, then (d (s) s) 0 for some s' squared Euclidean ground cost this is exactly the sufficient shielding condition stated in Section 5.2 of Schmitterr (2016). Therefore each unshielded d must satisfy
If l|d ( (s)ll
S(s). For the
LD_ || - mapa (+)]l S := 4
that is, (8+) C (R(mapa(>))
Assumption
5 gives (map;(s) (s))( & Ro for all ks which implies
(8)1(8) C BRo+ for all 0. Taking the union over all iterations,
U D C Dn BR, (m)p R20
By assumption 3 the cardinality is uniformly bounded as
U ( 120
[D BRo ((map,())) - pvoln (Bo(+BR+))) = Co.
which depends only the constants in assumptions 5 the dimension of the space, and is independent of the sizes of and Dr.
We now conclude the proof. Since each has at most Co unshielded targets, if at iteration k there is still some d (s)) then least one pair (s,d) d) must be added. As the total number of such unshielded pairs is bounded by Co. the process of adding new shielding edges can occur most Co times. Thus the sequence (Na) converges after at most Co iterations, which proves ccale-independent convergence.
Regarding sparsity, the shielding part contributes most Cols] pairs. Step 3 of Algorithmm can add at most K pairs per iteration, and since there are at most Co iterations, its total contribution is bounded by KCO. Therefore the final support satisfies
|N&l & Ca/S1 + BCoB3 (1 B)Cu|S)
Thus we can take := = (1 + B)Co, which depends only on the constants in assumptions on the space dimension, and on the chosen parameter B. At the fixed point 3y locally optimal with respect to Nr, and by construction Ng. is shielding neighbourhood. By the local-to-globe certification (Corollary 3.10 Schmitzer (2066)) a globally optimal once the process stabilizes. This completes the proof.
7
B.2 PROOF OF PROPOSITION
Proof 2 We first compute BBT;
diag( (r)
BBT
P
diag(c)
where E (0,1)mxn is the adjacency matrix of bipartite graph representing the non-zero entries in B. Next, we compute the rescaled matrix BB
BBT = D )
Im
In Q = D, 1/2 122//D 1/2
Q
= 2 Q
1
where D, = drag(r) and D. = diag(c).
I Q
For the matrix Q, it is known that the eigenvalues of the block matrix
are (Q))),
Q
thus:
1 7 + max (Q)).
1 + Tmax (Q)
Xmax. [BB
BB
2
We prove that Tmaxe (Q)
15

---

Published as a conference paper at ICLR 2026
Lower bound: Let u (yri..... (m)) and (vGi..... ........ vCn)
T. We have:
(Qv); = Σ
vg
Ti
vri Us
JEN()) VTici
Similarly, Quu u = v, so u and v are singular vectors, and the singular value is exactly Thus, mmax Q) 2 1.
Upper bound: For any a E IR", by the Cauchy-Schware inequality:
[Qz = -ー
BJ
Σ
ΣΣ Σ ag.
VS
(EN())
(EN()
JEN((
JEN(S)
Here, (i) denotes the indices the the non-zero elements in the adjacency matrix P. Thus, | S Therefore, max (Q) = 1,, and we conclude:
1+1
|
2
=1
2
C EXTENSION TO NON-GRID DATA
C.I CONSTRUCTION OF HIERARCHICAL STRUCTURE
For non-grid data (e.g., p clouds), primarily employ 2r-tree structure (e.g., Quadtrees in 2D, Octrees 3D) to construct the hierarchy via spatial partitioning. The procedure starts by defining an axis-aligned bounding hypercube containing all data points, which serves as the coarsest level. Finer levels are generated by subdividing each hypercube into 2 equal-sized sub-cubes, discarding any sub-cubes that contain no data. This subdivision continues recursively until a pre-determined depth is reached or the finest level (containing individual points) is achieved. At each level, the active nodes correspond these non non-empty sub-cubes, and the representative point is defined as the geometric center of the spatial region.
In higher dimensions, standard 2d-trees become intractable due to the 2d branching factor. To ad- dress this, sequential axis partitioning (similar to k-d trees) can be adopted to control the inter-level reduction ratio, ensuring hierarchy construction does not become computational bottleneck.
C.2 SHIELDING STRATEGY ON NON-GRID DATA
In the shielding-based active support update (Algorithm 3) the definitions of the local neighborhood
Rs)) and the unshielded D(s) require adaptation for non-grid domains. In the image setting (2D grids), R(s) consists of the 8 surrounding pixels. To maintain consistency and ensure extensibility high dimensions, we replace this with K-Nearest Neighbors (KNN), setting km 4 x d (e.g..
km = 12 for 3D point clouds). This choice aligns with the image setting (where 8 = 4 x 2) and keeps the neighborhood sparse.
Unlike grids where the spatial distribution is uniform, non-grid data exhibit irregular distributions. Consequently, single source point may correspond to a large number of unshielded candidates. For computational and memory efficiency, we enforce upper bound on the number of candidates added per iteration, denoted Umax-
Crucially, this budgeted strategy affects neither the convergence of HALO nor the theoretical result of Theorem [ Regarding convergence, hierarchical algorithms rely on the monotonic decrease of the objective function, which is guaranteed as long as N++1 5 supp(æx).. Regarding scale- independence, the complexity bound in Theorem D is derived from the total volume of potential unshielded targets (bounded by the constant Co in the proof B.I). While the theorem assumes all unshielded targets are identified, limiting the update size to Umax merely implies that these necessary neighbors are added over slightly more iterations, preserving the ecale-independen complexity.
Finally, we provide sensitivity analysis in Table and Table demonstrate the stability of our al- gorithm with respect to this hyperparameter. We tested Umax (10, 15, 20, 30, 40) ModelNet10
16

---

Published as a conference paper at ICLR 2026
(see Appendix The results indicate that the performance remains stable across varying scales, and we use U mmax 20 as a safe default setting, which is used throughout non non-grid experiments.
Table 8: Time (8) ablation on the budget size
Table 9: size U max for shielding on 3D point clouds.
Memory (GB) ablation on the budget
U Umax for shielding on 3D point clouds.
U max
2
2
26
2º7
2
max
2
2
2
2
2
16.18
28.32
50.62
92.82
265.8
10
0.54
0.62
0.74
1.06
1.72
15
15.70
26.83
48.92
88.15
236.6
15
0.54
0.63
0.76
1.10
1.76
15.54
26.25
47.42
88.51
229.7
20
0.54
0.64
0.77
1.12
1.83
30
15.85
27.58
50.18
91.30
224.7
30
0.55
0.65
0.79
1.17
1.92
40
16.04
26.09
51.81
103.7
244.2
40
0.56
0.61
0.83
1.21
1.95
C.3
NUMERICAL RESULTS ON 2D NON-GRID DATA
In addition to the 3D experiments presented Table 7 we provide results on 2D non-grid data constructed from ModelNet1O (see Appendix [ details). The results summarized in Table 10 Similar to the 3D setting, HALO demonstrates consistent superior performance. At n 21 it achieves 3.47x speedup, an 82.4% reduction in memory usage compared to HiRef, f, and 32.5% lower transport cost.
Table 10: Performance on 2D Non-Grid Data (ModelNe:1O-PCAA) gap denotes the relative objec- tive difference: for 214, the reference the exact solution computed by the standard EMD) solver Flamary al. (2021): for 2 where EMD solver is intractable, the reference is the solution of HALO.
Metric
Method
214
215
216
217
218
219
HALO
7.82
12.50
23.44
53.25
121.6
247.0
time
HiRef
23.40
46.40
94.10
190.4
422.0
OOM
Sinkhorn
28.50
OOM
OOM
OOM
OOM
OOM
HALO
0.55
0.62
0.74
1.08
1.92
3.97
memory
HiRef
0.90
1.02
1.65
3.60
10.89
OOM
Sinkhorn
10.60
OOM
OOM
OOM
OOM
OOM
HALO
+5.68E
gap
HiRef
+4.37E
+3.82E OOM
+4.04E-1
+3.00E-1
+3.25E-1
OOM
Sinkhorn
+4.02E
OOM
0OM
OOM
OOM
D ALGORITHM SETTINGS
For HALO, we employ the cuPDLPx Lu et al. (2025) solver with constant step-sizes unless oth- erwise specified. All experiments use the Pock Chambolle rescaling scheme, and the stopping cri- terion is set uniformly with primal fcasibility, dual feasibility, and objective gap thresholds of 10 In updateActive,, we choose = 0.25 S) for the operator Topx
For HOT Zhang et al. (2025). adopt the open-souree implementation with its default parameter
choices. The stopping criterion is fixed to 10-6 to ensure comparability with other baselines.
For ShortCut [Schmitzer (SchmittttSchmiteerr (2016), use the variant based the LEMON solver provided in the released code, with all default parameters left unchanged.
For M3S Chen et al. 2022 which is entropic regularization method, the official im- plementation with its predefined entropy-regularizztion coefficient and all other default parameter settings.
For Gurobi Gurobi Optimization, LLC 2024). we rely on the Barrier algorithm. Unless otherwise specified in the main text, the crossover procedure is disabled.
17

---

Published as a conference paper at ICLR 2026
For HiRef Halmos al. (225)) adopt the official open-source implementation and strictly adhere to the default parameter settings provided by the authors.
For Sinkhorn, employ the standard implementation provided by the POT library Flamary ct al. 2021). We set the regularization parameter to e = 10 to obtain an approximate solution with reasonable accuracy.
E DATASET
To comprchensively evaluate the scalability and generalization of HALO, we employed two distinct datasets covering both grid-based images and unstructured point clouds.
DOTmark (Image Data). To evaluate scalability across different problem sizes on grid data, we utilized the DOTmark (Schrieeer al. (2016) benchmark. In addition to the native resolutions (32 to 512), constructed additional variants by both downsampling and upsampling. For downsampling, each native-resolution image was resized to 16 16 using bilinear interpolation. The interpolated values were rescaled linearly to match the original intensity range, rounded integers, and clipped to avoid numerical overflow. For upsampling, we generated 1024 x 1024 images by bilinear interpo- lation, followed by rounding to integers. To guarantee exact consistency with the original data, each pixel at the original grid was enforced to coincide with its corresponding position in the enlarged image, and the resulting values were clipped the original intensity range before being stored. In this way, the 16 x 16 images provide small-scale test cases, while the 1024 x 1024 images serve as challenging large-scale benchmarks.
ModelNet10 (Non-Grid Data). To validate the performance of HALO on non-grid data, we con- structed benchmark using ModeINet10 Wu et al. (2015), widely used dataset for 3D point cloud analysis. We selected the top-3 samples from each of the 10 classes, generating a total of 30 pairs of point clouds for evaluation. The raw point clouds were normalized the unit hypercube, and varied the number of points n from 2 to 219 via random sampling to test scalability. To further evaluate the algorithm non-grid 2D data, we generated 2D counterparts of these 3D shapes using Principal Component Analysis (PCA), creating a non-grid 2D point cloud benchmark distinct from the regular grids DOTmark.
Qemcty
Classccmages
GiRFoiggh
LogGRF
LegitGRF
Microsorpy
Shapes
Whisenocse
Figure 4: Example images from the DOTmark benchmark.
F SENSITIVITY ANALYSIS OF HYPERPARAMETER IN DUAL- VIOLATION CORRECTION
key hyperparameter in HALO is 3, associated with the dual-violation correction step in Algo- rithm It controls the size of the candidate set added during the update by selecting the top K S pairs with the largest dual violations. To assess the sensitivity of HALO to this pa- rameter, evaluated the algorithm the DOTmark benchmark across wide range of values:
18

---

Published as a conference paper at ICLR 2026
Bathnab
Bed
Chair
Desk
Dreser
Night
Figure 5: Example 3D point clouds from the ModelNet10 benchmark. Each object is sampled with 8192 points.
E (2-4.2-3,22223 20). We also included the baseline case B 0, which actually disables the dual-violation correction module.
The numerical results are detailed in Table I1 The algorithm exhibits high stability regarding both memory and runtime: GPU memory usage remains virtually unaffected by the choice of 8; runtime also shows minimal fluctuation, with a maximum variation of approximately 17% at resolution 1024 x 1024. However, setting / 0 leads to sharp performance drop, which confirms that the violavion correction is crucial for the efficiency of the solver.
Based on these findings, we recommend 8 2- as safe default setting, which is used throughout our main experiments.
Table 11: Sensitivity analysis of 3 on DOTmark. Performance remains stable for E [2-3 20).
Setting 8 = 0 leads to significant slowdowns. Memory usage is consistent across all non-zero values.
8
y- r= 64
128
r
256
r
512
1024 r
20
0.77
1.18
2.79
9.24
28.86
2-1
0.78
1.12
2.70
8.76
25.19
2 2
0.78
1.11
2.62
8.84
25.42
time
2
0.72
1.11
2.71
8.72
24.67
2-4
1.73
2.42
4.44
11.90
29.64
0 (No dual)
1.90
2.99
5.47
15.23
52.85
memory
All
~ 0.39
0.48
0.76
2.10
~ 6.30
G COMPARISON WITH A GPU-BASED SHORTCUT IMPLEMENTATION
natural question arises regarding whether the ShortCut method could yield even better results than HALO if implemented on GPUs. To clarify this, we implemented ShortCut-CPP variant within the same framework as HALC, utilizing cuPDLPx x as the underlying solver. The primary difference lies in the active-support update strategy: ShortCut employs an aggressive pruning strategy that retains only the support the current coupling. whereas HALO uses conservative update (N++1 > N₂)) augmented with dual-violation conrection.
The comparison results on DOTmark are presented in Table While the GPU implementation brings efficiency gains ShortCut-GPU compared to CPU baselines, HALO still outperforms ShortCut-Coe -CPU significantly, achieving a 2.5 speedup at resolution 1024 x 1024.
This performance gap highlights critical algorithmic contribution of HALO. GPU-based first-order solvers typically yield solutions with lower precision compared to CPU-based classical methods. ShortCut' 's aggressive active-support update is highly sensitive this numerical noise, which
19

---

Published as a conference paper at ICLR 2026
leads to stagnation during the refinement process. In contrast, HALO's conservative update rule and the dual-violation correction provide the necessary stability. These designs make the hierarchical framework robust to the lower precision of GPU solvers, thereby unlocking the full potential of GPU acceleration.
Table 12: Runtime comparison between ShortCut-GPO and HALO DOTmark. Time is in seconds (s).
Resolution
512
1024
ShortCut-GPJ
21.40
68.33
HALO
11.17
27.73
H FLEXIBILITY WITH ALTERNATIVE SOLVERS
While HALO utilizes cuPDLPX Lu et al. (0055) as the default LP solver due to its state-of-the- art performance, the proposed hierarchical framework designed to be flexible and compatible with various first-order GPU solvers. To validate this flexibility, integrated an alternative solver, HPR-LP Chen et al. (024),, into HALO by replacing the backend of the solveRestricte component. The numerical results DOTmark are reported in Table I3
Although HPR-LP generally slower than cuPDLPx in this context, the HALO framework still effectively leverages it to solve large-scale instances efficiently. Specifically, resolution 1024 x 1024, HALO integrated with HPR-LP achieves a 6.1 speedup and a 68.8% reduction memory usage compared the state-of-the-art solver M3S (see Table
This experiment confirms that the efficiency of HALO is not solely dependent on a specific un- derlying LP solver; rather, the hierarchical active-support framework is critical component that significantly contributes to the overall performance. is worth noting that GPU-based LP solvers are a burgeoning field compared to mature CPU solvers. We believe that the continuous evolution of faster GPU-based solvers will further boost HALO's performance, solidifying our framework as highly promising direction for solving large-scale OT.
Table 13: Performance of HALO integrated with the alternative HPR-LP solver Chen al. 2024) on DOTmark. Time is in seconds (s) and Memory is in gigabytes (GB).
Resolution
64
128
256
512
1024
time
1.06
1.70
3.31
9.68
40.76
memory
0.72
0.78
1.00
1.90
6.61
gap
6.39E 6
2.36E-5
1.76E
infeas
4.09E-1
2.53E
1.33E-7
1.03E
6.08E-8
I ROBUSTNESS AND GENERALIZATION ACROSS DIVERSE DATASETS
To demonstrate generalization beyond DOTmark, we conducted additional evaluations on two high- resolution datasets: -ackk-Faces-HQ (FFHQ) Karras et al. (219)) and WikiArt Tan al. (019). Specifically, we selected the first 10 images from FFHQ and the first image from each of the top- 10 classes in WikiArt. For each dataset, we generated 45 instances by pairing these images, strictly following the DOTmark preprocessing pipeline.
The numerical results summarized in Table 14 HALO consistently achieves low memory usage and fast solving speeds on these diverse tasks, exhibiting performance metrics highly consistent with those on mOTmark. This confirms the broad applicability and robustness of HALO across diverse real-world datasets.
20

---

Published as a conference paper at ICLR 2026
Table 14: Generalization performance of HALO real-world datasets (FFHQ and WikiArt) com- pared to the standard DOTmark. Results are averaged over instances. Time is in seconds (s) and Memory is gigabytes (GB).
Resolution
256
512
1024
FFHQ
3.24
10.79
23.68
time
WikiArt
2.74
9.13
22.14
DOTmark
4.31
11.17
27.73
FFHQ
0.77
2.19
6.52
memory
WikiArt
0.74
1.97
6.10
DOTmark
0.76
2.07
6.25
21


Reviews:
[{"review_id": "0U7unVmYP1", "reviewer": "Reviewer_rnK4", "summary": "- The authors proposes HALO (Hierarchical Algorithm for Large-scale Optimal Transport), a novel GPU-friendly method for solving large discrete optimal transport (OT) problems efficiently. \n- The key idea is to build a coarse-to-fine multi-scale hierarchy and iteratively refine the transport plan while maintaining an active support set that captures potentially non-zero couplings. \n- Each refinement step solves a restricted OT problem via a primal-dual hybrid gradient (PDHG) method optimized for GPU computation.\n- Experiments on large-scale image OT tasks (DOTmark dataset) demonstrate significant speed and memory savings ove baselines like HOT, ShortCut, and M3S.", "questions": "Please mainly see the above weakness part", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- The paper successfully solves the scalability bottleneck of OT on high-resolution data with well-designed hierarchical + active-support framework and solid theoretical justification.\n- The experimental results are promising with great speed and memory advantage.", "weaknesses": "- Currently limited to 2D grid supports with squared Euclidean cost.\n- The active-support update relies on heuristic parameters (e.g., β) without sensitivity analysis.\n-Experiments only cover image data; no test on non-grid or higher-dimensional problems. Does it work on other distributions?", "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": "PQYhkxCaN7", "reviewer": "Reviewer_u6SB", "summary": "The paper introduces HALO (Hierarchical Algorithm for Large-scale Optimal Transport), a GPU-based, memory-efficient method for solving large-scale optimal transport (OT) problems in the plane with squared Euclidean cost. HALO combines a multiscale hierarchical framework with a sparse, active-support refinement scheme and a primal–dual hybrid gradient (PDHG) solver to overcome the severe memory and scalability limitations of existing OT solvers. By solving coarser OT problems first and using their solutions to warm-start finer levels, HALO efficiently refines the transport plan while keeping memory usage to  𝑂(r^2), where r is the number of pixels per dimension.\n\nThe method also includes a dual-violation augmentation step that improves robustness, and the authors establish a scale-independent iteration bound, proving that each refinement level requires only 𝑂(1) iterations. On the DOTmark benchmark, HALO outperforms state-of-the-art methods such as HOT, ShortCut, and M3S—achieving up to 8.9× speedup and 70% lower GPU memory use for 1024×1024 images, while maintaining comparable or superior accuracy. Overall, HALO demonstrates near-linear runtime scaling and high parallel efficiency, offering a theoretically grounded and practically scalable solution for large-scale OT computation.", "questions": "The following questions need to be addressed to further improve the quality:\n\n1. transportation cost: if the transportation cost is not the squared Euclidean distance, like L1 distance, will the method be applicable ?\n2. dimension : if the problem is not restricted on 2d, but general n dimensional Euclidean space, can the method work ?\n3. range: if the support of the target measure is not convex, but a concave Jordan domain, can the active-support refinement  algorithm still work ? In this situation, there will be complicated singularities in the domain, determining the singularity will be challenging.\n4. HALO is built on a first -order optimization framework, is it possible to use Newton's method ?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper’s main strengths lie in its combination of theoretical rigor and practical scalability. HALO introduces a hierarchical, GPU-friendly framework that reduces the memory requirement of large-scale optimal transport (OT) from 𝑂(𝑟^4) to 𝑂(𝑟^2), allowing it to handle very high-resolution problems that were previously infeasible. Its coarse-to-fine multiscale design and sparsity-based active support updates enable near-linear runtime scaling, while the factorization-free PDHG solver fully exploits GPU parallelism for efficient computation.\n\nEqually important, the paper provides strong theoretical guarantees and robust empirical validation. The authors prove a scale-independent iteration bound, ensuring that refinement at each level converges in a constant number of steps, and introduce dual-violation augmentation to improve robustness without sacrificing sparsity. Extensive experiments on the DOTmark benchmark confirm that HALO achieves up to 8.9× speedup and 70% memory savings compared to leading solvers, while maintaining high accuracy. Overall, the method is both theoretically elegant and practically impactful, setting a new benchmark for scalable OT computation.", "weaknesses": "The main weaknesses of the paper lie in its limited generality and empirical scope. HALO is tailored for 2D optimal transport problems with squared Euclidean cost on regular grids, and both its theoretical guarantees and hierarchical design rely on this structure. As a result, the method’s applicability to higher-dimensional settings, irregular domains, or non-Euclidean costs remains unclear. Moreover, while the algorithm performs impressively on the DOTmark benchmark, its evaluation is confined to this dataset, leaving open questions about robustness and generalization to more diverse or real-world applications.\n\nAnother limitation is the dependence on heuristic components such as dual-violation augmentation and Top-K active-set selection, whose performance may vary with parameter choices not deeply analyzed in the paper. The implementation is also technically complex—combining multiscale hierarchy, active-support refinement, and GPU-based PDHG—which could make reproduction or extension challenging for practitioners. Overall, while HALO is methodologically strong and achieves excellent performance, its restricted scope, heuristic tuning, and limited experimental diversity temper its broader applicability.", "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": "tYB99mmhbp", "reviewer": "Reviewer_a2S6", "summary": "The paper proposed a hierarchical algorithm to solve large-scare OT problems. By several techniques including hierarchical expansion, active support update, and GPU based PDHG. The algorithm justifies its efficiency by improving the state-of-the-art baselines and theoretical scale independent iteration complexity bound.", "questions": "Beyond the weakness that requires the author to address.\n1. The active support in Def. 1 also appears in [1] and even earlier (see references in [1]), will such randomized method improve the HALO?\n2. In experiment, the ShortCut is only implemented on CPU, given its low memory cost and rather low runtime in lower resolution, is there theoretical barrier for ShortCut to implement on GPU and yield an even better result?\n\n[1] Xie, Yue, Zhongjian Wang, and Zhiwen Zhang. \"Randomized methods for computing optimal transport without regularization and their convergence analysis.\" Journal of Scientific Computing 100.2 (2024): 37.", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "By carefully tuning the three approaches in discrete OT techniques, the algorithm justify itself by testing on common data sets.", "weaknesses": "1. A well known issue of the hierarchical approach, is the scalability to higher dimensions rather than 2D mesh like topology. The author shall consider a least a remark on extension to higher dimensions.\n2. Another issue of the manuscript is the theoretical justification relies on rather strong assumption (4,5), also as the authors mentioned. Since the complexity of HALO is only an upper bound under strong assumptions, the practical superiority of the HALO is not well-explained. The author may need further discuss the structural advantages comparing with other algorithm, for instance Multiscale-OT.", "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": "YaMOrH01Ct", "reviewer": "Reviewer_Tm5V", "summary": "This paper introduces HALO, a hierarchical and memory-efficient algorithm for large-scale optimal\ntransport (OT) problems on 2D supports with squared Euclidean cost. The method combines a coarseto-\nfine multi-scale framework with a GPU-optimized PDHG (Primal-Dual Hybrid Gradient) solver\nand an active-pruning mechanism for sparsity. The authors prove a scale-independent iteration\ncomplexity bound for the refinement phase and show strong empirical results on the DOTmark\nbenchmark, achieving significant speedup and memory usage saving compared to methods such as\nHOT, ShortCut, and M3S.", "questions": "1. HALO is currently limited to 2D support with squared Euclidean cost. The hierarchical design\nheavily relies on the specific problem (OT between two 2D images), while applications of OT\n(like in generative models) are on more general settings (Wasserstein GAN [1] and optimal flow\nmatching [2]) where the supports are high-dimensional. In this case, the HALO may not work well\nas in the image problems in DOTmark because the data locality property doesn't hold in these\nproblem. The limited problem setting is mentioned but not stated clearly in the introduction.\n\n2. Although the paper fairly acknowledges this, HALO’s performance heavily relies on a third-party\nGPU solver. A deeper analysis of how HALO interacts with other first-order solvers would make\nthe contribution more self-contained.\n\n3. The parameter beta in the dual-violation step appears to influence sparsity and runtime, but the\nsensitivity analyses are missing. Reporting how beta impacts convergence and memory would be\nvaluable.\n\n4. In DOTMark [3], there are various types of images. It would be interesting to see how HALO\nperforms on these classes of images respectively to have a clearer idea of the practicality of the\nassumptions and the influence of the data locality on the convergence of HALO.\n\nReferences:\n[1] Adler, Jonas, and Sebastian Lunz. \"Banach wasserstein gan.\" Advances in neural information\nprocessing systems 31 (2018).\n[2] Kornilov, Nikita, et al. \"Optimal flow matching: Learning straight trajectories in just one\nstep.\" Advances in Neural Information Processing Systems 37 (2024): 104180-104204.\n[3] Schrieber, Jörn, Dominic Schuhmacher, and Carsten Gottschlich. \"Dotmark–a benchmark for\ndiscrete optimal transport.\" IEEE Access 5 (2016): 271-282.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. Theorem 1 (scale-independent iteration complexity) provides a nontrivial and interpretable\nconvergence guarantee, addressing the gap in multi-scale OT methods (which often lack formal\ncomplexity bounds).\n2. HALO merges the hierarchical multi-scale structure with a GPU-based LP solver. This is a\nthoughtful contribution, given the growing interest in efficient OT solvers.\n3. The paper is clearly written and well-organized. Figures 1–3 effectively illustrate the hierarchical\nprocess and empirical trends. Proofs are included in the appendix and seem technically sound.", "weaknesses": "1. HALO is currently limited to 2D support with squared Euclidean cost. The hierarchical design\nheavily relies on the specific problem (OT between two 2D images), while applications of OT\n(like in generative models) are on more general settings (Wasserstein GAN [1] and optimal flow\nmatching [2]) where the supports are high-dimensional. In this case, the HALO may not work well\nas in the image problems in DOTmark because the data locality property doesn't hold in these\nproblem. The limited problem setting is mentioned but not stated clearly in the introduction.\n\n2. Although the paper fairly acknowledges this, HALO’s performance heavily relies on a third-party\nGPU solver. A deeper analysis of how HALO interacts with other first-order solvers would make\nthe contribution more self-contained.\n\n3. The parameter beta in the dual-violation step appears to influence sparsity and runtime, but the\nsensitivity analyses are missing. Reporting how beta impacts convergence and memory would be\nvaluable.\n\n4. In DOTMark [3], there are various types of images. It would be interesting to see how HALO\nperforms on these classes of images respectively to have a clearer idea of the practicality of the\nassumptions and the influence of the data locality on the convergence of HALO.\n\nReferences:\n[1] Adler, Jonas, and Sebastian Lunz. \"Banach wasserstein gan.\" Advances in neural information\nprocessing systems 31 (2018).\n[2] Kornilov, Nikita, et al. \"Optimal flow matching: Learning straight trajectories in just one\nstep.\" Advances in Neural Information Processing Systems 37 (2024): 104180-104204.\n[3] Schrieber, Jörn, Dominic Schuhmacher, and Carsten Gottschlich. \"Dotmark–a benchmark for\ndiscrete optimal transport.\" IEEE Access 5 (2016): 271-282.", "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)