{"id":"df2c07a2-d81d-499b-971a-5e7025d8b1e1","arxiv_id":"2608.12060","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":5,"one_line_summary":"MSIPRM integrates a multiscale OT hierarchy with an interior-point relaxation method, solving large OT problems through sparse adaptive active sets.","lead":"This paper combines a multiscale optimal transport hierarchy with a warm-startable interior-point relaxation solver. The resulting method solves sparse subproblems instead of the full transport plan, enabling problems whose full formulations have trillions of variables.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Global KKT residual at the trillion-variable scale is not certified: the paper gives no procedure for evaluating dual infeasibility over all O(mn) pairs without enumerating them.","rationale":"The paper is a serious algorithmic contribution with substantial theory, including nondegeneracy analysis, condition number estimates, and convergence results, and I do not dispute those components. My concern targets the practical certificate behind the headline empirical claim: the global KKT residual used for termination is over all transport variables, but the manuscript gives no method to compute or certify that quantity at the trillion-variable scale. This is not an ad hominem or a consensus disagreement; it is an internal gap between the stated stopping criterion and the described implementation. The reader's support-inheritance concern is real but concerns efficiency and graceful degradation; my concern is more fundamental because if the dual infeasibility check is only local, the reported residuals do not establish that the returned plan solves the full OT problem. I keep the CONDITIONAL verdict because the issue is potentially fixable by specifying and demonstrating an exact global screening procedure (or by weakening the claim to 'locally verified active-set solution') and by releasing code. If the proposed test shows that the true global dual residual exceeds the reported tolerance, the central claim should be rejected or substantially restated.","tokens_in":37005,"tokens_out":9814,"duration_ms":111058,"concrete_test":"Run MSIPRM on a 128×128 DOTmark image pair with tolout=1e-6. After termination, independently compute the exact global dual infeasibility δ = max_{i,j}(α_i+β_j−C_ij) over all 2.7e8 pairs, feasible in compiled code (or with an exact 2D additively weighted nearest-neighbor data structure for the squared-Euclidean cost). Compare δ with the reported ||(c−A^Tλ)_−||/(1+||c||). If δ > 1e-6 (or greater than 1e-5 allowing for rounding), the termination residual is local rather than global, and the active-set solution may be suboptimal. Repeat at 256×256 if feasible, and require the authors to release code so the same check can be performed on the claimed 2048×2048 instance.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Algorithm 4 terminates on KKTres,out (Section 6.1), whose dual-infeasibility term is ||(c−A^Tλ)_−||/(1+||c||) over all m·n transport variables. At 2048×2048, m·n≈1.8e13, so a single exact evaluation of c_ij−α_i−β_j for every pair is impossible on the reported hardware. Section 3.2 only says cost entries are 'generated on demand' to detect violated dual constraints, but the manuscript never specifies a certified procedure for finding the global maximum violation, e.g. an exact 2D additively weighted nearest-neighbor search exploiting C_ij=||u_i−v_j||^2. If the dual residual is checked only on the current active set or a geometrically local candidate set, then the reported KKT residual 7.5e-7 (Table 5) does not certify optimality of the full OT LP; a violated reduced cost outside the screened set could remain. This is load-bearing for the central claim because Proposition 2 requires knowledge of T_rel/T_abs over all pairs to guarantee finite termination, and the 'trillion variables solved' claim needs that global certificate. The support-inheritance concern raised by the reader is a performance risk; this is a correctness and verification gap.","agreement_with_reader":"disagree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes MSIPRM, a multiscale primal-dual interior-point relaxation method for discrete optimal transport. An outer multiscale hierarchy constructs coarse-to-fine OT problems, and at each level the solver restricts optimization to an adaptively refined active set initialized from the previous level's support. The inner solver, IPRM, is the authors' earlier barrier augmented Lagrangian method, which allows non-interior warm starts. The paper derives preprocessing and decomposition routines for sparse subproblems, proves generic primal constraint nondegeneracy, gives support-identification guarantees based on perturbed LP stability, analyzes Schur complement conditioning, and proves global and local convergence results. Numerical experiments compare MSIPRM against HOT, a CPLEX-based multiscale network simplex, Sp-Sinkhorn, and LEMON's network simplex on DOTmark and on synthetic 1024x1024 and 2048x2048 instances, reporting that the full 2048x2048 formulation contains about 1.8e13 variables while MSIPRM keeps about 7.1e7 active variables and reports a KKT residual of 7.5e-7.","tokens_in":37277,"tokens_out":10691,"duration_ms":110958,"significance":"If the reported numerical claims are certified, the paper would be a meaningful advance in practical large-scale discrete OT: it demonstrates that a single-machine solver can handle formulations with trillions of variables, and it outperforms several strong baselines on the tested benchmarks. The theoretical sections also contain solid contributions: Theorem 1 gives a measure-theoretic guarantee of generic primal constraint nondegeneracy, Lemma 4 and Theorem 2 provide support-stability and support-identification results, Proposition 4 gives explicit Schur complement condition-number bounds, and Theorem 4 establishes a global error bound. The proofs in the appendices are detailed and generally self-contained. The main caveats are that the global KKT residual reported at the largest scales is not certified by the manuscript, and the local quadratic convergence theorem assumes exact Newton solves while the implementation uses PCG; both issues are fixable but currently block full verification of the central claims.","major_comments":[{"comment":"The manuscript does not specify a certified procedure for evaluating the global dual-infeasibility term ||(c−Ā^Tλ)_−||/(1+||c||) over all m·n variables. Section 3.2 states only that cost entries are 'generated on demand' to detect violated dual constraints, and Section 6.4 mentions the overhead of evaluating c−Ā^Tλ, but no algorithm is provided for finding the global maximum violation over all pairs, e.g., an exact additively weighted nearest-neighbor query exploiting the separability of C_ij=||u_i−v_j||^2. Consequently, the KKTres,out values in Tables 4–5, including 7.5e-7 at 2048×2048, are not certified in the manuscript as residuals of the full LP; the paper does not establish that they were evaluated globally rather than only on the active set or a screened subset. This gap also affects Proposition 2, whose finite-termination proof requires the sets T_rel and T_abs in Algorithm 2 to be known over all pairs. Please describe the exact global screening procedure used, or explicitly restrict the optimality claims to certified active-set residuals.","section":"§3.2, §6.1 (Algorithm 2 and KKTres,out)"},{"comment":"Theorem 5 asserts local quadratic convergence for Algorithm 3, but the proof assumes exact Newton directions. In Section 4.4, however, the Schur complement system (29) is solved either by PCG with preconditioner U or by a direct Cholesky factorization. The PCG option is an inexact linear solve, and the paper gives no stopping criterion or inexact-Newton analysis that would make Theorem 5 applicable to the implemented algorithm. Please either add an inexact-Newton variant of Theorem 5 with explicit accuracy requirements for the PCG iterates, or state clearly that the quadratic-convergence theorem applies only to the direct-factorization variant.","section":"§4.4 and §5 (Theorem 5)"},{"comment":"The scalability of MSIPRM, and in particular the trillion-variable demonstration in Section 6.4, rests on the support-inheritance assumption encoded in Eq. (15): the fine-scale optimal plan is assumed to be supported essentially on the union of child cells of the coarse-scale support. The paper provides empirical evidence and single-level support-identification results (Lemma 4 and Theorem 2), but no theorem shows that the supports of the exact optimal plans at successive levels are nested under the random cost perturbations in Eq. (14). If this assumption fails, Algorithm 2 would need to introduce many new variables, and the memory and time advantages would shrink. I regard this as a performance risk rather than a correctness error in the reported experiments, but the paper should state it as an explicit assumption and discuss its limits.","section":"§4.1, Eq. (15)"}],"minor_comments":[{"comment":"Please specify how ||c|| is computed in the definition of KKTres,out when the full cost matrix is not stored; the scaling of this term affects the interpretation of the reported residuals.","section":"§6.1"},{"comment":"The K–W distance for the same pair at 1024×1024 is 121.61 in Table 4 and 243.13 in Table 5; the footnote explaining the factor-of-two scaling convention should appear at the first occurrence, and the convention should be stated in the main text.","section":"Tables 4 and 5"},{"comment":"The bound on ||A_P^{-1}A_Q|| would be easier to check if the paper explicitly stated that A_P^{-1} is integral because A_P is a nonsingular submatrix of a totally unimodular matrix.","section":"Appendix C, proof of Lemma 4"},{"comment":"In the rows for Test 1–Test 5 and Test 2–Test 4, the entries '31740.508.3E-07' and '15250.208.4E-07' are missing separators and should be corrected.","section":"Table 4"}],"recommendation":"major_revision","confidential_remarks":"The paper fits the journal's scope and the algorithmic idea is promising. My recommendation is driven by the certification gap in the global dual residual and the mismatch between Theorem 5 and the PCG implementation; both are fixable within the manuscript's scope. The support-inheritance issue is a performance risk that should be clearly labeled. I do not have concerns about novelty or attribution: the method is a natural integration of the authors' prior IPRM with a multiscale active-set framework, and the theoretical results on the OT side are new."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Matt — read this one. Sun et al. combine the multiscale outer loop from Schmitzer/Oberman-Ruan with their own IPRM inner solver and get real results: 2048x2048 OT instances, 7.1e7 active variables, KKT residual 7.5e-7, on a single machine. If that holds up, it's a genuine step forward for the field. The component-wise novelty is modest, but the integration and the support-identification theory are new. The nondegeneracy theorem (almost every marginal makes every optimal solution primal nondegenerate) is clean, and the Schur complement condition number bounds are useful.\n\nWhat the paper does well: the multiscale support inheritance works in practice, the active-set updates are simpler than MSSN's, and the warm-starting of IPRM from previous subproblem solutions is natural. The numerical comparisons against HOT, Sp-Sinkhorn, and network simplex are extensive and mostly fair; the 2048x2048 scaling experiment is the headline.\n\nNow the soft spots, in order of severity.\n\nFirst, the global optimality certificate. The termination criterion KKTres,out includes the dual infeasibility ||(c-A^T lambda)_-||/(1+||c||) over all m*n pairs. For a 2048x2048 grid that is 1.8e13 entries. The paper says cost entries are generated on demand and that the overhead is dominated by evaluating c-A^T lambda, but it never specifies how the maximum violation is found without enumerating everything. If the dual residual is only checked on the active set or a geometrically local candidate set, the reported 7.5e-7 does not certify optimality of the full LP. This matters directly for Proposition 2, which assumes complete knowledge of violated constraints to guarantee finite termination. The stress-test note has this right: it is a verification gap, not just a performance risk.\n\nSecond, the local quadratic convergence theorem assumes exact Newton solves and exact Schur complement solves. The implementation uses PCG. That gap is common in the IPM literature, but should be stated cleanly.\n\nThird, no code is released, and some tolerances and parameters (epsilon0, theta1, theta2, chi, IPRM's mu0, rho0, etc.) are not fully specified. That makes the headline experiments hard to reproduce.\n\nNone of these sink the paper. The central idea is sound and the experimental evidence is strong. But the trillion-variable claim should be phrased as 'solved to the accuracy certified by our screening procedure' until the global dual residual computation is specified and, ideally, implemented via an exact nearest-neighbor search (the squared Euclidean cost makes that feasible).\n\nRecommendation: send it to review. I'd like to see a revision that releases code, documents the residual evaluation, and reconciles the PCG practice with the exact-Newton theory. Worth citing as the current practical state of the art for large-scale discrete OT.","headline":"Strong algorithmic paper that likely makes discrete OT practical at new scale, but the trillion-variable claim rests on an uncertified global dual residual.","tokens_in":37828,"tokens_out":2934,"would_cite":true,"duration_ms":29275,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["15A03","65K05","90C06","90C08"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper claims that discrete optimal transport can be solved at grid resolutions whose full LP has trillions of variables, by solving only a chain of sparse subproblems.","keywords":["optimal transport","multiscale method","interior-point relaxation","support identification","large-scale linear programming","Schur complement","Wasserstein distance","warm start"],"falsifier":"Run MSIPRM on a cost and marginals where the fine-scale optimal plan must connect cells that are not children of any coarse-scale supporting edge, such as two narrow modes that only come into contact at fine resolution, and count active-set updates and memory growth as the grid is refined; if the active set grows to a large fraction of the full grid or the KKT residual stops improving, the support-inheritance assumption has failed.","tokens_in":36807,"feed_emoji":"🧮","tokens_out":8712,"duration_ms":82552,"temperature":0.7,"pith_summary":"Discrete optimal transport between two $2048\\times2048$ images has a full linear-programming formulation with about $18$ trillion transport variables, which no single machine can store. This paper claims that a multiscale outer loop can break that barrier: it solves coarse transport problems first, then uses the support of each coarse solution to predict a small set of active variables at the next finer scale, so each fine-scale problem is a sparse subproblem. An interior-point relaxation method serves as the inner solver, and because it accepts non-interior iterates it can warm start directly from the previous subproblem. The paper reports that a $2048\\times2048$ instance is solved with about $7.1\\times10^7$ active variables and a KKT residual of $7.5\\times10^{-7}$. If the sparse-support mechanism holds, high-resolution optimal transport ceases to be a memory-bound problem.","feed_headline":"Trillion-variable optimal transport solved on one machine","feed_subtitle":"A multiscale active-set scheme shrinks an 18-trillion-variable LP to 71 million variables at KKT residual 7.5e-7.","key_machinery":"Two mechanisms carry the argument. The first is the active-set inheritance rule $N^{[0]}_{\\ell} = \\bigcup_{(p,q)\\in\\operatorname{spt}(X^*_{\\ell+1})} \\operatorname{child}(u_{\\ell+1,p}) \\times \\operatorname{child}(v_{\\ell+1,q})$, which turns the coarse optimal support into the fine-scale variable set and is the reason the fine LP is small. The second is the primal-dual interior-point relaxation method, whose relaxed variables $z_j = \\tfrac12\\bigl(\\sqrt{(s_j/\\rho - x_j)^2 + 4\\mu/\\rho} - (s_j/\\rho - x_j)\\bigr)$ let iterates stay outside the interior while still tracing a perturbed central path; the Newton direction is computed through a Schur complement system whose condition number is shown to stay bounded as $\\mu\\to0$. The analysis connecting them is the support-identification theorem: for small enough $\\mu$, the sign of $z_j - \\sqrt{\\mu/\\rho}$ recovers the exact support of the optimal plan.","core_discovery":"The paper's central claim is that the optimal transport plan at fine resolution is sparse in a very specific, learnable way: its support is nearly contained in the union of child cells of the support of the coarser plan, so an active set built by that inheritance rule needs only a few enlargements before it contains an optimal solution. The paper argues that IPRM is the appropriate solver for the resulting subproblems, because its barrier augmented Lagrangian formulation replaces the interior-point strict positivity with a smooth relaxation of complementarity, and the perturbed primal-dual pair it traces makes warm starts and support thresholding rigorous. The theoretical part establishes that almost every marginal vector is primal nondegenerate for every connected active set, gives a threshold condition under which inexact IPRM iterates identify the exact support, proves bounded condition numbers for the Schur complement systems, and shows global convergence with error bound $O(\\sqrt{\\mu})$ plus local quadratic convergence. The numerical part reports that this combination solves instances whose full models contain trillions of variables, keeping the fine-scale active set to about $7.1\\times10^7$ variables on a $2048\\times2048$ grid.","pith_inferences":["Beyond the paper's own claims: if the 2048-by-2048 result is representative, the method makes Wasserstein-2 distances on high-resolution grids essentially routine on a single workstation, which would change how generative models and image barycenters are benchmarked.","The paper leaves implicit that the support-inheritance assumption is the true bottleneck: for costs whose fine-scale structure is invisible at the coarse scale, the active set would need many enlargements and the savings would erode; this can be tested by counting active-set updates on such costs.","The fine-scale overhead reported from evaluating $c - \\bar{A}^\\top\\lambda$ suggests that a parallel evaluation of that residual would directly multiply throughput, an engineering extension the paper does not pursue.","The support-identification threshold $\\sqrt{\\mu/\\rho}$ for $z_j$ is a cheap active-set predictor that could be reused as a warm-start heuristic for other LP solvers outside the multiscale setting."],"forward_implications":["The active-set update loop in Algorithm 2 is stated to terminate after finitely many enlargements, enumerating the supporting variables of a full optimal solution, with the objective value nonincreasing at each step.","Almost every marginal vector, in the measure-theoretic sense, makes every feasible connected sparse OT problem primal constraint nondegenerate, which implies uniqueness of the dual optimal solution and a rigorous basis for support identification.","IPRM enjoys a global convergence result with error bound $\\operatorname{dist}((x^{(k)},\\lambda^{(k)}),\\mathcal{S}) \\le C\\sqrt{\\mu^{(k)}}$ and quadratic local convergence under strict complementarity.","The Schur complement systems are provably well conditioned in the limit $\\mu\\to0$, so the linear algebra at the heart of the solver does not inherit classical IPM ill-conditioning.","Numerically, the method reports solved grid problems up to $2048\\times2048$ with full-model sizes around $1.8\\times10^{13}$ variables and active sets only a few tens of millions of variables."],"supporting_citations":[{"why":"supplies the primal-dual interior-point relaxation method and the barrier augmented Lagrangian function that IPRM is built on.","marker":"[20]"},{"why":"the multiscale semismooth Newton method this paper compares against and whose instability motivates the IPRM inner solver.","marker":"[21]"},{"why":"gives the multiscale support-transfer idea and the lemma used for warm-start initialization from coarser solutions.","marker":"[25]"},{"why":"the earlier sparse multiscale algorithm for dense optimal transport that this work extends.","marker":"[29]"},{"why":"the stabilized sparse Sinkhorn algorithm used as a multiscale baseline in the numerical comparisons.","marker":"[30]"},{"why":"the reduced-formulation solver used as a recent high-accuracy baseline.","marker":"[37]"},{"why":"supplies the smart-crossover lemma that supports identifying exact support from approximate optimal solutions.","marker":"[13]"},{"why":"provides the stability of optimal transport under perturbation that grounds the multiscale inheritance heuristic.","marker":"[34]"},{"why":"the Hoffman error bound used to prove IPRM's global error estimate.","marker":"[27]"},{"why":"active-set prediction under controlled perturbations, whose support-stability result this paper extends to nonpositive perturbations.","marker":"[4]"}],"fun_headline_variants":["Trillion-variable optimal transport solved via multiscale active sets","From trillions to millions: multiscale optimal transport solver","Support inheritance tames trillion-variable transport","Multiscale interior-point relaxation tackles trillion-variable optimal transport","Optimal transport with trillions of variables made tractable"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the fine-scale optimal transport plan is essentially supported on the child cells of the coarse-scale support, so the active set only needs a few enlargements; if that inheritance fails for a particular cost or marginals, the memory and time advantages of the method shrink.","fun_headline_variants_meta":{"raw":{"variants":["Trillion-variable optimal transport solved via multiscale active sets","From trillions to millions: multiscale optimal transport solver","Support inheritance tames trillion-variable transport","Multiscale interior-point relaxation tackles trillion-variable optimal transport","Optimal transport with trillions of variables made tractable"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000572,"raw_usage":{"total_tokens":2752,"prompt_tokens":1040,"completion_tokens":1712,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":656,"completion_tokens_details":{"reasoning_tokens":1628}},"tokens_in":656,"tokens_out":1712,"duration_ms":17186,"temperature":1.0,"reasoning_tokens":1628,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-16T00:20:02.403483+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run MSIPRM on a cost and marginals where the fine-scale optimal plan must connect cells that are not children of any coarse-scale supporting edge, such as two narrow modes that only come into contact at fine resolution, and count active-set updates and memory growth as the grid is refined; if the active set grows to a large fraction of the full grid or the KKT residual stops improving, the support-inheritance assumption has failed.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"supplies the primal-dual interior-point relaxation method and the barrier augmented Lagrangian function that IPRM is built on."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"the multiscale semismooth Newton method this paper compares against and whose instability motivates the IPRM inner solver."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"the earlier sparse multiscale algorithm for dense optimal transport that this work extends."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"the stabilized sparse Sinkhorn algorithm used as a multiscale baseline in the numerical comparisons."},{"cited_title":"IEEE Trans","cited_arxiv_id":null,"evidence_quote":"the reduced-formulation solver used as a recent high-accuracy baseline."},{"cited_title":"INFORMS J","cited_arxiv_id":null,"evidence_quote":"supplies the smart-crossover lemma that supports identifying exact support from approximate optimal solutions."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"provides the stability of optimal transport under perturbation that grounds the multiscale inheritance heuristic."},{"cited_title":"Linear Algebra Appl","cited_arxiv_id":null,"evidence_quote":"the Hoffman error bound used to prove IPRM's global error estimate."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"active-set prediction under controlled perturbations, whose support-stability result this paper extends to nonpositive perturbations."}],"review_version":1}