{"id":"7230a182-4542-4b84-9cd3-2d8fc2272e97","arxiv_id":"2608.09159","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":8,"one_line_summary":"PDHCG-CQP is a matrix-free GPU solver that provably converges linearly under strict complementarity for conic quadratic programs and scales to 440 million variables on 8 GPUs.","lead":"This paper presents PDHCG-CQP, a GPU-accelerated solver for large-scale conic quadratic programming, combining a restarted primal-dual hybrid gradient method with inexact inner solves and cone projections. It proves local linear convergence under strict complementarity and reports state-of-the-art robustness on standard benchmarks, scaling to 8 GPUs and hundreds of millions of variables.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Proposition 3.12's explicit dual-error constant is invalid: the bound on ∥M_ξ r_y∥ omits a √λ_ξ factor, so the stated κ_D and the αξ formula in (37) are not valid as written; the qualitative existence result survives only after correction.","rationale":"The reader's weakest assumption concerned the possibility of a face-transition boundary breaking the uniform calmness constants in Proposition 3.11. Under the paper's Corollary 3.13 hypothesis that every local KKT center satisfies x*+s* ∈ int K on an open neighborhood, the relevant faces are constant along the stratum via the cross-complementarity argument, and a compactness argument supplies a uniform positive margin; thus that specific concern does not land. The paper's central theorem is otherwise coherent: Theorem 2.5 is a valid perturbation argument, and the smoothed-gap machinery is standard. However, I identified a genuine, concrete error in the explicit constants of Proposition 3.12, which is part of the proof that Assumption 2.2 holds. This error does not invalidate the existence of a quadratic-growth constant, but it does mean the stated κ_D and the αξ in (37) are wrong as written; a corrected replacement (Lξ+κ_B)/√λξ is available. Additionally, the implemented solver replaces the Cesàro-averaged epoch with reflected-Halpern acceleration and adaptive restarts, which the convergence theorem does not cover; the paper is explicit about this scoping, but the title and abstract risk overstating the guarantee for the released code. In aggregate, the contributions are plausible and important, but the paper should be revised to correct the constant and to state clearly which implementation features are proven to preserve the rate. The reader's CONDITIONAL verdict remains appropriate; I do not see a reason to move it to ACCEPT or REJECT, so my recommendation is UNCHANGED.","tokens_in":35130,"tokens_out":33271,"duration_ms":317648,"concrete_test":"Verify the displayed inequality in Appendix B.7 with the scalar counterexample ξ=0.1, Q=0, Mξ=0.1I (so Lξ=0.1, λξ=0.1), κ_B=0, and any nonzero r. The claimed bound becomes 0.1∥r∥ ≤ 0.1√0.1 ∥r∥, which is false. Replacing the constant by (Lξ+κ_B)/√λξ fixes the bound; then recompute κ_D and the resulting αξ in (37) to confirm it remains positive, preserving the qualitative conclusion of Corollary 3.13.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The proof of Proposition 3.12 (Appendix B.7) bounds dist(s*(y), F_{x*}) by ∥Mξ r_y∥ + κ_B ∥r_y∥ and then claims this is at most (Lξ + κ_B √λξ) ∥r_y∥_{Mξ}, where Lξ = ∥Mξ∥ = ∥Q∥+ξ and λξ = λ_min(Mξ). This step is dimensionally inconsistent. Since ∥r_y∥ ≤ λξ^{-1/2} ∥r_y∥_{Mξ}, the correct upper bound is (Lξ/√λξ + κ_B/√λξ) ∥r_y∥_{Mξ} = (Lξ+κ_B)/√λξ · ∥r_y∥_{Mξ}. The paper's coefficient Lξ + κ_B√λξ can be smaller than the true coefficient when λξ < 1, making the displayed inequality false. For example, with ξ=0.1, Q=0, Mξ=0.1I, Lξ=0.1, λξ=0.1, and κ_B=0, the claimed inequality would assert ∥Mξ r∥ = 0.1∥r∥ ≤ 0.1·0.316∥r∥ = 0.0316∥r∥, which fails for every nonzero r. Consequently, the explicit κ_D in Proposition 3.12 and the αξ formula in (37)—which feed the contraction ratio in Theorem 2.5—are not established as stated. The existence of some finite κ_D is easily repaired by using the correct (Lξ+κ_B)/√λξ, so the central theorem is not void, but the paper's quantitative claim requires correction. This is a concrete flaw in the chain verifying the uniform local quadratic growth (Assumption 2.2), the essential bridge between strict complementarity and linear convergence.","agreement_with_reader":"disagree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes PDHCG-CQP, a matrix-free GPU solver for large-scale conic convex quadratic programs, and analyzes a restarted averaged primal-dual hybrid gradient method applied to the saddle formulation. The central theoretical claim is that, under a uniform local quadratic-growth condition on the smoothed duality gap, the restarted averaged method with exact or inexact primal proximal solves contracts the distance to the KKT set geometrically (Theorem 2.5). The paper then verifies the quadratic-growth condition for products of nonnegative, second-order, rotated second-order, and positive-semidefinite cones under blockwise strict complementarity, using a rotated-SOC epigraph lifting, a primal error bound, and a dual error bound derived from slack regularity and normal-cone calmness (Corollary 3.13). The practical solver adds reflected-Halpern acceleration, adaptive restarts, adaptive inner tolerances, Jacobi preconditioning, and a multi-GPU implementation, and the paper reports extensive experiments on QP, QCQP, SOCP, Lasso, and Fisher-equilibrium benchmarks, including instances with up to 4.4e8 stored primal coordinates.","tokens_in":35608,"tokens_out":16369,"duration_ms":151006,"significance":"If the theory is correct, the paper makes a valuable contribution: it extends the recent strict-complementarity-driven local linear convergence results from semidefinite programming to a broad class of conic quadratic programs, and it ships an open-source, matrix-free, single- and multi-GPU solver that is competitive with or better than existing first-order solvers on several large-scale benchmarks. The experiment section is unusually thorough, with public benchmarks, large synthetic instances, ablations, and distributed scaling results. The theoretical sections are structured around a clean sufficient condition (Assumption 2.2), and the appendices provide derivations for the main results. However, the paper currently contains a concrete quantitative error in the proof of the dual error bound, and the implemented algorithm differs from the analyzed one in ways that are not covered by the convergence theorems. These issues are repairable, and the qualitative local convergence claim appears defensible after correction.","major_comments":[{"comment":"The proof of Proposition 3.12 contains an invalid inequality. The text bounds dist(s*(y), F_x*) by ||M_xi r_y|| + kappa_B ||r_y|| and then claims this is at most (L_xi + kappa_B sqrt(lambda_xi)) ||r_y||_{M_xi}. Since ||r_y|| <= lambda_xi^{-1/2} ||r_y||_{M_xi}, a valid uniform bound from the preceding line is (L_xi + kappa_B)/sqrt(lambda_xi) ||r_y||_{M_xi}; the displayed coefficient L_xi + kappa_B sqrt(lambda_xi) can be smaller when lambda_xi < 1 and the inequality is false. For example, with Q=0, xi=0.1, M_xi=0.1I, L_xi=0.1, lambda_xi=0.1, and kappa_B=0, the claimed inequality would assert 0.1||r|| <= 0.0316||r||, which fails for every nonzero r. Consequently, the explicit constant kappa_D in Proposition 3.12 and the formula for alpha_xi in Eq. (37), which feed the contraction ratio in Theorem 2.5, are not established as written. The qualitative existence of a finite kappa_D is easily repaired with the corrected coefficient, so the central theorem is not void, but the displayed quantitative claims must be corrected.","section":"Appendix B.7 / Proposition 3.12 / Eq. (37)"},{"comment":"The convergence analysis applies to the restarted averaged scheme of Algorithm 1, but the implemented solver uses the reflected-Halpern recursion (46), the adaptive restart rule (61), the PID primal-weight update (58)-(59), and the adaptive inner tolerance (52)-(53). No theorem in Section 2 or Appendix A covers these recursions. The statement that the reflected-Halpern recursion is 'an implementation enhancement built on the same PDHG operator' does not transfer the local linear convergence guarantee. In particular, the relative proximal-error budget (18) is never shown to hold for the adaptive projected-BB inner loop; Proposition 2.6 analyzes fixed projected-gradient work only. The title and abstract promise local linear convergence of PDHCG-CQP, but the proved statement is for the simplified restarted averaged method. Please either extend the analysis to the implemented scheme or explicitly delimit the theoretical guarantee to the averaged method and present the enhanced solver as an empirical variant.","section":"Section 4.2.1 / Theorem 2.5"},{"comment":"The verification of uniform normal-cone calmness is presented as a sketch rather than a complete proof. For the PSD block, the argument introduces an unspecified absolute constant C, and the SOC boundary case is treated by saying it is 'the same computation' with lambda_min replaced by the distance from the cone vertex. Since Assumption 3.10 is one of the two regularity conditions needed for the dual error bound, and since Corollary 3.13 relies on the uniformity over the local KKT stratum, the proof should spell out the blockwise neighborhoods and moduli, including the boundary-ray case for SOC blocks and the product-cone assembly. This is a completeness issue rather than a demonstrated error, but it needs to be resolved before the quantitative verification of quadratic growth is complete.","section":"Appendix B.6 / Proposition 3.11"}],"minor_comments":[{"comment":"The statement of Assumption 3.10 uses a neighborhood W_B of (bar x, -bar s) but the paper does not explicitly note that the pairs (x*, -s*(y)) used in the dual error bound lie in that neighborhood after shrinking V_D and V_y; adding this sentence would make the uniformity argument easier to follow.","section":"Section 3.2.2 / Eq. (36)"},{"comment":"For reproducibility, please report the exact versions or commits of the baseline solvers (especially SCS, OSQP, PDCS, and the interior-point solvers) and the precise feasibility and optimality tolerances used by each solver, since first-order solver comparisons are sensitive to termination criteria.","section":"Section 5.1 / Experimental setup"},{"comment":"The projected Barzilai-Borwein method is described as the practical inner solver for general Q, but Proposition 2.6 analyzes projected gradient with a fixed stepsize and fixed number of steps. The paper should clarify that the BB variant and the adaptive stopping rule are empirical enhancements whose convergence properties are not covered by the theoretical results.","section":"Section 4.1 / Eq. (40)-(42)"},{"comment":"The SOCP benchmark comparison is reported at two tolerances only; reporting a third tolerance (for example 1e-8, as in the QCQP table) would give a more complete picture of the accuracy behavior on SOCP instances.","section":"Section 5.4 / Table 7"}],"recommendation":"major_revision","confidential_remarks":"The paper is competently written and the experimental contribution is substantial, but the proof error in Proposition 3.12 and the mismatch between the analyzed algorithm and the implemented one should be addressed before acceptance. The paper's self-citations to the PDHCG/PDQP/cuPDLP family are appropriate and do not appear to hide any conflict. The open-source release and the breadth of benchmarks are genuine strengths."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"PDHCG-CQP is worth knowing about: it is the first matrix-free GPU solver I know that keeps a native quadratic objective over general cones, and it scales to hundreds of millions of variables on eight GPUs. The theory is a genuine extension of the restarted-PDHG plus quadratic-growth framework to conic QP with inexact proximal solves, and the experiments on Maros-Meszaros, Mittelmann QP/QCQP/SOCP, and Fisher equilibria are extensive. Credit where due: the paper is open source, the ablation in Appendix E is honest, and the multi-GPU scaling numbers are impressive.\n\nThe soft spots, in order of seriousness. First, the stress-test is right: Proposition 3.12 as written contains a dimensionally inconsistent inequality. The proof bounds dist(s*(y), F_x*) by norm of M_xi r_y plus kappa_B norm r_y, and then claims a coefficient L_xi + kappa_B sqrt(lambda_xi) in the M_xi-norm. Since norm r_y <= lambda_xi^{-1/2} norm r_y_{M_xi}, the correct coefficient is (L_xi + kappa_B)/sqrt(lambda_xi). The stated kappa_D and the alpha_xi formula in (37) are therefore not established. The existence of a finite kappa_D is easily repaired, so the qualitative theorem survives, but the quantitative claim needs correction. A referee should require this fix.\n\nSecond, the uniformity in Proposition 3.11 is sketched, not proved, for SOC and PSD blocks, with 'absolute constants' left unspecified. The reader's worry about face-transition boundaries is legitimate: if a KKT center sits at a face transition, the calmness modulus can blow up and the dual error bound fails. The paper should either prove the uniform bound with explicit constants or restrict the statement to the constant-face regime.\n\nThird, the implemented solver uses reflected-Halpern, adaptive restart, PID stepsize control, and adaptive inner tolerances; the convergence theorem covers restarted averaged PDHG with a fixed epoch length and a relative-error budget. The paper is candid that these are enhancements, but it never quantifies the gap. That is normal for a solver paper, but the claims should be trimmed accordingly.\n\nFourth, benchmark reporting lacks statistical discipline: single seeds, no variances, and 'f' lumps time/memory/numerical failures. For the headline comparisons this matters less than for the timing claims.\n\nWho should read it: anyone working on first-order methods for conic optimization or GPU solvers. It deserves a serious referee; the theory is largely sound after the constant correction, and the engineering contribution is real.","headline":"A serious GPU conic-QP solver with a genuine local linear convergence theorem, but the explicit contraction constant in Proposition 3.12 is wrong as written and the implemented acceleration is outside the theory.","tokens_in":36119,"tokens_out":3018,"would_cite":true,"duration_ms":27253,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["90C20","90C25","90C06","65K05"],"pacs":[],"model":"deepseek-v4-flash","headline":"Restarted averaged PDHG converges locally linearly for conic quadratic programs under strict complementarity, and a GPU implementation scales to hundreds of millions of variables.","keywords":["conic quadratic programming","primal-dual hybrid gradient","local linear convergence","strict complementarity","GPU acceleration","matrix-free optimization","smoothed duality gap","large-scale quadratic programming"],"falsifier":"Run PDHCG-CQP on a family of small conic QPs whose optimal solution moves continuously from a strictly complementary interior-face point to a non-strictly-complementary boundary point, such as an SOCP where the optimal slack's smallest nonzero eigenvalue tends to zero. If the per-epoch distance-to-KKT residual stays geometric throughout the boundary family, the strict-complementarity regime is too narrow; if decay degrades to sublinear exactly at the transition, the central claim is corroborated. A direct algebraic check is to track the constants in the slack-regularity and normal-cone calmness assumptions along the family, since the smallest nonzero eigenvalue or boundary-ray distance should approach zero precisely at the transition.","tokens_in":34944,"feed_emoji":"⚡","tokens_out":7966,"duration_ms":71750,"temperature":0.7,"pith_summary":"This paper establishes that a restarted averaged primal-dual hybrid gradient method, PDHCG-CQP, converges locally linearly for conic convex quadratic programming whenever the smoothed primal-dual gap grows quadratically near the solution set. It then proves that this quadratic growth holds under strict complementarity for products of nonnegative, second-order, rotated second-order, and positive-semidefinite cones. The convergence proof covers inexact proximal solves, so each iteration needs only matrix-vector products and cone projections, not factorizations. The result is a matrix-free GPU solver with a proven geometric contraction rate, demonstrated on instances with up to 4.4×$10^{8}$ stored primal coordinates and on multiple GPUs. This matters because previous first-order conic solvers either reformulate the quadratic objective at extra cost or lack a local linear convergence guarantee under strict complementarity.","feed_headline":"Conic QP on GPUs now has a proven local linear rate","feed_subtitle":"No factorizations: restarted averaged PDHG shrinks error geometrically and scales to 440 million variables.","key_machinery":"The load-bearing object is the smoothed duality gap $G_\\xi(z; z_\\star)$ of Definition 2.1, which penalizes the distance of the comparison point to a fixed center. One epoch's Cesàro averaging converts PDHG energy inequalities into an $O(1/T)$ bound on this gap, and uniform quadratic growth of the gap, Assumption 2.2, converts that bound into a geometric contraction of the distance to the KKT set. The paper verifies Assumption 2.2 by decomposing the smoothed gap into primal and dual error bounds, proving the primal bound through a rotated second-order-cone epigraph lifting of $Q = B^*B$ combined with complementarity error bounds under primal strict complementarity, and proving the dual bound through slack-set linear regularity and normal-cone calmness under dual strict complementarity. Strict complementarity makes the relevant complementary faces constant along the local KKT stratum, which supplies the uniform constants the convergence theorem needs.","core_discovery":"The central claim is Theorem 2.5: under uniform local quadratic growth of the smoothed gap, each epoch of restarted averaged PDHG contracts the distance to the KKT set by a fixed factor, $d_n \\le q_\\eta^n d_0$, even when the primal proximal subproblems are solved only inexactly within a relative error budget. Corollary 3.13 verifies the required quadratic growth: if every local KKT center in a neighborhood is strictly complementary and the cone is a finite product of nonnegative, Lorentz (including rotated Lorentz), and positive-semidefinite cones, then the smoothed gap is bounded below by a constant times the squared distance to the KKT set. Together these statements extend the principle of global convergence followed by strict-complementarity-driven local linear convergence, previously established for semidefinite programming, to conic quadratic programs with a native quadratic objective.","pith_inferences":["A testable prediction is that on instances where the limiting KKT point sits at a face-transition boundary and strict complementarity fails, the residual decay should slow from geometric to sublinear, revealing exactly where the local regularity assumptions break.","The smoothed-gap quadratic-growth certificate is solver-independent, so the same Assumption 2.2 could be checked for other first-order conic methods to give them the same local rate whenever their epoch averages satisfy an analogous gap bound.","The paper's explicit relation between inner work and the condition number of the proximal subproblem suggests a practical diagnostic: measuring proximal error decay as a function of projected-gradient steps should match the predicted geometric factor, letting users choose inner iterations without tuning.","Because the convergence guarantee is local, the practical outcome depends on whether global PDHG dynamics drive iterates into a strictly complementary stratum; warm-starting from a cheap first-order run could be the principled way to reach that regime."],"forward_implications":["The method is a matrix-free, GPU-resident conic QP solver with a proven local linear rate, so users get the scalability of PDHG without giving up a convergence guarantee near strictly complementary solutions.","Quadratic objectives are kept natively, avoiding the rotated-second-order-cone epigraph reformulation that linear-objective conic solvers need; the paper's comparisons show that reformulation roughly triples the outer iteration count and adds per-iteration lifted projections.","Inexact proximal solves preserve the geometric rate: a fixed number of warm-started projected-gradient steps per subproblem, logarithmic in the restart length, is enough to satisfy the error budget.","The theory covers products of nonnegative, second-order, rotated second-order, and positive-semidefinite cones, while the implementation additionally supports exponential and power cones with the same algorithm.","On the tested benchmarks the solver solves more standard QP instances than the other tested first-order solvers and is the only solver reported as reaching optimality on the largest market-equilibrium instances, which have over 4.4×10^8 stored primal coordinates."],"supporting_citations":[{"why":"Supplies the blockwise complementarity error bounds used to verify the lifted residual error bound in Proposition 3.6.","marker":"[15]"},{"why":"Establishes the smoothed-gap quadratic error bound framework that turns ergodic PDHG estimates into restart contraction.","marker":"[18]"},{"why":"Shows local linear convergence of PDHG for semidefinite programming under strict complementarity, the principle this paper extends to conic QP.","marker":"[26]"},{"why":"Provides the analogous ADMM strict-complementarity local linear result for semidefinite programming that motivates the conic extension.","marker":"[27]"},{"why":"The closest PDHG-based conic solver, used as baseline and as the source of rescaled-cone projection formulas for exponential and power cones.","marker":"[30]"},{"why":"Predecessor PDHCG method for large-scale QP whose proximal-subproblem design is carried over to the conic setting.","marker":"[25]"},{"why":"Gives C^2-cone reducibility conditions used to verify normal-cone calmness for second-order and semidefinite blocks.","marker":"[33]"},{"why":"Provides polyhedral multifunction calmness used for nonnegative cone blocks in the uniform normal-cone calmness verification.","marker":"[50]"}],"fun_headline_variants":["GPU conic QP solver: linear rate proven","PDHCG: GPU solver with local linear convergence for conic QP","Scaling conic QP to 440M variables with proven linear rate","First-order GPU solver for conic QP: linear convergence under strict complementarity","Restarted averaged PDHG on GPUs: linear rate for large-scale conic QP"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The linear rate relies on the assumption that every KKT point near the solution is strictly complementary and that the resulting regularity constants stay bounded uniformly; if the iterates approach a solution on a face-transition boundary where strict complementarity fails, the uniformity breaks and the contraction proof no longer applies.","fun_headline_variants_meta":{"raw":{"variants":["GPU conic QP solver: linear rate proven","PDHCG: GPU solver with local linear convergence for conic QP","Scaling conic QP to 440M variables with proven linear rate","First-order GPU solver for conic QP: linear convergence under strict complementarity","Restarted averaged PDHG on GPUs: linear rate for large-scale conic QP"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000784,"raw_usage":{"total_tokens":3495,"prompt_tokens":1015,"completion_tokens":2480,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":631,"completion_tokens_details":{"reasoning_tokens":2380}},"tokens_in":631,"tokens_out":2480,"duration_ms":15078,"temperature":1.0,"reasoning_tokens":2380,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-11T22:20:28.182598+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run PDHCG-CQP on a family of small conic QPs whose optimal solution moves continuously from a strictly complementary interior-face point to a non-strictly-complementary boundary point, such as an SOCP where the optimal slack's smallest nonzero eigenvalue tends to zero. If the per-epoch distance-to-KKT residual stays geometric throughout the boundary family, the strict-complementarity regime is too narrow; if decay degrades to sublinear exactly at the transition, the central claim is corroborated. A direct algebraic check is to track the constants in the slack-regularity and normal-cone calmness assumptions along the family, since the smallest nonzero eigenvalue or boundary-ray distance should approach zero precisely at the transition.","supporting_citations":[{"cited_title":"A strict complementarity approach to error bound and sensitivity of solution of conic programs.Optimization Letters, 17(7):1551–1574, 2023","cited_arxiv_id":null,"evidence_quote":"Supplies the blockwise complementarity error bounds used to verify the lifted residual error bound in Proposition 3.6."},{"cited_title":"Local Linear Convergence of the Primal-Dual Hybrid Gradient Method for Semidefinite Programming","cited_arxiv_id":"2607.08035","evidence_quote":"Shows local linear convergence of PDHG for semidefinite programming under strict complementarity, the principle this paper extends to conic QP."},{"cited_title":"A restarted primal-dual hybrid conjugate gradient method for large-scale quadratic programming.INFORMS Journal on Computing, 2025","cited_arxiv_id":null,"evidence_quote":"Predecessor PDHCG method for large-scale QP whose proximal-subproblem design is carried over to the conic setting."},{"cited_title":"Computation of graphical derivatives of normal cone maps to a class of conic constraint sets.Set-Valued and Variational Analysis, 27(3):783–806, 2019","cited_arxiv_id":null,"evidence_quote":"Gives C^2-cone reducibility conditions used to verify normal-cone calmness for second-order and semidefinite blocks."},{"cited_title":"Robinson","cited_arxiv_id":null,"evidence_quote":"Provides polyhedral multifunction calmness used for nonnegative cone blocks in the uniform normal-cone calmness verification."}],"review_version":1}