GPU-Accelerated Conic Quadratic Programming with Local Linear Convergence under Strict Complementarity

GPU-ACCELERATED C ONIC Q UADRATIC P ROGRAMMING WITH L OCAL L INEAR C ONVERGENCE UNDER S TRICT C OMPLEMENTARITY

Hongpei Li1 , Yicheng Huang2 , Huikang Liu3 , Dongdong Ge3 , and Yinyu Ye3,4 1

arXiv:2608.09159v1 [math.OC] 10 Aug 2026

2

Northwestern University Shanghai University of Finance and Economics 3 Shanghai Jiao Tong University 4 Stanford University

A BSTRACT We present PDHCG-CQP, a GPU-accelerated first-order solver for large-scale conic convex quadratic programming. PDHCG-CQP supports affine constraints and Cartesian products of nonnegative, second-order, rotated second-order, exponential, and three-dimensional power cones. At its core is a restarted averaged primal-dual hybrid gradient (PDHG) method, whose primal update is computed inexactly by solving a conic quadratic proximal subproblem with projected gradient iterations. We establish local linear convergence of the restarted averaged scheme with both exact and inexact primal proximal evaluations under a uniform local quadratic-growth condition on the smoothed primal-dual gap. We further show that this condition holds under strict complementarity by exploiting a rotated second-order-cone lifting together with local primal and dual regularity conditions. Our C/CUDA implementation combines matrix-free linear algebra, batched cone projections, adaptive inner solves, reflected-Halpern acceleration, and fully device-resident KKT residual computations. It also supports multi-GPU execution through a two-dimensional partitioning of the problem data. Extensive experiments on standard and large-scale quadratic programming (QP), convex quadratically constrained quadratic programming (QCQP), second-order cone programming (SOCP), and quasilinear Fisher equilibrium benchmarks demonstrate that PDHCG-CQP achieves state-of-the-art robustness among first-order solvers while scaling efficiently to 8 GPUs and instances with up to 4.4 × 108 stored primal coordinates. PDHCG-CQP is open source and available at https://github.com/Lhongpei/PDHCG.

1

Introduction

Conic convex quadratic programming (CQP) is a fundamental optimization problem class with applications in machine learning [53], control [42], signal processing [34], finance [40], and market design [16]. A quadratic objective represents curvature, regularization, or risk, while conic constraints model geometric structures such as norms, perspectives, exponential relations, and power laws. Formally, we consider the conic convex quadratic program 1 ⟨x, Qx⟩ + ⟨c, x⟩ 2 s.t. Ax = b, x ∈ K,

min f (x) := x∈E

(1)

P REPRINT. U NDER REVIEW.

where E and Y are finite-dimensional Euclidean spaces, A : E → Y is linear, Q : E → E is self-adjoint and positive semidefinite, c ∈ E, b ∈ Y , and K ⊆ E is a nonempty closed convex cone. The associated equality-conic saddle function is 1 L(x, y) := ⟨x, Qx⟩ + ⟨c, x⟩ + ⟨y, b − Ax⟩, x ∈ K, y ∈ Y, (2) 2 and a saddle point of (2) recovers a primal-dual solution of (1). Classical conic and quadratic-programming solvers are dominated by interior-point and active-set methods. Modern interior-point solvers such as Clarabel [23], CuClarabel [10], and QOCO-GPU [8] provide reliable high-accuracy solutions, with the latter two accelerating sparse factorizations on GPUs. Their scalability, however, is limited by repeated Newton–KKT solves, whose sparse factorizations can become expensive in both time and memory because of fill-in and are difficult to distribute efficiently across multiple accelerators. First-order splitting methods alleviate this bottleneck: SCS [47, 48], ABIP+ [14, 29], and OSQP [52] reuse a single factorization or employ indirect linear solves, but the former can still be prohibitive at scale while the latter requires multiple matrix–vector products per iteration. Fully matrix-free restarted primal-dual hybrid gradient (PDHG) methods eliminate linear-system solves altogether and have demonstrated strong scalability for LP in PDLP [2] and cuPDLP [38], and for box-constrained convex QP in PDQP [36], HPR-QP [9], and PDHCG [25]. For conic feasible sets, this matrix-free toolbox is still incomplete. PDCS [30], the closest PDHG-based conic solver, treats a linear objective and handles quadratic objectives through an epigraph reformulation, which requires a factorization Q = B ∗ B, roughly triples the outer iteration count, and adds the projection of a rotated second-order-cone (RSOC) block of dimension at least rank(Q) + 2 to every iteration; Section 2 quantifies these effects. Moreover, the current PDCS implementation does not yet convert the matrix-free design into a practical advantage: in the conic benchmarks of Section 5, it is outperformed by a wide margin — even by factorization-based solvers on the large instances for which first-order methods are intended. No available solver combines a native quadratic objective, general vector-level cones, and matrix-free GPU execution. We therefore seek a method with the following properties. • Matrix-free computation. Each outer iteration uses only applications of A, A∗ , and Q, together with productcone projections; no factorization is required. • Native conic quadratic modeling. The quadratic objective is kept explicitly, and the principal vector-level cones are supported without an objective-epigraph reformulation. • Linear convergence with implementable proximal steps. The convergence theory covers both the exact primal proximal map and practical inexact solutions of the conic quadratic subproblems. • Scalable implementation. The computational kernels map naturally to a single GPU and admit a distributed extension when the problem exceeds the memory of one device. This paper develops PDHCG-CQP, a matrix-free solver for (1) with these four properties. The algorithm is restarted averaged PDHG applied directly to the saddle problem (2): each iteration performs a primal-first PDHG update whose primal step is a strongly convex conic quadratic proximal subproblem, and each epoch restarts from the Cesàro average of its iterates. The proximal subproblem generally has no closed form; it is solved by direct weighted cone projections when a separable oracle is available and by matrix-free projected-gradient iterations otherwise, and the convergence analysis covers both exact and inexact proximal evaluations, with the latter satisfying a stated error budget. The analysis rests on a single local regularity property: uniform quadratic growth of the smoothed duality gap near the solution set. A single PDHG iterate provides no usable gap estimate, but the one-step energy inequalities telescope into an O(1/T ) ergodic bound at the epoch average, and quadratic growth converts this bound into a geometric contraction of the distance to the KKT set; a perturbation argument extends the contraction to inexact proximal solves. Quadratic growth, in turn, is not automatic for nonpolyhedral cones. We characterize it exactly through a primal and a dual error bound and verify both under strict complementarity, so that the “global convergence plus strict-complementarity-driven

2

P REPRINT. U NDER REVIEW.

local linear convergence” principle, recently established for semidefinite programming (SDP) [26, 27], extends to the general conic program with a quadratic objective. The contributions of this paper can be summarized as follows. • Quadratic growth for conic convex QP under strict complementarity. We characterize uniform local quadratic growth of the smoothed primal-dual gap exactly through a primal and a dual error bound, and we verify both bounds — through a rotated-SOC lifting, the error bound of [15], slack regularity, and normal-cone calmness — when K is a finite product of nonnegative, second-order, and positive-semidefinite cones whose local KKT centers are strictly complementary. • Local linear convergence of restarted averaged PDHG with an inexact inner solver. Under this quadratic growth, restarting from epoch averages contracts the distance to the KKT set geometrically. A solverindependent relative proximal-error budget preserves the rate, and a fixed number of warm-started projectedgradient steps per subproblem satisfies the budget. • A state-of-the-art GPU solver. Our matrix-free CUDA implementation combines structured quadratic products, batched cone projections, adaptive inner tolerances, and a multi-GPU extension, delivering state-ofthe-art performance on the tested large-scale conic QP benchmarks. The experiments in Section 5 support these claims. On the Maros-Mészáros benchmark [41], PDHCG-CQP solves 126 instances at tolerance 10−6 , the most among the tested solvers, with the best average runtime; on the Mittelmann QP benchmark [43] it solves 17 instances at both 10−6 and 10−8 , more than any other tested solver. On the public convex quadratically constrained quadratic-programming (QCQP) benchmark it is the only first-order solver that solves every instance through 10−6 , and on the Mittelmann second-order cone-programming (SOCP) benchmark [44] it attains the best solved counts among the tested first-order methods. On large sparse Lasso QPs it is the fastest solver on seven of the nine instances. On large-scale quasilinear Fisher equilibrium instances, the distributed implementation achieves up to a 7.97× speedup on eight GPUs and is the only solver reported as reaching optimality at n = 107 buyers; the large Lasso instances similarly benefit from the multi-GPU implementation of PDHCG-CQP. 1.1

Related literature

Conic and quadratic-programming solvers. Interior-point and active-set methods remain the standard approaches for obtaining high-accuracy solutions of conic and quadratic programs. Mature commercial solvers, including MOSEK [45], Gurobi [24], and COPT [21], provide highly optimized interior-point implementations for conic and convex quadratic programs and serve as high-accuracy baselines in our experiments. Clarabel [23] provides an interior-point method for quadratic objectives over general convex cones, while CuClarabel [10], QOCO-GPU [8], and NVIDIA’s cuOpt [46] investigate GPU acceleration of sparse interior-point linear algebra. These approaches are complementary to the matrix-free regime studied here. First-order and splitting-based solvers. SCS [47, 48] and ABIP+ [14, 29] solve homogeneous self-dual formulations using ADMM, while OSQP [52] specializes ADMM to convex QPs. Their direct implementations rely on a reusable factorization of a fixed linear system; indirect variants replace the factorization by iterative linear solves but typically require several matrix-vector products per outer iteration. PDCS [30] instead provides a matrix-free, GPU-oriented PDHG method for linear-objective conic programs and handles convex QPs through conic reformulation. Our method retains the quadratic objective in the primal PDHG proximal step and therefore avoids the corresponding objectiveepigraph lifting. Restarted primal-dual methods. Restarted PDHG underlies the PDLP family [2] for large-scale LP, and GPU implementations such as cuPDLP [35] and cuPDLP-C[38] demonstrate the scalability of this approach. PDQP [36] extends related first-order ideas to convex QP, while HPR-QP [9] and the earlier PDHCG method [25] provide further restarted primal-dual QP algorithms. Restarted Halpern PDHG and its reflected variant [37] have also been analyzed for 3

P REPRINT. U NDER REVIEW.

LP. Concurrent work by [32] studies a broader reflected-Halpern framework under fixed-point sharpness. The present paper instead analyzes restarted averaged PDHG for the conic saddle problem through the smoothed duality gap, with explicit perturbation guarantees for inexact quadratic proximal solves, and establishes the required quadratic growth from conic primal and dual geometry. 1.2

Notation

The adjoint of A is denoted by A∗ , the dual cone of K by K ∗ , and the normal cone of K at x by NK (x). For a primal-dual pair (x, y), define the dual slack s(x, y) := Qx + c − A∗ y. A point (x⋆ , y ⋆ ) is a saddle point if and only if, with s⋆ = s(x⋆ , y ⋆ ), Ax⋆ = b, x⋆ ∈ K, s⋆ ∈ K ∗ , ⟨x⋆ , s⋆ ⟩ = 0. (3) We denote the saddle-point set and its primal and dual projections by Z ⋆ := {(x, y) : (x, y) satisfies (3)},

X ⋆ := projE Z ⋆ ,

Y ⋆ := projY Z ⋆ .

Throughout the convergence analysis, we assume Z ⋆ ̸= ∅ and fix a reference KKT point z̄ ∈ Z ⋆ . All unqualified norms are Euclidean, and product spaces use the corresponding product norm. For a nonempty set S, let dist(u, S) := inf ∥u − v∥,

ΠS (u) := argmin ∥u − v∥

v∈S

v∈S

whenever the projection is single-valued. For r > 0 and z ∈ E × Y , define Br (z) := {w ∈ E × Y : ∥w − z∥ < r},

2

B r (z) := {w ∈ E × Y : ∥w − z∥ ≤ r}.

PDHCG for Conic Quadratic Programming

In this section, we present PDHCG-CQP, which directly applies the restarted averaged PDHG method to the conic saddle-point problem (2). Algorithm 1 PDHCG-CQP: restarted averaged PDHG for the original conic QP 1: Input: z 0,0 = (x0,0 , y 0,0 ), stepsizes τ, σ > 0, an epoch length T . 2: for n = 0, 1, 2, . . . do 3: Initialize the running average z̄ n,0 = z n,0 . 4: for k = 0, 1, . . . , T − 1 do 5: Compute exactly or approximately

xn,k+1 ≈ arg min x∈K

6: 7:




2 1 1 x − xn,k + τ A∗ y n,k . ⟨x, Qx⟩ + ⟨c, x⟩ + 2 2τ

Set xn,k+1 = 2xn,k+1 − xn,k and y n,k+1 = y n,k + σ(b − Axn,k+1 ). r r Update the running average: z̄ n,k+1 =

1 k z̄ n,k + z n,k+1 . k+1 k+1

8: end for 9: Restart from the epoch average: z n+1,0 = z̄ n,T . 10: end for

Algorithm 1 has two nested loops. Within an epoch, the inner loop performs primal-first PDHG updates and maintains the running Cesàro average of its iterates; the outer loop restarts each epoch from the previous epoch average. Since every exact or feasible inexact primal iterate lies in K and K is convex, the averaged primal component remains conically feasible. The displayed minimization is either exact, when it admits a closed form, or inexact, computed by an arbitrary inner solver. The theory below covers both cases.

4

P REPRINT. U NDER REVIEW.

Benefit of PDHCG. The predecessor of the present method, PDHCG [25], was proposed for large-scale convex QP as a restarted primal-dual method that solves the strongly convex primal proximal subproblem accurately by conjugate-gradient-type inner iterations, rather than taking a single forward gradient step. Treating the quadratic term through its proximal subproblem moves the curvature of Q from the outer loop into cheap inner iterations and thereby reduces the number of outer PDHG iterations substantially. The same design principle underlies PDHCG-CQP. Retaining the quadratic objective natively is also preferable to eliminating it. If Q = B ∗ B, problem (1) admits the rotated second-order-cone (RSOC) epigraph reformulation min t + ⟨c, x⟩

s.t.

x,t

Ax = b,

x ∈ K,

(t, 1, Bx) ∈ Qr ,

where Qr := {(u, v, w) : u, v ≥ 0, 2uv ≥ ∥w∥2 }, after which any conic solver with a linear objective applies. Three considerations argue against this route. First, unless a factor B is supplied, the reformulation requires a potentially expensive and fill-inducing factorization of Q. Second, the reformulation degrades the iteration path itself: the outer iteration count grows — roughly threefold in Table 1 — and every outer iteration must additionally project onto an RSOC block of dimension at least rank(Q) + 2, which also adds storage and communication costs. These per-iteration lifted projections are typically more expensive than the product-set projections performed inside the projected-gradient inner solves of the native formulation. Third, the native formulation exposes the curvature of Q directly through the proximal Hessian Q + τ −1 I, which could improve practical convergence when this curvature is informative, although the benefit is problem-dependent. The test problems in Table 1 are drawn from the QPLIB-QCQP benchmark collection maintained by Hans Mittelmann [43]. On these instances, the native formulation is approximately three times faster and uses about one third as many iterations as the RSOC reformulation. Table 1: Native quadratic objective versus the RSOC reformulation on the Mittelmann QCQP benchmark instances. Target accuracy 10−4 10−6 10−8

Native quadratic

RSOC reformulation

Time

Iterations

Time

Iterations

4.67 37.83 82.24

25,811 310,074 687,301

13.28 118.04 252.46

76,836 868,239 1,911,390

The rest of this section proves eventual local linear convergence of PDHCG-CQP. Its local convergence analysis first uses the exact PDHG operator to establish the basic epoch contraction and then passes to an implementable inexact PDHG operator, in which the primal proximal subproblem may be approximated by any inner method satisfying the error conditions stated below. The conic lifting introduced later in Section 3 is used only to verify a quadratic-growth property. 2.1

Preliminaries

For saddle-point problems, the primal-dual gap is a standard merit function: it compares the Lagrangian at the current primal and dual variables against a primal-dual comparison point. Ergodic convergence estimates for PDHG are naturally expressed through this two-point gap [5]. For z = (x, y) ∈ K × Y and zb = (b x, yb) ∈ K × Y , define Q(z, zb) := L(x, yb) − L(b x, y).

(4)

Maximizing this comparison function over zb gives the ordinary duality gap. Although it is a natural primal-dual progress measure, this supremum can be infinite when the primal or dual domain is unbounded. Following the smoothed-gap framework of [18], we instead penalize the distance from the comparison point to a prescribed center. The resulting quantity remains finite and can satisfy a quadratic error bound, thereby linking the ergodic PDHG estimate to restart

5

P REPRINT. U NDER REVIEW.

contraction. This framework has also been used in restarted first-order methods for convex QP, including rAPDHG and the earlier PDHCG method [25, 36]. Definition 2.1 (smoothed duality gap). For ξ > 0, z = (x, y) ∈ K × Y , and ż = (ẋ, ẏ) ∈ E × Y , define   ξ ξ 2 2 Gξ (z; ż) := sup Q(z, zb) − ∥b x − ẋ∥ − ∥b y − ẏ∥ . 2 2 x b∈K, y b∈Y

(5)

The quadratic penalty makes the supremum finite and is also the bridge between an ergodic gap estimate and distance to the KKT set. Assumption 2.2 (uniform local quadratic growth). There are a KKT point z̄ ∈ Z ⋆ , constants R > 0, ξ > 0, and αξ > 0 such that Gξ (z; z ⋆ ) ≥ αξ dist2 (z, Z ⋆ ) (6) for every z ∈ (K × Y ) ∩ BR (z̄),

z ⋆ ∈ Z ⋆ ∩ BR (z̄).

The center z ⋆ in (6) is allowed to vary over the local solution stratum, while the constant αξ and the effective neighborhood remain fixed. Section 3 verifies this assumption from primal and dual geometric error bounds. Thus Assumption 2.2 is the only problem-dependent local regularity property used in the convergence proof below. The proof has two stages. Stage I assumes that every primal proximal subproblem is solved exactly. A one-step energy inequality gives an O(1/T ) gap bound for one Cesàro-averaged epoch, and local quadratic growth turns this bound into a strict contraction. Stage II treats the finite accuracy of a general inexact proximal oracle as a perturbation of the exact epoch and gives method-independent conditions under which the same local linear rate is retained. Projected gradient is then presented only as one concrete oracle satisfying those conditions. 2.2

Exact averaged epochs

We first analyze the exact scheme, in which every primal proximal subproblem is solved exactly; it is the reference trajectory against which the inexact method of Stage II is compared. Throughout this stage the stepsizes and the epoch length are fixed. Set 1 F (x) := ⟨x, Qx⟩ + ⟨c, x⟩ + δK (x), 2 and fix τ, σ > 0 such that τ σ∥A∥2 < 1. (7) One primal-first PDHG step from z k = (xk , y k ) is xk+1 = proxτ F (xk + τ A∗ y k ),

(8a)

xk+1 = 2xk+1 − xk , r k+1 k

(8b)

y

=y

+ σ(b − Axk+1 ), r

(8c)

denoted T (z k ) = z k+1 . Because the quadratic term is treated proximally, the stability condition (7) involves A but not ∥Q∥. During restart epoch n, perform T steps z n,k+1 = T (z n,k ) and restart from the Cesàro average T −1

z̄ n,T :=

1 X n,k+1 z , T

z n+1,0 := z̄ n,T .

(9)

k=0

Averaging is what produces the required gap estimate: the one-step energy inequalities telescope into an O(1/T ) ergodic bound at z̄ n,T , whose primal component moreover remains in K by convexity. Joint primal-dual progress is measured in the metric of the symmetric preconditioner " # τ −1 I A∗ P := , A σ −1 I 6

(10)

P REPRINT. U NDER REVIEW.

which is positive definite under (7). In this metric the exact update is a preconditioned resolvent of the KKT operator and hence Fejér monotone with respect to Z ⋆ (Lemma A.1). With p λ− := λmin (P ), λ+ := λmax (P ), κP := λ+ /λ− , (11) write, for the exact restart sequence, s dn := dist(z

n,0

⋆

, Z ),

q0 :=

λ+ , αξ T

R0 :=

1 + κP d0 . 1 − e−1

Theorem 2.3 (local linear convergence of exact restarted PDHG). Suppose d0 > 0 and Assumption 2.2 holds on BR (z̄), with B R0 (z 0,0 ) ⊆ BR (z̄). (12) If the restart length satisfies  T ≥ max

2λ+ e2 λ+ , ξ αξ

 ,

(13)

n ≥ 0.

(14)

then the distance to the KKT set decreases geometrically: dist(z n,0 , Z ⋆ ) ≤ e−n d0 ,

Of the two lower bounds in (13), the first permits the conversion of the ordinary gap into the smoothed gap, and the second gives q0 ≤ e−1 , so that each exact epoch contracts by at least e−1 . The complete proof is deferred to Appendix A.1. The local initialization condition in Theorem 2.3 should be understood as an eventual condition rather than as a requirement that a user provide a warm start. For the exact reference scheme, standard PDHG convergence theory ensures convergence from an arbitrary initialization to some saddle point z ∞ ∈ Z ⋆ under the usual stepsize condition [5, 6]; finite Cesàro averaging and restart preserve the saddle-point fixed set. Once the iterates enter a sufficiently small neighborhood of z ∞ , that iterate may be reindexed as z 0,0 in Theorem 2.3. Hence the substantive local requirement is that the saddle point selected by the global dynamics admit uniform local quadratic growth. Under the sufficient conic regularity conditions developed in Section 3, strict complementarity at the limiting KKT point is the central mechanism that yields this property. Accordingly, within the strict-complementarity regime covered by Section 3, the local initialization hypothesis amounts, after a finite burn-in, to assuming that the limiting KKT point is strictly complementary. Strict complementarity. This interpretation parallels recent SDP results: PDHG is eventually R-linearly convergent when its limiting KKT point satisfies strict complementarity [26], and ADMM is locally linearly convergent when its limiting primal-dual solution is strictly complementary [27]. Thus, subject to the conic regularity conditions in Section 3, our result extends the same “global convergence followed by strict-complementarity-driven local linear convergence” principle from SDP to conic convex QP with a quadratic objective. 2.3

Inexact proximal solves

We now separate the convergence argument from the choice of inner solver. For an input v ∈ E, define p(v) := proxτ F (v) as the unique minimizer over K of the strongly convex function Hv (u) :=

1 1 ⟨u, Qu⟩ + ⟨c, u⟩ + ∥u − v∥2 . 2 2τ

(15)

Let zek = (e xk , yek ) be an inexact orbit and set vek := x ek + τ A∗ yek , 7

pk+1 v k ). sh := p(e

(16)

P REPRINT. U NDER REVIEW.

At step k, an arbitrary inner solver returns a point x ek+1 ∈ K. Its proximal error is δk := ∥e xk+1 − pk+1 sh ∥.

(17)

The extrapolated primal point and dual update are then x ek+1 = 2e xk+1 − x ek , r

yek+1 = yek + σ(b − Ae xk+1 ). r

For a relative-error parameter η ≥ 0, set εn :=

T −1 X

δn,k ,

qη := q0 + Cep η,

βη := 1 + κP + Cep η,

k=0

Rη :=

βη d0 , 1 − qη

where the constant Cep is supplied in Lemma 2.4 in Appendix A.2. Theorem 2.5 (local linear convergence with inexact proximal steps). Suppose Assumption 2.2 holds on BR (z̄), let T ≥ 2λ+ /ξ, and suppose εn ≤ ηdn , qη < 1. (18) If the initial point z 0,0 satisfies B Rη (z 0,0 ) ⊆ BR (z̄),

(19)

then the distance to the KKT set decreases geometrically: dn ≤ qηn d0 ,

n ≥ 0.

(20)

The complete proof, together with the perturbation lemma on which it relies, is deferred to Appendix A.2. The relative-error condition (18) is not directly implementable: εn aggregates distances to unknown exact proximal points, and dn is the distance to the unknown KKT set. The next subsection replaces both by computable certificates and a movement-based inner tolerance. 2.3.1

Accuracy certificates

The error δk involves the unknown exact point pk+1 sh , but it is controlled by computable certificates. Set Hτ := Q + τ −1 I,

µτ := λmin (Hτ ) ≥ τ −1 ,

Lτ := λmax (Hτ ) = ∥Q∥ + τ −1 .

If the inner algorithm returns x e ∈ K for the proximal input v together with a stationarity residual r ∈ Qe x + c + τ −1 (e x − v) + NK (e x),

(21)

then, because the operator u 7→ Qu + c + τ −1 (u − v) + NK (u) is µτ -strongly monotone and contains 0 at p(v), ∥e x − p(v)∥ ≤ µ−1 τ ∥r∥.

(22) p Alternatively, an objective gap Hv (e x) − Hv (p(v)) ≤ ∆ and µτ -strong convexity give ∥e x − p(v)∥ ≤ 2∆/µτ . Any feasible inner method whose certificates obey a geometric or relative budget therefore satisfies Theorem 2.5, regardless of how the inexact point is produced. 2.3.2

An illustrative example: Projected Gradient

Projected gradient (PG) is one convenient way to produce the required inexact proximal point, but it is not required by the preceding theorem. Choose ηPG :=

2 , Lτ + µτ

ρPG :=

8

Lτ − µτ ∈ [0, 1). Lτ + µτ

(23)

P REPRINT. U NDER REVIEW.

Proposition 2.6 (fixed projected-gradient work). Suppose each proximal subproblem is warm-started at un,k,0 = x en,k and solved by J PG steps with stepsize ηPG . There is a local constant Cw > 0, independent of J ≥ 1, such that εn ≤ T Cw ρJPG dn .

(24)

Consequently, any fixed J satisfying T Cw ρJPG < 1 gives local Q-linear convergence. The proof is deferred to Appendix A.3. Proposition 2.6 suggests that we should choose J in the order of logρPG T −1 =

log T log T .  −1 ≤ log ρPG log 1 + 2 τ ∥Q∥

Thus, a logarithmic number of the restart length T is sufficient to reduce the proximal error to the accuracy required by the convergence analysis. In particular, the inner work grows only mildly with T , while its constant depends on the contraction factor ρPG , or equivalently on the conditioning of the quadratic proximal subproblem through τ |Q|.

3

Uniform Local Quadratic Growth

Section 2 showed that uniform local quadratic growth of the smoothed gap is sufficient for local linear convergence. This section characterizes that property through explicit primal and dual error bounds. The proof first establishes the product geometry of the KKT set and an exact decomposition of the smoothed gap, which identifies separate primal and dual error-bound certificates. It then uses a RSOC lifting to verify the primal error bound and a weighted projection residual to verify the dual error bound. Combining these two bounds yields uniform local quadratic growth. The lifting is used only to justify a primal error bound. Algorithm 1 continues to operate directly on the original conic QP (1). All proofs of this section are deferred to Appendix B. Throughout this section, ξ > 0 is fixed. 3.1

Exact smoothed-gap decomposition

For a fixed x⋆ ∈ X ⋆ , introduce the affine slack set and the complementary face DQ := Qx⋆ + c − range A∗ ,

Fx⋆ := K ∗ ∩ (x⋆ )⊥ .

Lemma 3.1 (geometry of the KKT solution set). Every x⋆ ∈ X ⋆ can be paired with every y ⋆ ∈ Y ⋆ . Consequently, Z ⋆ = X ⋆ × Y ⋆,

dist2 ((x, y), Z ⋆ ) = dist2 (x, X ⋆ ) + dist2 (y, Y ⋆ ).

(25)

Then DQ is independent of the chosen x⋆ , and the optimal slack set is S ⋆ = DQ ∩ Fx⋆ .

(26)

Fix a KKT center z ⋆ = (x⋆ , y ⋆ ) and set s⋆ = Qx⋆ + c − A∗ y ⋆ . Define Pξ (x; z ⋆ ) := Gξ ((x, y ⋆ ); z ⋆ ),

Dξ (y; z ⋆ ) := Gξ ((x⋆ , y); z ⋆ ).

These are the primal and dual components of the standard smoothed-gap decomposition [36, Lemma 5], i.e., Gξ (z; z ⋆ ) = Pξ (x; z ⋆ ) + Dξ (y; z ⋆ ). The exact identities needed below are also derived in the proof of Theorem 3.4. Fix a reference point z̄ = (x̄, ȳ) ∈ Z ⋆ . All local error bounds and regularity assumptions below are imposed near this reference point, and their constants are required to be uniform over all KKT centers in the local solution stratum. For the primal component, define the complementarity-curvature error EP (x; z ⋆ ) :=

1 ∥x − x⋆ ∥2Q + ⟨s⋆ , x⟩. 2

Both terms are nonnegative for x ∈ K.

9

(27)

P REPRINT. U NDER REVIEW.

Definition 3.2 (uniform local primal error bound). We say that the primal error bound holds uniformly locally around z̄ if there are neighborhoods VP of z̄ and Vx of x̄, and constants aP , bP > 0, such that, for every z ⋆ = (x⋆ , y ⋆ ) ∈ Z ⋆ ∩ VP and every x ∈ K ∩ Vx , dist2 (x, X ⋆ ) ≤ aP EP (x; z ⋆ ) + bP ∥Ax − b∥2 . (28) For the dual component, set Mξ := Q + ξI ≻ 0,

s⋆ (y) := Qx⋆ + c − A∗ y,

and let

1 M ΠK ξ (v) := arg min ∥u − v∥2Mξ . u∈K 2 The corresponding projection point and residual are   M uy := ΠK ξ x⋆ − Mξ−1 s⋆ (y) , ry := x⋆ − uy =: Rξ (y; z ⋆ ). Definition 3.3 (uniform local dual error bound). We say that the dual error bound holds uniformly locally around z̄ if there are neighborhoods VD of z̄ and Vy of ȳ, and a constant κD > 0, such that, for every z ⋆ = (x⋆ , y ⋆ ) ∈ Z ⋆ ∩ VD and every y ∈ Y ∩ Vy , dist(y, Y ⋆ ) ≤ κD ∥Rξ (y; z ⋆ )∥Mξ . (29) Theorem 3.4 (characterization of uniform local quadratic growth). Uniform local quadratic growth of the smoothed gap holds around z̄; that is, there are a neighborhood V of z̄ and a constant αξ > 0 such that Gξ (z; z ⋆ ) ≥ αξ dist2 (z, Z ⋆ )

(30)

for every z ∈ (K × Y ) ∩ V and z ⋆ ∈ Z ⋆ ∩ V , if and only if the uniform local primal and dual error bounds in Definitions 3.2 and 3.3 both hold around z̄. Theorem 3.4 reduces uniform local quadratic growth exactly to the two error bounds. The following subsection verifies them under the corresponding one-sided strict-complementarity conditions. 3.2

Primal and dual error bound under strict complementarity

Primal and dual strict complementarity. Let (x⋆ , s⋆ ) ∈ K × K ∗ be a complementary pair. We say that it satisfies, respectively,
primal strict complementarity (P-SC) : x⋆ ∈ ri K ∩ (s⋆ )⊥ ,
dual strict complementarity (D-SC) : s⋆ ∈ ri K ∗ ∩ (x⋆ )⊥ . Two-sided, or facial, strict complementarity means that both conditions hold; for product cones, these definitions are understood blockwise. Our terminology names the component required to lie in a relative interior. The terminology of [15, Definition 2] is reversed: their “dual strict complementarity” is P-SC above. We therefore state the relevant side explicitly whenever that result is invoked. For K = Rn+ , complementarity makes P-SC and D-SC equivalent to x⋆i + s⋆i > 0 for every i. The Goldman-Tucker theorem [22] guarantees the existence of such an optimal pair for LP whenever both optimal sets are nonempty, but not that every optimal pair is strictly complementary. More generally, D-SC implies P-SC when the minimal face of K containing x⋆ is exposed; the reverse implication follows from the analogous property of K ∗ . Hence the two sides coincide when both K and K ∗ are facially exposed [11]. In particular, for symmetric cones they are equivalent to x⋆ + s⋆ ∈ int K. This covers nonnegative, second-order-cone (SOC) and rotated-SOC, and PSD cones. For cones with nonexposed faces, P-SC and D-SC may differ; the exponential cone is a standard example [31]. P-SC and D-SC concern a particular complementary pair. The uniform bounds below require the relevant condition over a local certificate family: the primal error bound is verified under P-SC at every local KKT center, whereas Corollary 3.9 requires each local primal center to admit a D-SC optimal slack. These conditions are distinct from primal and dual nondegeneracy [1]. 10

P REPRINT. U NDER REVIEW.

3.2.1

Primal error bound from a rotated-SOC lifting

Rotated-SOC lifting. We verify the primal error bound in Definition 3.2 through an epigraph lifting. Choose B such that Q = B ∗ B, and define the rotated SOC as Qr := {(u, v, w) : u ≥ 0, v ≥ 0, 2uv ≥ ∥w∥2 }. Since u ≥ 12 ∥Bx∥2 is equivalent to (u, 1, Bx) ∈ Qr , problem (1) is equivalent to min

u + ⟨c, x⟩

s.t.

Ax = b,

x ∈ K,

v = 1,

w = Bx,

x,u,v,w

(31)

(u, v, w) ∈ Qr . Write χ := (x, u, v, w), denote the optimal solution set of (31) by XL⋆ , and collect the lifted equality constraints as AL (x, u, v, w) := (Ax, v, w − Bx),

eb := (b, 1, 0).

At every optimum the epigraph inequality is tight. Attaching multipliers (y, α, t) to the equality residuals b − Ax, 1 − v, and Bx − w shows that the lifted dual slacks have the form p = (1, −α, −t) ∈ Qr , with lifted dual objective ⟨b, y⟩ + α. Finally, for every x ∈ K, define the lifted comparison point   1 2 q(x) := ∥Bx∥ , 1, Bx , χ b(x) := (x, q(x)). 2 Transfer of strict complementarity. For any KKT center z ⋆ = (x⋆ , y ⋆ ), set s⋆ := Qx⋆ + c − A∗ y ⋆ , and define the canonical lift 1 q ⋆ := (u⋆ , 1, Bx⋆ ), t⋆ := Bx⋆ , α⋆ := −u⋆ , p⋆ := (1, u⋆ , −Bx⋆ ), u⋆ := ∥Bx⋆ ∥2 , 2 together with χ⋆ := (x⋆ , q ⋆ ), λ⋆ := (y ⋆ , α⋆ , t⋆ ), ζ ⋆ := (s⋆ , p⋆ ). The point χ⋆ is primal feasible, ζ ⋆ is the lifted slack at λ⋆ and is dual feasible, and ⟨χ⋆ , ζ ⋆ ⟩ = ⟨x⋆ , s⋆ ⟩ + 2u⋆ − ∥Bx⋆ ∥2 = 0, so (χ⋆ , λ⋆ , ζ ⋆ ) is a lifted KKT certificate. Besides, the added RSOC block is always strictly complementary because q ⋆ + p⋆ = (u⋆ + 1, u⋆ + 1, 0) ∈ ri Qr . Consequently, the lifted center satisfies P-SC exactly when the original center does since


χ⋆ = (x⋆ , q ⋆ ) ∈ ri K ∩ (s⋆ )⊥ × ri Qr ∩ (p⋆ )⊥ = ri (K × Qr ) ∩ (s⋆ , p⋆ )⊥ , where the last equality holds because relative interiors commute with products. Similarly, D-SC also transfers blockwise between the original and the lifted certificates. The verification is organized around the following residual error bound for the lifted problem. Assumption 3.5 (uniform lifted residual error bound). There are a neighborhood VSC of z̄, a common open neighborhood N of χ b(x̄), and constants κL , γL > 0 such that, for every z ⋆ ∈ Z ⋆ ∩ VSC with canonical lifted slack ζ ⋆ = (s⋆ , p⋆ ), dist(χ, XL⋆ ) ≤ κL ⟨ζ ⋆ , χ⟩1/2 + γL ∥AL χ − eb∥, χ ∈ (K × Qr ) ∩ N . (32)

11

P REPRINT. U NDER REVIEW.

Proposition 3.6 (verification of the lifted residual error bound). Let K be a finite Cartesian product of nonnegative cones, second-order cones (including rotated second-order cones), and positive-semidefinite cones. If every z ⋆ ∈ Z ⋆ ∩ VSC satisfies P-SC, for some neighborhood VSC of z̄, then Assumption 3.5 holds. Proposition 3.6 covers the rotated block Qr because the rotated second-order cone is the image of the standard √ √
second-order cone {(t, ξ) : t ≥ ∥ξ∥} under the orthogonal map (t, ξ1 , ξ2 ) 7→ (t + ξ1 )/ 2, (t − ξ1 )/ 2, ξ2 , a 45◦ rotation of the (t, ξ1 )-plane. Orthogonal maps preserve distances, faces, and relative interiors, so properties of second-order cones—in particular the error bounds established in [15]—hold verbatim, with identical constants, for rotated second-order-cone blocks. The corollary below therefore only needs to compare the original and lifted solution-set distances. Corollary 3.7 (verification of the uniform local primal error bound). Under Assumption 3.5, the uniform local primal error bound in Definition 3.2 holds with aP = 2κ2L , bP = 2γL2 . The linear lifted feasibility residual in (32) is essential: a square-root residual would yield ∥Ax − b∥, rather than ∥Ax − b∥2 , after squaring and would not imply Definition 3.2. 3.2.2

Dual error bound

We next verify the dual error bound in Definition 3.3. The weighted projection defining Rξ (y; z ⋆ ) yields vy := Mξ ry − s⋆ (y) ∈ NK (uy ). Normal-cone calmness and vy + s⋆ (y) = Mξ ry control the distance to the complementary face; linear regularity then gives the distance to the optimal slack set, and a pseudoinverse of A∗ lifts this estimate to the multiplier space. The key identity is S ⋆ = DQ ∩ Fx⋆ so we first state the two local regularity properties used in this chain. Linear regularity. Let C1 and C2 be nonempty closed convex sets and let s̄ ∈ C1 ∩ C2 . The pair {C1 , C2 } is locally linearly regular at s̄ if there are a neighborhood W of s̄ and a constant κ > 0 such that dist(s, C1 ∩ C2 ) ≤ κ [dist(s, C1 ) + dist(s, C2 )] ,

s ∈ W.

(33)

Uniform local linear regularity means that the neighborhood and modulus can be chosen independently of the parameter. Applied to the two sets defining S ⋆ , it gives the following assumption. Assumption 3.8 (uniform local slack regularity). Set s̄ := Qx̄ + c − A∗ ȳ. There are neighborhoods VA of z̄ and WA of s̄, and a constant κA > 0, such that, for every z ⋆ = (x⋆ , y ⋆ ) ∈ Z ⋆ ∩ VA and every s ∈ WA , dist(s, S ⋆ ) ≤ κA [dist(s, DQ ) + dist(s, Fx⋆ )] .

(34)

This is the uniform local linear regularity of {DQ , Fx⋆ } along the local KKT stratum. Since s⋆ (y) ∈ DQ , the first distance in (34) vanishes along the residual path. D-SC gives a direct verification of this property. Corollary 3.9 (dual strict complementarity implies slack regularity). Fix a KKT center z ⋆ = (x⋆ , y ⋆ ). If its optimal slack satisfies D-SC, s⋆ ∈ ri Fx⋆ , (35) then {DQ , Fx⋆ } is locally linearly regular at s⋆ , so the corresponding pointwise slack-regularity bound holds. If every x⋆ ∈ X ⋆ ∩ VX,SC , for some neighborhood VX,SC of x̄, admits a D-SC optimal slack, then Assumption 3.8 holds uniformly on the corresponding local KKT stratum. We next impose the normal-cone calmness used in this projection estimate.

12

P REPRINT. U NDER REVIEW.

Assumption 3.10 (uniform normal-cone calmness). There are neighborhoods VB of z̄ and WB of (x̄, −s̄), and a constant κB > 0, such that, for every z ⋆ = (x⋆ , y ⋆ ) ∈ Z ⋆ ∩ VB and every (u, v) ∈ gph NK ∩ WB , dist(v, NK (x⋆ )) ≤ κB ∥u − x⋆ ∥.

(36)

A standard pointwise sufficient condition is C 2 -cone reducibility. A closed convex set Ω is C 2 -cone reducible at x̄ ∈ Ω if, locally, it has the representation Ω = {x : Ξ(x) ∈ C}, where Ξ is twice continuously differentiable, Ξ(x̄) = 0, DΞ(x̄) is surjective, and C is a pointed closed convex cone. If K is C 2 -cone reducible at x⋆ , then NK is calm at every normal graph point, including (x⋆ , −s⋆ ) [33, Theorem 2.1]. This class includes polyhedral, Lorentz, positive-semidefinite, and Ky Fan k-norm epigraph cones, and is preserved by finite products. For polyhedral cones, calmness also follows directly from piecewise polyhedrality [50]; for exponential and power cones it must be verified separately unless a suitable reducibility result is available. Table 2 summarizes these pointwise verification routes; rotated SOCs are included with Lorentz cones. Table 2: Typical pointwise verification routes for the two dual regularity conditions. Cone class

Slack-set linear regularity

Normal-cone calmness

Polyhedral cones, including Rn + Lorentz/SOC and rotated SOC Positive-semidefinite cone Finite Cartesian products General closed convex cone

Hoffman bound; often global.

Polyhedral multifunction [50].

D-SC: s⋆ ∈ ri Fx⋆ .

C 2 -cone reducibility [33].

D-SC, equivalently the rank condition.

C 2 -cone reducibility [33].

Blockwise D-SC. Verify linear regularity directly.

Reducibility preserved by products [33]. Verify calmness directly.

The table gives pointwise conditions. If each local primal center admits a D-SC optimal slack, Corollary 3.9 makes slack regularity uniform. The next proposition verifies uniform normal-cone calmness in the same regime for the common cone classes. Proposition 3.11 (verification of uniform normal-cone calmness). Let K be a finite Cartesian product of nonnegative cones, second-order cones (including rotated second-order cones), and positive-semidefinite cones. If every x⋆ ∈ X ⋆ ∩ VX,SC , for some neighborhood VX,SC of x̄, admits a blockwise D-SC optimal slack, then Assumption 3.10 holds. Without the constant-face structure provided by D-SC, the calmness modulus can blow up as centers approach a smaller face, so the stratification is what makes uniformity available. We now combine the two properties. Proposition 3.12 (verification of the uniform local dual error bound). Under Assumptions 3.8 and 3.10, the uniform local dual error bound in Definition 3.3 holds with constant Lξ + κB κD = ∥(A∗ )† ∥ κA p , λξ where Lξ := ∥Mξ ∥ = ∥Q∥ + ξ, λξ := λmin (Mξ ) ≥ ξ, and the pseudoinverse norm is restricted to range A∗ . 3.3

Quadratic growth for common-cone QPs under strict complementarity

The preceding verification routes yield the following direct consequence for the cone classes most commonly used in conic optimization. The uniformity clause is stated explicitly because the corresponding result at a single strictly complementary KKT point is only pointwise. Corollary 3.13 (common-cone conic QPs). Fix ξ > 0, and suppose that K is a finite Cartesian product of nonnegative cones, Lorentz cones (including rotated Lorentz cones), and positive-semidefinite cones. Let z̄ ∈ Z ⋆ , and assume strict 13

P REPRINT. U NDER REVIEW.

complementarity along the local KKT stratum: there is a neighborhood VSC of z̄ such that every z ⋆ = (x⋆ , y ⋆ ) ∈ Z ⋆ ∩ VSC satisfies x⋆ + s⋆ ∈ int K. Then, there are a neighborhood V of z̄ and a constant αξ > 0 such that Gξ (z; z ⋆ ) ≥ αξ dist2 (z, Z ⋆ ) for all z ∈ (K × Y ) ∩ V and z ⋆ ∈ Z ⋆ ∩ V . With the constants κL and γL from Assumption 3.5, κA from Assumption 3.8, and κB from Assumption 3.10, one may take   1 1 ξ αξ = min , , . (37) 2κ2L 4ξγL2 2∥(A∗ )† ∥2 κ2A (∥Q∥ + ξ + κB )2 Theorem 3.4 verifies Assumption 2.2; it is not a separate convergence argument. The exact averaged PDHG output lies in K × Y because every exact primal iterate lies in K and K is convex. Therefore the theorem applies directly to the epoch average whenever it lies in the stated neighborhood.

4

Practical Implementation

This section describes the numerical components used to implement and accelerate the PDHG iteration of Algorithm 1. We first discuss the primal proximal problem, then present the reflected-Halpern scheme, problem rescaling, inner-solver enhancements, adaptive restart, and the primal-weight update. Finally, we describe the distributed implementation for multi-GPU acceleration. The implementation accepts problems in the following standard form: 1 ⊤ x Qx + c⊤ x 2 s.t. la ≤ Ax ≤ ua ,

min

x∈Rn

(38)

xI0 ∈ X0 := [lx , ux ], xIj ∈ Kj , j = 1, . . . , J, QJ Here A ∈ Rm×n , Q ∈ Sn+ , c ∈ Rn , and X := X0 × j=1 Kj . Let R := R ∪ {−∞, +∞}. The affine bounds satisfy |I0 |

m

la , ua ∈ R with la ≤ ua , and the box bounds satisfy lx , ux ∈ R with lx ≤ ux . Thus equal finite endpoints represent equalities, while infinite endpoints represent one-sided or absent bounds. The index sets {Ij }Jj=0 form a partition of {1, . . . , n}, so every primal coordinate belongs to exactly one block. The implementation supports nonnegative blocks and the vector cones Ksoc := {(v, t) : ∥v∥2 ≤ t}, Krsoc := {(v, s, t) : ∥v∥22 ≤ 2st, s, t ≥ 0}, Kexp := cl{(r, s, t) : s > 0, s exp(r/s) ≤ t}, α Kpow := {(r, s, t) : r, s ≥ 0, rα s1−α ≥ |t|},

α ∈ (0, 1).

PSD cone blocks are excluded because every projection would require an eigenvalue decomposition. For a closed convex set C and a positive diagonal matrix M ≻ 0, define the projection in the M -norm by 1 2 ΠM C (r) := arg min ∥u − r∥M , u∈C 2

∥z∥2M := ⟨z, M z⟩.

(39)

We write ΠC := ΠIC for the Euclidean projection. Because X is a Cartesian product and M is diagonal, ΠM X decomposes into independent box and cone projections. Table 12 summarizes the projection methods for the box and cone blocks, and Appendix C gives the corresponding formulas under the Euclidean norm and the M -norm.

14

P REPRINT. U NDER REVIEW.

4.1

Primal Proximal Problem Solving

Each PDHG iteration requires solving the following strongly convex primal proximal subproblem, where τn is the current stepsize and v := x + τn A∗ y:   1 1 p+ = arg min ⟨u, Qu⟩ + ⟨c, u⟩ + ∥u − v∥2 , (40) u∈X 2 2τn whose smooth Hessian Q + τn−1 I ≻ 0 makes the solution unique even for singular Q. The solution method depends on whether Q is diagonal or a general positive-semidefinite operator. General Quadratic Objective. For a general matrix-free Q, we solve (40) by a projected Barzilai-Borwein method [12]. Let g j := Quj + c + τn−1 (uj − v). (41) For j ≥ 1, define sj := uj − uj−1 and δgj := g j − g j−1 . The BB1 stepsize and projected update are αjBB :=

∥sj ∥22 , ⟨sj , δgj ⟩

uj+1 = ΠX uj − αjBB g j .

(42)

Each inner iteration requires one application of Q and one product-set projection. Products with a structured operator may be evaluated as Q = P + R⊤ DR, where D is symmetric and Q ⪰ 0, without materializing Q. Diagonal Quadratic Objective. Suppose Q = diag(q) with qi ≥ 0. More generally, let Msc ≻ 0 denote the positive diagonal matrix induced by the cone scaling described in Section 4.2.2; its cone-block components are defined in (51). When no cone scaling is applied, Msc = I. The resulting projection metric and center are
Mn := Msc + τn Q, x e := Mn−1 Msc x − τn (c − A∗ y) . (43) The proximal subproblem is therefore the single metric projection n p+ = ΠM x). X (e

(44)

Thus diagonal quadratic curvature and diagonal cone scaling are handled by the same projection definition (39); no inner projected-gradient iteration is required. 4.2

Algorithm Enhancement

The following enhancements improve practical solution efficiency through reflected-Halpern acceleration, problem rescaling, reduced inner work, timely restarts, and primal-dual balance. 4.2.1

Reflected-Halpern Acceleration

The practical solver applies an anchored reflected-Halpern acceleration to the PDHG operator. This acceleration has demonstrated practical effectiveness in cuPDLPx [39]; its convergence and nonergodic O(1/k) residual bounds in a more general setting are established in [32]. During epoch n, let Tn denote the PDHG update operator with fixed τn , σn , set the anchor an = z n,0 , and store the candidate wn,k+1 := Tn (z n,k ).

(45)

For a reflection coefficient ρ ∈ [1/2, 1] and k = 0, . . . , Nn − 2, the anchored update is z n,k+1 =

k+1 1 2ρwn,k+1 + (1 − 2ρ)z n,k + an . k+2 k+2

(46)

The choice ρ = 1/2 gives standard Halpern iteration on Tn , whereas ρ = 1 gives Halpern iteration on the reflected operator 2Tn − I. A reflected-Halpern state need not have a primal component in the cone. Consequently, KKT residuals

15

P REPRINT. U NDER REVIEW.

are evaluated at the stored PDHG candidates, and an accepted epoch restarts from the final candidate z n+1,0 := wn,Nn = Tn (z n,Nn −1 ),

(47)

or from its inexact counterpart in the implementation. Section 2 analyzes the restarted averaged scheme with inexact primal proximal evaluations; the reflected-Halpern recursion is used here as an implementation enhancement built on the same PDHG operator. 4.2.2

Problem Rescaling

Following the preconditioned-instance convention of cuPDLPx [39], we combine ℓ∞ Ruiz rescaling [51], PockChambolle rescaling [49], and bound-objective rescaling. To describe the accumulated diagonal transformation, let Sx ≻ 0 and Sy ≻ 0 be the variable and constraint scalings and set x e = Sx x. The main problem data become e = Sy AS −1 , A x

e = S −⊤ QS −1 , Q x x

e c = Sx−⊤ c,

(48)

with all bounds transformed consistently. Bound-objective rescaling adds positive scalar normalizations of the bound and objective data. For primal cone blocks, we support two treatments of the variable scaling. Let Ij index a cone block and, at the current rescaling stage, let the coordinatewise Ruiz and Pock-Chambolle candidates be !1/2  1/2 X 2−α eℓi | eℓi | ri := max |A , pi := |A . (49) ℓ

ℓ

The cone-preserving block scaling mode replaces these candidates on Ij by  1/2 X 1 rIj := max ri , pIj :=  p2i  , i∈Ij |Ij |

(50)

i∈Ij

respectively. Thus Ruiz retains its blockwise ℓ∞ safeguard, whereas Pock-Chambolle aggregates the column energies by their root mean square. All coordinates of a cone block receive the resulting common factor, so Sj = dj I and Sj Kj = Kj . Nonconic variable scalings and all constraint scalings remain coordinatewise. Alternatively, the rescaled-cone mode retains the candidates in (49) coordinatewise. A cone block then becomes Sj Kj , where Sj = diag((Sx )Ij ). This applies the Ruiz and Pock-Chambolle scalings without blockwise aggregation and invokes the rescaled-cone projections described by [30]. Equivalently, in the original cone coordinates, M

Sj−1 ΠSj Kj (Sj r) = ΠKjsc,j (r),

Msc,j := Sj⊤ Sj .

(51)

Constraint scalings use the same coordinatewise rule in both modes. As in cuPDLPx, termination criteria are evaluated after recovering the candidate for the original, unpreconditioned problem, so rescaling does not affect the reported accuracy. The default rescaling settings are reported in Table 3. 4.2.3

Adaptive Inner Accuracy

For a general Q, let ∆un,k denote the final projected update of the inner solve at outer step (n, k). For k ≥ 1, the inner tolerance is    ∥xn,k − xn,k−1 ∥ in ϵin = min ϵ , max γ , ϵ , (52) in min n,k n,k−1 τn and the recorded stopping test is ∥∆un,k ∥2 ≤ ϵin (53) n,k . By projection optimality, the final update supplies a stationarity residual proportional to ∥∆un,k ∥2 , so (22) converts this displacement into a bound on the proximal error of the returned point. The first tolerance of an epoch is inherited from the preceding epoch, with ϵin 0 used at initialization. 16

P REPRINT. U NDER REVIEW.

4.2.4

Jacobi Inner Preconditioning

For a general Q, we accelerate the Euclidean projected-BB method in (42) with the Jacobi metric
Jn := diag Q + τn−1 I = diag(Q) + τn−1 I.

(54)

The corresponding diagonally scaled projected-BB step is αjJBB =

⟨sj , Jn sj ⟩ , ⟨sj , δgj ⟩

uj+1 = ΠJXn uj − αjJBB Jn−1 g j ,

(55)

where g j denotes the gradient of the smooth part of (40) at uj , sj := uj − uj−1 , and δgj := g j − g j−1 . This is a scaled gradient-projection step with a spectral steplength in the sense of [4]. Since Jn is diagonal, its inverse is applied elementwise, and the metric projection is evaluated using the formulas in Appendix C. The Jacobi metric is used only to define the inner direction and projection; the stopping test is restored to the c n,k := Jn1/2 ∆un,k , unpreconditioned Euclidean coordinates. Indeed, for ∆u c n,k ∥2 = ∥∆un,k ∥J , ∥∆u n

c n,k ∥2 . ∥∆un,k ∥2 = ∥Jn−1/2 ∆u

(56)

The second quantity is the one used in (53). Appendix E reports an ablation of the adaptive inner-accuracy rule and Jacobi inner preconditioning. 4.2.5

Adaptive Stepsize and Restart

Following cuPDLPx [39], we define the primal and dual stepsizes using a step size η and a primal weight ωn : τn =

η , ωn

σn = ηωn ,

η :=

0.998 . ∥A∥2

(57)

The primal weight ωn is updated when a restart is triggered. Primal Weight Update. For the accepted endpoint wn,Nn , the proportional-integral-derivative (PID) update uses the movement ratio and exponentially discounted integral ! ∥wyn,Nn − y n,0 ∥ , In = βI In−1 + en . (58) en = log ωn ∥wxn,Nn − xn,0 ∥ We update the primal weight by log ωn+1 = log ωn + KP en + KI In + KD (en − en−1 ). Adaptive Restart.

(59)

For the stored PDHG candidate, define n,k dn,k − wxn,k+1 , x := x

dn,k := y n,k − wyn,k+1 . y

The fixed-point error is 2 2 −1 n,k 2 n,k n,k En,k = ωn ∥dn,k x ∥ + ωn ∥dy ∥ + 2η⟨Adx , dy ⟩.

(60)

Using our thresholds, a restart occurs under sufficient decay, necessary decay with no local progress, or an artificial iteration limit: En,k ≤ βs En,0 , or En,k ≤ βn En,0 and En,k > En,k−1 ,

or

k ≥ βa Ktot , where Ktot is the cumulative number of outer iterations. Table 3 gives the default constants.

17

(61)

P REPRINT. U NDER REVIEW.

Table 3: Default parameters for problem rescaling and algorithm enhancements.

4.3

Component

Parameters

Default values

Ruiz rescaling Pock-Chambolle rescaling Cone scaling Bound-objective rescaling

Iterations α Mode Status

10 1.00 Cone-preserving Enabled

Inner accuracy Jacobi preconditioning Reflected-Halpern acceleration Restart Primal weight

(ϵin 0 , γin , ϵmin ) Status ρ (βs , βn , βa ) (βI , KP , KI , KD )

(10−3 , 5×10−4 , 10−9 ) Enabled 1.00 (0.20, 0.80, 0.36) (0.30, 0.99, 0.01, 0.00)

Multi-GPU Implementation

PDHCG-CQP distributes the products with A, A⊤ , and Q across multiple GPUs. The constraint operator uses the two-dimensional partition of D-PDLP [28], while the quadratic operator follows the same primal-vector partition. This increases aggregate memory capacity without fully replicating the primal and dual vectors. 4.3.1

Distributed Storage

The constraint matrix A is sharded into a Pr × Pc grid of blocks. Writing   A11 · · · A1Pc  . ..  x = (x1 , . . . , xPc ), A =  .. . , APr 1

···

y = (y1 , . . . , yPr ),

APr Pc

rank (i, j) stores Aij , the corresponding slices xj and yi , and their local work vectors. Thus xj is shared only within process column j and yi only within process row i. The quadratic operator uses the same partition x = (x1 , . . . , xPc ). For the structured representation Q = P + R⊤ DR, write h i h i P = P:1 · · · P:Pc , R = R1 · · · RPc , (62) where P:j and Rj contain the columns multiplying xj . Rank (i, j) stores these two column shards; because they depend only on j, they are replicated within process column j. The middle matrix D and the rank-dimensional work vectors are replicated. For diagonal Q, rank (i, j) stores only the local diagonal block qj . 4.3.2

Distributed Computation

The two constraint products are assembled from local sparse products as (Ax)i =

Pc X

⊤

Aij xj ,

(A y)j =

j=1

Pr X

A⊤ ij yi .

(63)

i=1

The first sum is reduced across each process row and the second across each process column. For the partition in (62), the quadratic product is assembled as p :=

Pc X

P:ℓ xℓ ,

ℓ=1

h :=

Pc X

Rℓ xℓ ,

ℓ=1

(Qx)j = pj + Rj⊤ Dh,

(64)

j = 1, . . . , Pc .

The sums defining p and h are reduced across each process row; every rank then extracts the jth primal slice pj and applies Rj⊤ . The sparse and low-rank terms may be used independently by setting the other term to zero. For diagonal 18

P REPRINT. U NDER REVIEW.

Q, the product is local, (Qx)j = qj ⊙ xj . Global scalar quantities in Section 4.2 are obtained by summing local inner-product and squared-norm contributions. The Fisher study in Section 5.6 reports the resulting end-to-end scaling.

5

Numerical Experiments

We evaluate PDHCG-CQP on standard and large-scale QP, public convex QCQP, Mittelmann SOCP, and large-scale quasilinear Fisher equilibrium benchmarks. We study multi-GPU scaling on both large-scale QP and conic instances. 5.1

Experimental Setup

Selected baselines. We consider a broad set of candidate solvers spanning first-order and interior-point methods, and compare PDHCG-CQP with the selected baselines. The selected first-order baselines are PDQP [36], PDHCG [25], HPR-QP and HPR-SOCP [9], SCS [47, 48], PDCS [30], and OSQP [52]. The selected interior-point baselines are MOSEK [45], Gurobi [24], COPT [21], and Clarabel [23]. All solvers are compared at the same target accuracy ϵ. Computing environment. Experiments are conducted on a server with eight NVIDIA H100 GPUs, each with 80GB HBM3, an Intel Xeon Platinum 8469C CPU at 2.60GHz, and 512GB RAM. Unless otherwise stated, each run uses a single GPU. Evaluation metrics. We report the number of solved instances and shifted geometric mean (SGM) runtime with shift 10; the QP tables also report arithmetic mean runtime: ! N 1 X log(ti + 10) − 10. (65) SGM10 = exp N i=1 Any unsuccessful run is charged the time limit when it is included in an aggregate runtime metric. 5.2

Standard Convex QP Benchmarks

We first evaluate the QP specialization of PDHCG-CQP on two standard benchmark repositories. • Maros-Mészáros benchmark [41]: 134 convex quadratic programming instances, tested at tolerance ϵ = 10−6 with a 1000 second time limit. • Mittelmann benchmark [43]: 21 continuous convex QP instances, tested at tolerances ϵ = 10−6 and ϵ = 10−8 with a 3600 second time limit. Table 4: Performance on 134 Maros-Mészáros instances with a 1000-second time limit and tolerance ϵ = 10−6 . The best results are marked in bold, and the second best are underlined. Solver PDQP PDHCG HPR-QP PDHCG-CQP

Total

Solved

SGM10 (s)

Average (s)

134

118 111 124 126

28.53 33.51 10.56 10.95

160.70 210.53 93.32 83.49

Table 4 shows that PDHCG-CQP improves robustness over the earlier PDHCG implementation and is competitive with HPR-QP: with Jacobi preconditioning on the non-diagonal subset it solves 126 instances, attains the best arithmetic mean runtime, and is second best in shifted geometric mean. On the Mittelmann benchmark, PDHCG-CQP solves 17 of the 21 instances at both tolerances, the most among the tested solvers, and tightening the tolerance from 10−6 to 10−8 increases its runtime only moderately; the aggregate comparison does not isolate the individual effects of reflection, restart, and inner accuracy.

19

P REPRINT. U NDER REVIEW.

Table 5: Performance on 21 Mittelmann instances with a 3600-second time limit. The best results are marked in bold, and the second best are underlined. Solver

Total

Solved

SGM10 (s)

Average (s)

−6

Tolerance ϵ = 10 PDQP PDHCG HPR-QP PDHCG-CQP

21

11 10 14 17

410.34 486.81 126.88 72.74

1942.89 1987.11 1282.34 766.94

Tolerance ϵ = 10−8 HPR-QP PDHCG-CQP

5.3

21

14 17

137.42 91.22

1318.61 894.32

Public Convex QCQP Benchmark

We next consider 21 public convex QCQP instances assembled from QPLIB [19] and the Mittelmann benchmark collection [43]. Write an instance as min

x∈X

s.t.

1 ⊤ x Q0 x + c⊤ 0x 2 1 ⊤ x Qi x + c⊤ i x ≤ di , 2

(66) i = 1, . . . , p,

where X collects the linear constraints and variable bounds and Qi ⪰ 0. Choosing factors Qi = Ri⊤ Ri and using the rotated second-order cone from Section 4.1, each quadratic constraint is equivalent to (di − c⊤ i x, 1, Ri x) ∈ Krsoc .

(67)

Retaining the quadratic objective in (66) and replacing only the quadratic constraints by (67) gives a quadratic SOCP (QSOCP). A solver requiring a linear objective receives the equivalent SOCP min

t + c⊤ 0x

s.t.

(t, 1, R0 x) ∈ Krsoc ,

x∈X , t

(di − c⊤ i x, 1, Ri x) ∈ Krsoc ,

(68) i = 1, . . . , p.

Thus the SOCP path adds an objective-epigraph cone, whereas the QSOCP path retains Q0 explicitly. The formulation used by each solver is recorded in Table 6; COPT GPU is evaluated through its native QCQP interface. Each run has a 3600 second limit. The table reports the number of runs declared solved at the target tolerance and SGM10 over all 21 instances, with every unsuccessful run charged the time limit. PDHCG-CQP solves all 21 instances at 10−4 and 10−6 and 20 at 10−8 . It remains the most robust and fastest first-order solver in aggregate at every target tolerance. Against HPR-SOCP, PDHCG-CQP records 21 wins and no losses, 19 wins and 2 losses, and 17 wins, 3 losses, and 1 tie at 10−4 , 10−6 , and 10−8 , respectively. At 10−4 , PDHCG-CQP attains the second-lowest SGM10 overall, behind only Gurobi; the commercial QCQP interfaces retain an aggregate runtime advantage at the tighter tolerances. 5.4

Public SOCP Benchmark

The SOCP study uses the 18 public instances in the Mittelmann collection [44]. Both PDHCG-CQP and HPR-SOCP solve all 18 instances at both tolerances. HPR-SOCP has the lower aggregate SGM10 , although the two solvers are close at 10−4 , with values of 18.11 and 18.76. Moreover, PDHCG-CQP remains competitive at the instance level, with 8 wins and 10 losses at 10−4 and 5 wins and 13 losses at 10−6 . Thus the updated 20

P REPRINT. U NDER REVIEW.

Table 6: Aggregate performance on 21 public convex QCQP instances. The formulation column gives the actual input form used in each run. Times are in seconds and the time limit is 3600 seconds. Within each metric column, the best result is marked in bold and the second best is underlined; tied values share the same marking. ϵ = 10−4 Solver

Formulation

ϵ = 10−6

ϵ = 10−8

Solved

SGM10

Solved

SGM10

Solved

SGM10

First-order methods PDHCG-CQP QSOCP HPR-SOCP QSOCP SCS indirect SOCP SCS direct SOCP PDCS SOCP

21 19 13 16 13

8.11 36.00 209.18 62.70 152.94

21 19 10 14 13

21.64 54.56 1107.13 215.37 249.69

20 19 5 11 12

42.63 82.06 2744.60 625.47 664.11

Interior-point methods MOSEK QCQP Gurobi QCQP COPT GPU QCQP Clarabel SOCP

20 21 20 21

16.40 6.06 9.68 14.66

20 20 20 21

16.54 11.74 9.78 19.97

20 20 20 17

16.76 14.71 9.94 50.54

Table 7: Aggregate performance on 18 public SOCP instances. Times are in seconds and the time limit is 3600 seconds. Within each metric column, the best result is marked in bold and the second best is underlined; tied values share the same marking. ϵ = 10−4 Solver

Solved

ϵ = 10−6

SGM10

Solved

SGM10

18 14 14 4 18

18.76 613.98 475.45 2479.75 18.11

18 12 9 3 18

69.25 614.03 1166.90 2644.17 51.19

Interior-point methods COPT GPU 18 Gurobi 18 MOSEK 18 Clarabel 14

19.00 22.27 30.97 601.98

18 18 18 14

19.01 19.85 32.44 602.04

First-order methods PDHCG-CQP PDCS SCS direct SCS indirect HPR-SOCP

SOCP results show matching robustness but favor HPR-SOCP in aggregate runtime. In contrast, the QCQP results indicate that PDHCG-CQP benefits more clearly when the quadratic objective is retained and handled directly. 5.5

Real-World Large-Scale Sparse QPs

We also evaluate large-scale Lasso instances derived from LIBSVM [7] and the SuiteSparse Matrix Collection [13]. The Lasso problem min ∥Ax − b∥22 + λ∥x∥1 (69) x

is solved through its standard QP reformulation min y ⊤ y + λ1⊤ t x,y,t

s.t.

y = Ax − b, − t ≤ x ≤ t.

21

(70)

P REPRINT. U NDER REVIEW.

We set λ = 0.01∥A⊤ b∥∞ and use tolerance ϵ = 10−6 . The dimensions and densities of the nine instances are listed in Table 13 of Appendix D. Table 8: Single-GPU solve times on large-scale Lasso-derived QPs. Times are in seconds, with a 7200-second limit. The best results are marked in bold, the second best are underlined, and “f” denotes an unsuccessful run. First-order methods Problem

Interior-point

PDHCG-CQP

HPR-QP

PDHCG

PDQP

SCS GPU

OSQP

COPT

1.96 2.37 1337.05 217.70 1642.69 11.59 9.90 703.36 341.16

2.09 6.21 4911.41 753.65 3213.38 26.87 29.36 1971.24 842.02

3.35 7.12 1377.54 1429.55 4224.95 11.10 61.00 387.00 2705.94

7.30 19.54 5124.82 5557.97 f 46.49 148.01 1715.50 f

345.09 f f f f 490.66 f f f

80.32 f f f f 255.81 f f f

88.21 f f f f 69.57 f f f

SLS rcv1_test avazu-site.tr avazu-app avazu-site kddb2010_test kdda2010_test kddb2010_train kdda2010_train

Table 8 shows that PDHCG-CQP is the fastest solver on seven of the nine instances and is especially effective at the largest scales: on avazu-app it is 3.46× faster than HPR-QP, while the remaining solvers are slower still or fail. Table 9: Multi-GPU scaling of PDHCG-CQP on the Lasso-derived QPs. Times are in seconds, and Sp = t1 /tp is the speedup over one GPU. The fastest time in each row is marked in bold, and the second fastest is underlined. Solve time (s)

Speedup

Problem

1 GPU

2 GPUs

4 GPUs

8 GPUs

S2

S4

S8

SLS rcv1_test avazu-site.tr avazu-app avazu-site kddb2010_test kdda2010_test kddb2010_train kdda2010_train

1.96 2.37 1337.05 217.70 1642.69 11.59 9.90 703.36 341.16

1.10 1.48 953.02 116.92 993.62 11.18 8.50 466.28 237.19

0.56 1.06 547.28 65.86 765.46 8.79 6.77 317.69 167.62

0.50 0.81 376.90 46.30 503.99 7.07 5.61 233.44 124.89

1.78 1.60 1.40 1.86 1.65 1.04 1.16 1.51 1.44

3.50 2.24 2.44 3.31 2.15 1.32 1.46 2.21 2.04

3.92 2.93 3.55 4.70 3.26 1.64 1.76 3.01 2.73

All nine Lasso-derived QPs benefit from multi-GPU execution. On avazu-app, eight GPUs reduce the solve time from 217.70 to 46.30 seconds, a 4.70× speedup; the corresponding speedup on avazu-site.tr is 3.55×. 5.6

Large-Scale Fisher Equilibrium

Consider a quasilinear Fisher market with n buyers and m divisible goods. Buyer i has budget wi > 0, good j has supply bj > 0, and uij is buyer i’s value for good j. Let E := {(i, j) : uij > 0} be the sparse valuation graph. The quasilinear extension of the Eisenberg-Gale program [20] is   n n X X X max wi log uij xij + δi  − δi x,δ

s.t.

i=1

X

i=1

j:(i,j)∈E

xij = bj ,

j = 1, . . . , m,

i:(i,j)∈E

xij ≥ 0,

(i, j) ∈ E,

22

δi ≥ 0,

i = 1, . . . , n,

(71)

P REPRINT. U NDER REVIEW.

where δi represents the money retained by buyer i. Introducing zi and ti gives the equivalent exponential-cone formulation used in our experiments: n X (δi − wi zi )

min

x,δ,z,t

i=1

X

s.t.

xij = bj ,

j = 1, . . . , m,

i:(i,j)∈E

ti =

(72)

X

uij xij + δi ,

i = 1, . . . , n,

j:(i,j)∈E

(zi , 1, ti ) ∈ Kexp , xij ≥ 0,

i = 1, . . . , n,

(i, j) ∈ E,

δi ≥ 0,

i = 1, . . . , n.

Let e := |E| and ρ := e/(nm). We store allocations only on E and represent each (zi , 1, ti ) by an exponential-cone block whose middle coordinate is fixed at one. Thus (72) has e + 4n stored primal coordinates, m + n affine equations, n exponential cones, and e + n nonnegative coordinates. We evaluate five instances with 103 ≤ n ≤ 107 buyers, 400 ≤ m ≤ 4000 goods, and valuation density 0.01 ≤ ρ ≤ 0.20. Nonzero valuations uij and budgets wi are sampled uniformly from [0.1, 1.1], every good has supply bj = 0.2n, and all instances use random seed 1. The exact values of (n, m, ρ) and the resulting conic dimensions are reported in Table 14 of Appendix D. The solver time limit is 3600 seconds. Table 10: Solver comparison on the quasilinear Fisher equilibrium instances, with each GPU solver using one GPU. Times are in seconds. The best successful time in each row is marked in bold, the second best is underlined, and “f” denotes an unsuccessful or omitted run under the 3600-second resource limit. First-order methods n

PDHCG-CQP

PDCS

Interior-point methods

SCS direct

COPT GPU

MOSEK

Gurobi

Clarabel

0.67 6.61 64.69 f f

0.60 24.19 68.36 687.52 f

157.13 f f f f

1.30 80.82 1353.31 f f

0.50 7.16 66.10 f f

0.48 22.77 83.41 873.38 f

f f f f f

1.81 123.97 1929.71 f f

−4

Tolerance ϵ = 10 1,000 10,000 100,000 1,000,000 10,000,000

0.36 0.69 4.43 46.65 231.94

30.32 96.22 96.02 1565.44 f

7.45 356.51 f f f Tolerance ϵ = 10−6

1,000 10,000 100,000 1,000,000 10,000,000

0.85 0.77 4.77 48.32 241.53

30.96 502.63 1315.10 f f

238.16 f f f f

The symbol “f” covers time-limit, memory-limit, and numerical terminations. After a solver failed at a given size and tolerance, larger or tighter runs with the same solver were omitted. The smallest instances do not amortize collective latency, but scaling becomes effective from n = 105 : at n = 106 , eight GPUs reduce the solve time from 46.65 to 5.85 seconds at 10−4 and from 48.32 to 6.07 seconds at 10−6 , corresponding to speedups of 7.97× and 7.96×, and at n = 107 , PDHCG-CQP is the only solver reported as reaching optimality within the time and memory limits.

23

P REPRINT. U NDER REVIEW.

Table 11: Multi-GPU scaling of PDHCG-CQP on the quasilinear Fisher equilibrium instances. Times are in seconds, and Sp = t1 /tp is the speedup over one GPU. The fastest time in each row is marked in bold, and the second fastest is underlined. Solve time (s) n

1 GPU

2 GPUs

Speedup

4 GPUs

S2

S4

S8

0.49 0.55 1.05 5.85 34.28

0.86 1.15 1.81 2.04 1.81

0.80 1.38 3.03 4.18 3.52

0.73 1.25 4.22 7.97 6.77

1.11 0.60 1.12 6.07 35.69

0.88 1.17 1.83 2.04 1.81

0.84 1.38 3.08 4.19 3.52

0.77 1.28 4.26 7.96 6.77

8 GPUs −4

Tolerance ϵ = 10 1,000 10,000 100,000 1,000,000 10,000,000

0.36 0.69 4.43 46.65 231.94

0.42 0.60 2.45 22.87 128.38

0.45 0.50 1.46 11.16 65.93

Tolerance ϵ = 10−6 1,000 10,000 100,000 1,000,000 10,000,000

6

0.85 0.77 4.77 48.32 241.53

0.97 0.66 2.61 23.74 133.57

1.01 0.56 1.55 11.54 68.68

Conclusion

We introduced PDHCG-CQP, a matrix-free restarted primal-dual solver for large-scale conic convex quadratic programming on GPUs. The practical method combines primal-first PDHG updates with an anchored reflected-Halpern acceleration, matrix-free quadratic proximal solves, specialized cone projections. PDHCG-CQP supports both singleand multi-GPU execution. For the underlying restarted averaged scheme, we proved local linear convergence under uniform quadratic growth of the smoothed duality gap, with exact and inexact proximal evaluations, including a fixed amount of projected-gradient work per subproblem; the quadratic-growth property itself was verified under strict complementarity. The experiments cover standard and large-scale sparse QPs, public convex QCQP and SOCP benchmarks, and largescale quasilinear Fisher equilibrium instances. The Fisher study shows strong scaling through eight GPUs and robustness at 107 buyers, corresponding to more than 4.4 × 108 stored primal coordinates, and the Lasso-derived QPs confirm multi-GPU acceleration on large sparse diagonal-Q instances.

References [1] Farid Alizadeh, Jean-Pierre A. Haeberly, and Michael L. Overton. Complementarity and nondegeneracy in semidefinite programming. Mathematical Programming, 77(1):111–128, 1997. doi: 10.1007/BF02614432. [2] David Applegate, Mateo Díaz, Oliver Hinder, Haihao Lu, Miles Lubin, Brendan 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. [3] Heinz H. Bauschke, Jonathan M. Borwein, and Wu Li. Strong conical hull intersection property, bounded linear regularity, Jameson’s property (G), and error bounds in convex optimization. Mathematical Programming, 86(1): 135–160, 1999. doi: 10.1007/s101079900086. [4] Silvia Bonettini, Riccardo Zanella, and Luca Zanni. A scaled gradient projection method for constrained image deblurring. Inverse Problems, 25(1):015002, 2009. [5] Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40:120–145, 2011.

24

P REPRINT. U NDER REVIEW.

[6] Antonin Chambolle and Thomas Pock. On the ergodic convergence rates of a first-order primal-dual algorithm. Mathematical Programming, 159(1):253–287, 2016. [7] Chih-Chung Chang and Chih-Jen Lin. LIBSVM: a library for support vector machines. ACM Transactions on Intelligent Systems and Technology (TIST), 2(3):1–27, 2011. [8] Govind M. Chari and Behçet Açıkmeşe. QOCO-GPU: A quadratic objective conic optimizer with GPU acceleration. arXiv preprint arXiv:2603.29197, 2026. [9] Kaihuang Chen, Defeng Sun, Yancheng Yuan, Guojun Zhang, and Xinyuan Zhao. Hpr-qp: A dual halpern peaceman-rachford method for solving large-scale convex composite quadratic programming. arXiv preprint arXiv:2507.02470, 2025. [10] Yuwen Chen, Danny Tse, Parth Nobel, Paul Goulart, and Stephen Boyd. CuClarabel: GPU acceleration for a conic optimization solver. arXiv preprint arXiv:2412.19027, 2024. [11] Chek Beng Chua and Levent Tunçel. Invariance and efficiency of convex representations. Mathematical Programming, 111(1–2):113–140, 2008. doi: 10.1007/s10107-006-0072-6. [12] Yu-Hong Dai and Roger Fletcher. Projected Barzilai-Borwein methods for large-scale box-constrained quadratic programming. Numerische Mathematik, 100(1):21–47, 2005. [13] Timothy A Davis and Yifan Hu. The university of florida sparse matrix collection. ACM Transactions on Mathematical Software (TOMS), 38(1):1–25, 2011. [14] Qi Deng, Qing Feng, Wenzhi Gao, Dongdong Ge, Bo Jiang, Yuntian Jiang, Jingsong Liu, Tianhao Liu, Chenyu Xue, Yinyu Ye, et al. An enhanced alternating direction method of multipliers-based interior point method for linear and conic optimization. INFORMS Journal on Computing, 37(2):338–359, 2025. [15] Lin Ding and Madeleine Udell. A strict complementarity approach to error bound and sensitivity of solution of conic programs. Optimization Letters, 17(7):1551–1574, 2023. [16] Edmund Eisenberg and David Gale. Consensus of subjective probabilities: The pari-mutuel method. The Annals of Mathematical Statistics, 30(1):165–168, 1959. [17] Jacques Faraut and Adam Korányi. Analysis on Symmetric Cones. Oxford Mathematical Monographs. Clarendon Press, Oxford, 1994. [18] Olivier Fercoq. Quadratic error bound of the smoothed gap and the restarted averaged primal-dual hybrid gradient. arXiv preprint arXiv:2206.03041, 2022. [19] Fabio Furini, Emiliano Traversi, Pietro Belotti, Antonio Frangioni, Ambros Gleixner, Nick Gould, Leo Liberti, Andrea Lodi, Ruth Misener, Hans Mittelmann, et al. QPLIB: a library of quadratic programming instances. Mathematical Programming Computation, 11:237–265, 2019. [20] Yuan Gao and Christian Kroer. Infinite-dimensional fisher markets and tractable fair division. Operations Research, 71(2):688–707, 2023. doi: 10.1287/opre.2022.2344. [21] Dongdong Ge, Qi Huangfu, Zizhuo Wang, Jian Wu, and Yinyu Ye. Cardinal optimizer (copt) user guide. arXiv preprint arXiv:2208.14314, 2022. [22] Alan J. Goldman and Albert W. Tucker. Theory of linear programming. In Harold W. Kuhn and Albert W. Tucker, editors, Linear Inequalities and Related Systems, number 38 in Annals of Mathematics Studies, pages 53–97. Princeton University Press, Princeton, NJ, 1956. [23] Paul J Goulart and Yuwen Chen. Clarabel: An interior-point solver for conic programs with quadratic objectives. arXiv preprint arXiv:2405.12762, 2024. [24] Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, Version 12.0, 2024. https://www.gurobi.com.

25

P REPRINT. U NDER REVIEW.

[25] Yicheng Huang, Wanyu Zhang, Hongpei Li, Dongdong Ge, Huikang Liu, and Yinyu Ye. A restarted primal-dual hybrid conjugate gradient method for large-scale quadratic programming. INFORMS Journal on Computing, 2025. [26] Xin Jiang. Local linear convergence of the primal-dual hybrid gradient method for semidefinite programming. arXiv preprint arXiv:2607.08035, 2026. [27] Shucheng Kang, Xin Jiang, and Heng Yang. Local linear convergence of the alternating direction method of multipliers for semidefinite programming under strict complementarity. arXiv preprint arXiv:2503.20142, 2025. [28] Hongpei Li, Yicheng Huang, Huikang Liu, Dongdong Ge, and Yinyu Ye. D-PDLP: Scaling PDLP to distributed multi-GPU systems. arXiv preprint arXiv:2601.07628, 2026. [29] Tianyi Lin, Shiqian Ma, Yinyu Ye, and Shuzhong Zhang. An ADMM-based interior-point method for large-scale linear programming. Optimization Methods and Software, 36(2-3):389–424, 2021. [30] Zhenwei Lin, Zikai Xiong, Dongdong Ge, and Yinyu Ye. A practical GPU-enhanced matrix-free primal-dual method for large-scale conic programs. arXiv preprint arXiv:2505.00311, 2025. [31] Scott B. Lindstrom, Bruno F. Lourenço, and Ting Kei Pong. Error bounds, facial residual functions and applications to the exponential cone. Mathematical Programming, 200(1):229–278, 2023. doi: 10.1007/s10107-022-01883-8. [32] Benqi Liu, Ju Cao, Wotao Yin, and Zaiwen Wen. Restarted reflected halpern acceleration for augmented primaldual methods. arXiv preprint arXiv:2606.16552, 2026. [33] Yulan Liu, Ying Sun, and Shaohua Pan. 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. doi: 10.1007/s11228-018-0494-3. [34] Miguel Sousa Lobo, Lieven Vandenberghe, Stephen Boyd, and Hervé Lebret. Applications of second-order cone programming. Linear Algebra and its Applications, 284(1-3):193–228, 1998. [35] Haihao Lu and Jinwen Yang. cuPDLP. jl: A GPU implementation of restarted primal-dual hybrid gradient for linear programming in julia. arXiv preprint arXiv:2311.12180, 2023. [36] Haihao Lu and Jinwen Yang. A practical and optimal first-order method for large-scale convex quadratic programming. arXiv preprint arXiv:2311.07710, 2023. [37] Haihao Lu and Jinwen Yang. Restarted halpern pdhg for linear programming. arXiv preprint arXiv:2407.16144, 2024. [38] 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 programming by c language. arXiv preprint arXiv:2312.14832, 2023. [39] 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. [40] Harry Markowitz. Portfolio selection*. The Journal of Finance, 7(1):77–91, 1952. doi: https://doi.org/10.1111/j. 1540-6261.1952.tb01525.x. URL https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1540-6261.1952.tb01525.x. [41] István Maros and Csaba Mészáros. A repository of convex quadratic programming problems. Optimization Methods and Software, 11(1-4):671–681, 1999. [42] David Q. Mayne, James B. Rawlings, Christopher V. Rao, and Pierre O. M. Scokaert. Constrained model predictive control: Stability and optimality. Automatica, 36(6):789–814, 2000. [43] Hans Mittelmann. Decision tree for optimization software. http://plato.asu.edu/guide.html, 2021. [44] Hans Mittelmann. Benchmark of second-order cone programming solvers. https://plato.asu.edu/ftp/socp.html, 2024. [45] MOSEK ApS. The MOSEK Optimizer API for C Manual, Version 11.2, 2025. https://docs.mosek.com. 26

P REPRINT. U NDER REVIEW.

[46] NVIDIA Corporation. NVIDIA cuOpt: GPU-Accelerated Decision Optimization, 2025. https://github.com/ NVIDIA/cuopt. [47] Brendan O’Donoghue, Eric Chu, Neal Parikh, and Stephen Boyd. SCS: Splitting conic solver, version 3.2.4. https://github.com/cvxgrp/scs, November 2023. [48] Brendan O’Donoghue, Eric Chu, Neal Parikh, and Stephen Boyd. Conic optimization via operator splitting and homogeneous self-dual embedding. Journal of Optimization Theory and Applications, 169:1042–1068, 2016. [49] Thomas Pock and Antonin Chambolle. Diagonal preconditioning for first order primal-dual algorithms in convex optimization. In 2011 International Conference on Computer Vision, pages 1762–1769. IEEE, 2011. [50] Stephen M. Robinson. Some continuity properties of polyhedral multifunctions. Mathematical Programming Study, 14:206–214, 1981. [51] Daniel Ruiz. A scaling algorithm to equilibrate both rows and columns norms in matrices. Technical report, CM-P00040415, 2001. [52] Bartolomeo Stellato, Goran Banjac, Paul Goulart, Alberto Bemporad, and Stephen Boyd. Osqp: An operator splitting solver for quadratic programs. Mathematical Programming Computation, 12(4):637–672, 2020. [53] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288, 1996.

27

P REPRINT. U NDER REVIEW.

A

Proof in Section 2

A.1

Proof of Theorem 2.3

The exact contraction is obtained in four steps. We first establish a one-step energy inequality. Summing this inequality over one epoch and using convexity gives the ordinary gap bound for the Cesàro average. The ordinary gap estimate is then converted into a smoothed-gap estimate. Finally, Assumption 2.2 converts the smoothed-gap bound into a contraction of the distance to the KKT set. Define the KKT operator "

# ∂F (x) − A∗ y FKKT (x, y) := . Ax − b

(73)

−1 Then Z ⋆ = FKKT (0).

Lemma A.1 (one-step energy inequality). For every comparison point zb = (b x, yb) ∈ K × Y , Q(z k+1 , zb) ≤

1 k 1 1 ∥z − zb∥2P − ∥z k+1 − zb∥2P − ∥z k+1 − z k ∥2P . 2 2 2

(74)

In particular, for every z ⋆ ∈ Z ⋆ , ∥z k+1 − z ⋆ ∥P ≤ ∥z k − z ⋆ ∥P ,

∥z k − z ⋆ ∥ ≤ κP ∥z 0 − z ⋆ ∥.

(75)

Proof. Optimality of the primal proximal step gives τ −1 (xk − xk+1 ) ∈ ∂F (xk+1 ) − A∗ y k . After adding A∗ (y k − y k+1 ), τ −1 (xk − xk+1 ) + A∗ (y k − y k+1 ) ∈ ∂F (xk+1 ) − A∗ y k+1 . The dual step similarly gives A(xk − xk+1 ) + σ −1 (y k − y k+1 ) = A(xk − xk+1 ) + A(2xk+1 − xk ) − b = Axk+1 − b. Together, these two relations show that g k+1 := P (z k − z k+1 ) ∈ FKKT (z k+1 ). Convexity of F implies Q(z k+1 , zb) ≤ ⟨g k+1 , z k+1 − zb⟩. The three-point identity in the P -inner product now gives (74). At a saddle point, Q(z k+1 , z ⋆ ) ≥ 0. Dropping this term and the final nonnegative squared norm proves Fejér monotonicity in (75). Finally, −1/2

∥z k − z ⋆ ∥ ≤ λ−

−1/2

∥z k − z ⋆ ∥P ≤ λ−

∥z 0 − z ⋆ ∥P ≤ κP ∥z 0 − z ⋆ ∥.

Proposition A.2 (exact one-epoch estimates). Let T −1

z̄ T :=

1 X k+1 z . T k=0

⋆

⋆

Then, for every z ∈ Z , ∥z k − z ⋆ ∥ ≤ κP ∥z 0 − z ⋆ ∥,

∥z̄ T − z ⋆ ∥ ≤ κP ∥z 0 − z ⋆ ∥.

28

(76)

P REPRINT. U NDER REVIEW.

If, in addition, ξ > 0 and T ≥ 2λ+ /ξ, then, for every center ż ∈ E × Y , Gξ (z̄ T ; ż) ≤

λ+ 0 ∥z − ż∥2 . T

(77)

Proof. Summing (74) from k = 0 to T − 1 and discarding nonpositive terminal terms gives T −1 X

Q(z k+1 , zb) ≤

k=0

1 0 ∥z − zb∥2P . 2

The function L(·, yb) is convex and L(b x, ·) is affine. Then, for every zb ∈ K × Y , Jensen’s inequality shows that λ+ 0 1 0 ∥z − zb∥2P ≤ ∥z − zb∥2 , (78) 2T 2T where the second inequality follows from P ⪯ λ+ I. The first estimate in (76) follows from Lemma A.1; convexity of the norm and averaging give the second. Q(z̄ T , zb) ≤

For the smoothed-gap estimate, set a = λ+ /(2T ) and bξ = ξ/2. Then  Gξ (z̄ T ; ż) ≤ sup a∥z 0 − zb∥2 − bξ ∥b z − ż∥2 .

(79)

z b∈K×Y

Relaxing x b ∈ K to x b ∈ E gives an upper bound. For a < bξ , completion of the square gives the unrestricted supremum abξ ∥z 0 − ż∥2 . bξ − a Condition T ≥ 2λ+ /ξ gives a ≤ bξ /2 and hence bξ /(bξ − a) ≤ 2. Substitution into (79) proves (77). Proof of Theorem 2.3. We prove by induction that the restart points remain in the retention ball and that (14) holds. Both claims are immediate for n = 0. Suppose they hold through outer iteration D. For every t ≤ D, choose zt⋆ ∈ ΠZ ⋆ (z t,0 ). Proposition A.2 gives ∥z D+1,0 − z 0,0 ∥ ≤

D X

∥z t+1,0 − z t,0 ∥

t=0

≤

D X

∥z t+1,0 − zt⋆ ∥ + ∥zt⋆ − z t,0 ∥

t=0

≤ (1 + κP )

D X

dt

t=0

≤ (1 + κP )

D X

e−t d0 ≤ R0 .

t=0

Thus the new restart point stays in the ball on which quadratic growth is valid. The same estimate with the last term ⋆ omitted shows that the moving center zD also belongs to that ball. Applying Assumption 2.2 and Proposition A.2 gives ⋆ d2D+1 ≤ αξ−1 Gξ (z D+1,0 ; zD )≤

λ+ D,0 λ+ 2 ⋆ 2 ∥z − zD ∥ = d ≤ e−2 d2D ≤ e−2(D+1) d20 , αξ T αξ T D

where the penultimate inequality is precisely the second condition in (13). This closes both parts of the induction. The proof follows a deliberate order: boundedness is established first, so that the local quadratic-growth constant is valid at the new epoch output; only then is quadratic growth used to obtain contraction. This avoids assuming local retention as an unstated premise.

29

P REPRINT. U NDER REVIEW.

A.2

Proof of Theorem 2.5

We first establish the finite-epoch perturbation estimate used to compare an inexact epoch with its exact shadow orbit, and then prove the theorem by combining this estimate with the exact one-epoch contraction. Lemma 2.4 (perturbation of one finite epoch). For fixed T , τ , and σ, there is Cep < ∞, independent of the epoch number, the inner algorithm, and the error values, such that ∥EeT (z) − ET (z)∥ ≤ Cep

T −1 X

δk ,

(80)

k=0

where ET and EeT are the exact and inexact epoch averages starting from the same point. Proof. Let zek = (xke , yek ) be the exact orbit from the same initial point, and set ak := ∥e xk − xke ∥,

bk := ∥e y k − yek ∥,

a0 = b0 = 0.

The definition of δk and nonexpansiveness of the exact proximal map give ak+1 ≤ ak + τ ∥A∥bk + δk .

(81)

The extrapolated-primal difference is bounded by 2ak+1 + ak , and hence bk+1 ≤ bk + σ∥A∥(2ak+1 + ak ).

(82)

With hk = ak + bk , these inequalities imply hk+1 ≤ Lrec hk + cδ δk for constants depending only on τ , σ, and ∥A∥. Expanding this recursion over the fixed number T bounds every hk , PT −1 and therefore the difference of the epoch averages, by a constant times k=0 δk . Proof of Theorem 2.5. For each n, choose zn⋆ ∈ ΠZ ⋆ (z n,0 ) and define the exact shadow average zen+1 := ET (z n,0 ). When zen+1 and zn⋆ lie in the quadratic-growth neighborhood, Assumption 2.2 and Proposition A.2 give dist2 (zen+1 , Z ⋆ ) ≤ αξ−1 Gξ (zen+1 ; zn⋆ ) ≤

λ+ n,0 ∥z − zn⋆ ∥2 = q02 d2n . αξ T

Consequently, the triangle inequality and Lemma 2.4 yield dn+1 ≤ ∥z n+1,0 − zen+1 ∥ + dist(zen+1 , Z ⋆ ) ≤ Cep εn + q0 dn , so, whenever the local estimate is applicable, dn+1 ≤ q0 dn + Cep εn .

(83)

We next prove local retention under (18)–(19). At n = 0, the rate bound is an equality and the displacement bound below holds because its defining sum is empty. Suppose inductively that dt ≤ qηt d0 through index n and that ∥z n,0 − z 0,0 ∥ ≤ βη

n−1 X t=0

30

dt .

P REPRINT. U NDER REVIEW.

The center and the exact shadow average satisfy ∥zn⋆ − z 0,0 ∥ ≤ βη

n−1 X

dt + dn ,

t=0

and ∥zen+1 − z 0,0 ∥ ≤ βη

n−1 X

dt + (1 + κP )dn .

t=0

Here the second estimate uses Proposition A.2. Since βη ≥ 1 + κP , both right-hand sides are at most βη

n X

qηt d0 ≤ Rη .

t=0

Thus (19) permits the use of (83), which gives dn+1 ≤ (q0 + Cep η)dn = qη dn . Furthermore, ∥z n+1,0 − z n,0 ∥ ≤ ∥z n+1,0 − zen+1 ∥ + ∥zen+1 − zn⋆ ∥ + ∥zn⋆ − z n,0 ∥ ≤ Cep εn + (1 + κP )dn ≤ βη dn . Adding this estimate to the preceding displacement bound gives ∥z n+1,0 − z 0,0 ∥ ≤ βη

n X

dt .

t=0

This closes the retention argument and proves (20). A.3

Proof of Proposition 2.6

Proof. For a proximal input v, let u0 ∈ K and define the projected-gradient iteration
uj+1 = ΠK uj − ηPG [Quj + c + τ −1 (uj − v)] . Projection onto K is nonexpansive and ∥I − ηPG Hτ ∥ =

max λ∈[µτ ,Lτ ]

|1 − ηPG λ| = ρPG .

The exact proximal point is the unique fixed point of the PG map, and hence ∥uJ − p(v)∥ ≤ ρJPG ∥u0 − p(v)∥.

(84)

If TvJ (u) denotes J such steps, the same nonexpansiveness argument gives ∥TvJ (u) − TvJ′ (u′ )∥ ≤ ρJPG ∥u − u′ ∥ + ηPG τ −1

1 − ρJPG ∥v − v ′ ∥. 1 − ρPG

(85)

The difference between consecutive PG iterates also supplies a computable stationarity certificate. Let ∆uj := uj+1 −uj and define −1 rj+1 := ∇Hv (uj+1 ) + ηPG (uj − uj+1 ) − ∇Hv (uj ). The projection step implies −1 ηPG (uj − uj+1 ) − ∇Hv (uj ) ∈ NK (uj+1 ),

and therefore rj+1 ∈ ∇Hv (uj+1 ) + NK (uj+1 ),

31

−1 ∥rj+1 ∥ ≤ (Lτ + ηPG )∥∆uj ∥.

(86)

P REPRINT. U NDER REVIEW.

Combining this estimate with (22) yields j+1 ∥uj+1 − p(v)∥ ≤ µ−1 ∥. τ ∥r

It remains to control the warm start uniformly over a local epoch. At every KKT point the warm start is the exact proximal point, so the approximate PG step fixes the KKT point. The affine dual update and extrapolation are Lipschitz, while (85) has a parameter-Lipschitz constant bounded uniformly in J. Induction over the T steps of one epoch therefore gives ∥e z n,k − zn⋆ ∥ ≤ Cin dn , zn⋆ ∈ ΠZ ⋆ (z n,0 ). Nonexpansiveness of the exact proximal map then yields ∥e xn,k − pn,k+1 ∥ ≤ (2 + τ ∥A∥)Cin dn =: Cw dn . sh Applying (84) to each proximal subproblem gives δn,k ≤ Cw ρJPG dn . Summing over the T inner steps proves εn ≤ T Cw ρJPG dn , which is (24). Theorem 2.5 with η = T Cw ρJPG has contraction factor qJ := q0 + Cep T Cw ρJPG . Thus qJ < 1 gives local Q-linear convergence. For 0 < ρPG < 1, the explicit sufficient condition J>

log(Cep T Cw /(1 − q0 )) − log ρPG

(87)

ensures qJ < 1. If ρPG = 0, one PG step is exact.

B

Proofs in Section 3

B.1

Proof of Lemma 3.1

Proof. Take x⋆1 , x⋆2 ∈ X ⋆ and set d = x⋆1 − x⋆2 . Since the feasible set is convex and both points are optimal, the objective is constant on the segment joining them. The second directional derivative of f along this segment is ⟨d, Qd⟩, and hence ⟨d, Qd⟩ = 0. Because Q ⪰ 0, one has Q1/2 d = 0 and therefore Qd = 0, proving that Qx⋆1 = Qx⋆2 . Take two KKT pairs (x⋆1 , y1⋆ ) and (x⋆2 , y2⋆ ). The invariance just established gives s(x⋆2 , y1⋆ ) = Qx⋆2 + c − A∗ y1⋆ = s(x⋆1 , y1⋆ ) ∈ K ∗ . Let d = x⋆2 − x⋆1 . Feasibility gives Ad = 0, while equality of the two objective values and Qd = 0 give 0 = f (x⋆2 ) − f (x⋆1 ) = ⟨Qx⋆1 + c, d⟩ = ⟨s(x⋆1 , y1⋆ ), d⟩. Together with complementarity at x⋆1 , this proves complementarity at x⋆2 . Hence (x⋆2 , y1⋆ ) is also a KKT pair. Crosspairing follows, and the distance identity follows from the Euclidean product norm. The same invariance makes DQ independent of x⋆ . Every optimal slack belongs to the intersection in (26). Conversely, if s = Qx⋆ + c − A∗ y belongs to this intersection, then primal feasibility, dual feasibility, and complementarity hold. Thus (x⋆ , y) is a KKT pair and s is an optimal slack. B.2

Proof of Theorem 3.4

Proof. We first derive the component identities used in both directions. Insert and subtract L(x⋆ , y ⋆ ) in (5). The maximizations over yb and x b separate, and the unused part of each maximization has value zero at the corresponding ⋆ component of z . Hence Gξ (z; z ⋆ ) = Pξ (x; z ⋆ ) + Dξ (y; z ⋆ ). (88) 32

P REPRINT. U NDER REVIEW.

Completing the square in yb gives Pξ (x; z ⋆ ) = f (x) − f (x⋆ ) + ⟨y ⋆ , b − Ax⟩ +

1 ∥Ax − b∥2 . 2ξ

Using Ax⋆ = b, s⋆ = Qx⋆ + c − A∗ y ⋆ , and ⟨x⋆ , s⋆ ⟩ = 0, we obtain Pξ (x; z ⋆ ) = EP (x; z ⋆ ) +

1 ∥Ax − b∥2 . 2ξ

(89)

For the dual component, writing x b = x⋆ + d yields ⋆



Dξ (y; z ) =

sup x⋆ +d∈K

1 −⟨s (y), d⟩ − ∥d∥2Mξ 2 ⋆

 .

The unique maximizer is uy , and therefore 1 Dξ (y; z ⋆ ) = ⟨s⋆ (y), ry ⟩ − ∥ry ∥2Mξ . 2

(90)

Weighted projection optimality gives Mξ ry − s⋆ (y) ∈ NK (uy ). Testing this normal-cone relation with x⋆ ∈ K shows that ⟨s⋆ (y), ry ⟩ ≥ ∥ry ∥2Mξ , and hence Dξ (y; z ⋆ ) ≥

1 ∥ry ∥2Mξ . 2

(91)

Moreover, Dξ (·; z ⋆ ) is convex and differentiable. Danskin’s theorem and the uniqueness of uy give ∇Dξ (y; z ⋆ ) = A(uy − x⋆ ) = −Ary .

(92)

The same projection condition identifies the zero set of the residual: Rξ (y; z ⋆ ) = 0 ⇐⇒ −s⋆ (y) ∈ NK (x⋆ ) ⇐⇒ s⋆ (y) ∈ DQ ∩ Fx⋆ = S ⋆ ⇐⇒ y ∈ Y ⋆ . In particular, (90) gives Dξ (y; z ⋆ ) = 0 for every y ∈ Y ⋆ . Suppose first that the two uniform local error bounds hold. Choose witness neighborhoods VP , Vx , VD , Vy from Definitions 3.2 and 3.3, and choose V contained in VP ∩ VD ∩ (Vx × Vy ), shrinking it if necessary so that both bounds apply to every z ∈ (K × Y ) ∩ V and every z ⋆ ∈ Z ⋆ ∩ V . From (28) and (89),   1 1 Pξ (x; z ⋆ ) ≥ min , dist2 (x, X ⋆ ). aP 2ξbP Similarly, (29) and (91) give Dξ (y; z ⋆ ) ≥

1 dist2 (y, Y ⋆ ). 2κ2D

Thus (88) and (25) prove (30) with  αξ = min

1 1 1 , , aP 2ξbP 2κ2D

 .

Conversely, suppose (30) holds on a neighborhood VQ with constant αξ > 0. Choose product neighborhoods Ux of x̄ and Uy of ȳ such that Ux × Uy ⊆ VQ , and restrict the KKT centers to Z ⋆ ∩ (Ux × Uy ). For x ∈ K ∩ Ux , apply (30) to z = (x, y ⋆ ). By (25), (89), and the definition of Pξ , dist2 (x, X ⋆ ) ≤

1 1 1 Pξ (x; z ⋆ ) = EP (x; z ⋆ ) + ∥Ax − b∥2 . αξ αξ 2ξαξ

This is the uniform local primal error bound with aP = αξ−1 and bP = (2ξαξ )−1 .

33

P REPRINT. U NDER REVIEW.

For y ∈ Y ∩ Uy , apply (30) to z = (x⋆ , y) and let yb ∈ ΠY ⋆ (y). By Lemma 3.1, (x⋆ , yb) is a KKT pair, and the zero-set characterization above gives Dξ (b y ; z ⋆ ) = 0. Convexity and (92) therefore yield −1/2

αξ dist2 (y, Y ⋆ ) ≤ Dξ (y; z ⋆ ) ≤ ⟨∇Dξ (y; z ⋆ ), y − yb⟩ ≤ ∥AMξ

∥ ∥ry ∥Mξ dist(y, Y ⋆ ).

After dividing when the distance is nonzero, we obtain the uniform local dual error bound with any positive constant satisfying −1/2 ∥AMξ ∥ κD ≥ . αξ When the distance is zero the bound is immediate. This proves the reverse implication. The center neighborhood Ux × Uy , together with Ux and Uy for the primal and dual variables, provides the common witness neighborhoods required in Definitions 3.2 and 3.3. B.3

Proof of Proposition 3.6

Proof. The lifted cone K × Qr is again a finite product of the listed classes, and the canonical certificates constructed above exhibit lifted KKT points, so the lifted problem (31) satisfies the standing strong-duality hypotheses of [15]. By the transfer of strict complementarity, P-SC of every center in the stratum becomes lifted P-SC of its canonical certificate, which is the strict-complementarity hypothesis of [15, Corollary 1]. Fix z ⋆ ∈ Z ⋆ ∩ VSC . Remark 4 of [15] assembles the violation-of-complementarity function of the product cone K × Qr as the sum of the blockwise functions of their Lemma 3; each summand is linear or of square-root type in ⟨ζ ⋆ , χ⟩, with coefficients determined by the blockwise spectra of the components of ζ ⋆ (smallest nonzero entries or eigenvalues; the RSOC block enters through p⋆ = (1, u⋆ , −Bx⋆ ), whose norm is at least one). Remark 1 bounds the norm arguments of these functions on the bounded neighborhood N , and [15, Corollary 1] then yields, for every χ ∈ (K × Qr ) ∩ N , dist(χ, XL⋆ ) ≤ c1 (z ⋆ )⟨ζ ⋆ , χ⟩ + c2 (z ⋆ )⟨ζ ⋆ , χ⟩1/2 + c3 (z ⋆ )∥AL χ − eb∥. Since ⟨ζ ⋆ , χ⟩ is bounded on N , the linear term is absorbed into the square-root term, which is (32) at the fixed center. It remains to choose the constants uniformly. The mirror image of the cross-complementarity argument of Corollary 3.9 shows that P-SC at every center fixes the primal complementary faces K ∩ (s⋆ )⊥ along the stratum, and with them the blockwise supports, boundary rays, and ranks of the slacks s⋆ . The blockwise smallest nonzero entries and eigenvalues therefore vary continuously with z ⋆ and remain bounded below after shrinking VSC , while the complementary subspaces entering c3 , the map AL , and the norm bound on N are fixed. Hence supz⋆ ∈Z ⋆ ∩VSC ci (z ⋆ ) < ∞ for i = 1, 2, 3, and the constants in (32) can be chosen uniformly. B.4

Proof of Corollary 3.7

Proof. Lemma 3.1 and Q = B ∗ B imply that Bx⋆ is constant on X ⋆ . Denote this value by w̄ and set ū := 12 ∥w̄∥2 . Tightness of the lifted epigraph gives XL⋆ = X ⋆ × {(ū, 1, w̄)}, (93) and hence dist2 (b χ(x), XL⋆ ) = dist2 (x, X ⋆ ) + ∥q(x) − (ū, 1, w̄)∥2 ≥ dist2 (x, X ⋆ ).

(94)

The map χ b is continuous, so Vx := χ b−1 (N ) is a common neighborhood of x̄. For x ∈ K ∩ Vx , the lifted point is cone feasible and AL χ b(x) − eb = (Ax − b, 0, 0). For every x ∈ K, ⟨p⋆ , q(x)⟩ =

1 1 1 ∥Bx∥2 + ∥Bx⋆ ∥2 − ⟨Bx⋆ , Bx⟩ = ∥x − x⋆ ∥2Q , 2 2 2

(95)

and hence ⟨ζ ⋆ , χ b(x)⟩ = ⟨s⋆ , x⟩ + ⟨p⋆ , q(x)⟩ = EP (x; z ⋆ ). 34

(96)

P REPRINT. U NDER REVIEW.

Combining (96), (94), and (32) gives dist(x, X ⋆ ) ≤ κL

p EP (x; z ⋆ ) + γL ∥Ax − b∥.

Squaring and using (a + b)2 ≤ 2a2 + 2b2 yields dist2 (x, X ⋆ ) ≤ 2κ2L EP (x; z ⋆ ) + 2γL2 ∥Ax − b∥2 . Taking VP := VSC and Vx as above proves the stated constants uniformly over the local stratum. B.5

Proof of Corollary 3.9

Proof. Every optimal slack lies in the affine set DQ , whose relative interior is itself. Thus (35) gives ri DQ ∩ ri Fx⋆ ̸= ∅, which implies bounded, and hence local, linear regularity [3]. For the uniform claim, take two primal centers x⋆1 , x⋆2 and D-SC optimal slacks s⋆i ∈ ri Fx⋆i . The cross-complementarity in Lemma 3.1 gives s⋆1 ∈ Fx⋆2 and s⋆2 ∈ Fx⋆1 . Since each slack is a relative-interior point of its face, these inclusions imply Fx⋆1 ⊆ Fx⋆2 , Fx⋆2 ⊆ Fx⋆1 . Hence all complementary faces on the stratum coincide with a fixed face F. Bounded linear regularity of the fixed pair {DQ , F} supplies common WA and κA , proving the uniform claim. B.6

Proof of Proposition 3.11

Proof. The graph of NK is the product of the blockwise graphs and squared distances add across blocks, so it suffices to produce common neighborhoods and moduli for each block; the product then takes the smallest radius and the largest modulus. For a polyhedral block, gph NKj is a finite union of polyhedra, so the upper-Lipschitz modulus of the polyhedral multifunction NKj is determined by its finitely many affine pieces and is valid, with a common neighborhood, at every reference point [50].
For the remaining blocks, the target set does not move: for a cone, NKj (x⋆j ) = − Kj∗ ∩ (x⋆j )⊥ = −Fx⋆j depends on x⋆j only through its complementary face, and the cross-complementarity argument of Corollary 3.9 shows that blockwise D-SC fixes these faces along the stratum: all local centers share one face Fj per block. Blocks with x̄j interior to Kj , or with s̄j interior to Kj∗ , are trivial: in the first case nearby graph points satisfy vj = 0 ∈ NKj (x⋆j ), and in the second case x⋆j = 0, so vj ∈ NKj (uj ) ⊆ −Kj∗ = NKj (x⋆j ); in both cases the left-hand side of (36) vanishes. Consider a PSD block on the boundary. Constancy of Fj fixes ker x⋆j and hence the rank, so the smallest nonzero ⋆ eigenvalue λ+ min (xj ) is bounded below by some λj > 0 after shrinking the stratum neighborhood. A graph point (uj , vj ) satisfies uj ⪰ 0, vj ⪯ 0, and uj vj = 0. Decompose both matrices in blocks adapted to range(x⋆j ) ⊕ ker(x⋆j ) and let vj11 , vj12 , vj22 denote the blocks of vj . Zeroing vj11 and vj12 produces a matrix of −Fj , so dist(vj , −Fj ) ≤ ∥vj11 ∥ + 2∥vj12 ∥. For ∥uj − x⋆j ∥ ≤ λj /2, the restriction of uj to range(x⋆j ) is invertible with inverse bounded by 2/λj , and the first block row of the identity uj vj = 0 give ∥vj11 ∥ + 2∥vj12 ∥ ≤

C∥vj ∥ ∥uj − x⋆j ∥ λj

for an absolute constant C. Since ∥vj ∥ is bounded on the graph neighborhood of (x̄j , −s̄j ), this is (36) with a modulus ⋆ depending only on λj and the neighborhood. A boundary SOC block is the same computation with λ+ min (xj ) replaced ⋆ by the distance of xj from the cone vertex, which is bounded below because the fixed face pins the boundary ray carrying x⋆j ; a rotated second-order cone reduces to this case through a linear isometry. Taking the smallest blockwise radius and the largest blockwise modulus completes the proof.

35

P REPRINT. U NDER REVIEW.

B.7

Proof of Proposition 3.12

Proof. For the projection point and residual defined above, recall that vy = Mξ ry − s⋆ (y). The maps (z ⋆ , y) 7→ s⋆ (y) and (z ⋆ , y) 7→ (uy , vy ) are jointly continuous and, at y = y ⋆ , equal s⋆ and (x⋆ , −s⋆ ), respectively. After shrinking common neighborhoods VD of z̄ and Vy of ȳ, both assumptions therefore apply to every z ⋆ ∈ Z ⋆ ∩ VD and y ∈ Y ∩ Vy . Weighted projection optimality gives vy ∈ NK (uy ). Since NK (x⋆ ) = −Fx⋆ , dist(s⋆ (y), Fx⋆ ) = dist(−s⋆ (y), NK (x⋆ )) ≤ ∥Mξ ry ∥ + κB ∥ry ∥ ≤

Lξ + κB p ∥ry ∥Mξ . λξ

The vector s⋆ (y) already lies in DQ . Therefore Assumption 3.8 gives Lξ + κB dist(s⋆ (y), S ⋆ ) ≤ κA p ∥ry ∥Mξ . λξ Let sb ∈ ΠS ⋆ (s⋆ (y)). Both slacks belong to DQ , so their difference lies in range A∗ . Define d := (A∗ )† (s⋆ (y) − sb),

yb := y + d.

Then s⋆ (b y ) = sb, and hence yb ∈ Y ⋆ . Consequently, dist(y, Y ⋆ ) ≤ ∥d∥ ≤ ∥(A∗ )† ∥ dist(s⋆ (y), S ⋆ ), which proves (29). B.8

Proof of Corollary 3.13

Proof. The listed blocks are symmetric cones, up to the linear isometry carrying rotated Lorentz cones onto Lorentz cones, so P-SC and D-SC are each equivalent, blockwise, to the displayed condition x⋆ + s⋆ ∈ int K [11, 17]. Hence every local KKT center satisfies P-SC, and every local primal center admits a blockwise D-SC optimal slack: by Lemma 3.1, each x⋆ ∈ X ⋆ near x̄ pairs with ȳ to form a KKT center in the stratum, whose optimal slack is then D-SC. The first property and Proposition 3.6 supply the lifted residual error bound of Assumption 3.5, and Corollary 3.7 converts it into the uniform primal error bound with aP = 2κ2L and bP = 2γL2 . The second property, Corollary 3.9, and Proposition 3.11 verify Assumptions 3.8 and 3.10. Proposition 3.12 then gives the uniform dual error bound, and Theorem 3.4 proves the asserted growth estimate. Substituting the primal constants and the value of κD from Proposition 3.12 gives (37).

C

Cone Projections under the Euclidean Norm and the M -Norm

This section gives the projections under the Euclidean norm and the M -norm used in Section 4.1. Let M = diag(m) ≻ 0. For any closed convex cone K, the membership conditions r∈K

− M r ∈ K∗

and

M give ΠM K (r) = r and ΠK (r) = 0, respectively. In all remaining cases, the projected point lies on the boundary of K. For a diagonal metric, the change of variables w = M 1/2 u gives −1/2 ΠM ΠM 1/2 K (M 1/2 r). K (r) = M

(97)

Consequently, the M -norm projections onto the second-order and exponential cones are precisely the Euclidean projections onto the diagonally rescaled cones analyzed in PDCS [30]. We record their formulas below in the original cone coordinates. For the power cone, we derive the Euclidean and M -norm scalar equations and their bisection procedure explicitly.

36

P REPRINT. U NDER REVIEW.

Table 12: Projection methods for the box and cone blocks under the Euclidean norm and the M -norm. Euclidean projection

M -norm projection

Closed form Closed form Closed form Closed form Scalar Root-finding Scalar Root-finding

Closed form Closed form Scalar Root-finding Scalar Root-finding Scalar Root-finding Scalar Root-finding

Set or cone Box constraints Nonnegative cone Second-order cone Rotated second-order cone Exponential cone Three-dimensional power cone

C.1

Box and Nonnegative Cone

The Euclidean and M -norm projection problems both separate by coordinate and have the same solution: [ΠM Rd (r)]i = max{ri , 0}, +

[ΠM [ℓ,u] (r)]i = min{ui , max{ℓi , ri }}.

(98)

Thus their M -norm and Euclidean projections coincide. C.2

Second-Order Cone

Euclidean projection.

For r = (rv , rt ) ∈ Rd × R, let a := ∥rv ∥2 . The Euclidean projection is   (rv , rt ), a ≤ rt ,    0, a ≤ −rt , ΠKsoc (rv , rt ) =     a + r a + r t t   rv , , otherwise. 2a 2

(99)

M -norm projection. For M = diag(mv,1 , . . . , mv,d , mt ), the nontrivial boundary case is determined by a scalar λ ≥ 0. In the original cone coordinates, the rescaled-cone characterization of PDCS [30] becomes vi (λ) =

mv,i (rv )i , mv,i + λ

t(λ) =

mt rt , mt − λ

∥v(λ)∥22 − t(λ)2 = 0.

(100)

If rt > 0, the admissible root lies in (0, mt ); if rt < 0, it lies in (mt , ∞). When rt = 0, the singular case λ = mt gives vi = mv,i (rv )i /(mv,i + mt ) and t = ∥v∥2 . The scalar equation is solved by bisection; the endpoint signs and the resulting bisection intervals are established in [30]. C.3

Rotated Second-Order Cone

Euclidean projection.

Define the orthogonal map  U (v, s, t) :=

s−t v, √ 2



s+t , √ 2

 .

It maps Krsoc onto a second-order cone, and hence ΠKrsoc (r) = U −1 ΠKsoc (U r).

37

(101)

P REPRINT. U NDER REVIEW.

M -norm projection. For r = (rv , rs , rt ) and M = diag(mv,1 , . . . , mv,d , ms , mt ), a nonzero boundary projection is determined by ξ ≥ 0 through

"

ms −2ξ

mv,i (rv )i vi (ξ) = , mv,i + 2ξ # " # #" ms rs −2ξ s(ξ) = , t(ξ) mt rt mt

(102)

0 = ∥v(ξ)∥22 − 2s(ξ)t(ξ). √ √ √ The 2×2 system is singular at 2ξ = q := ms mt . When its right-hand side is consistent, namely, ms rs + mt rt = 0, the singular solution is obtained explicitly. Define q √ √ mv,i (rv )i δ := ms rs = − mt rt , vi := , χ := δ 2 + 2q∥v∥22 . mv,i + q Then

δ+χ s= √ , 2 ms

t=

−δ + χ √ 2 mt

(103)

satisfies ∥v∥22 = 2st and also covers projections onto either coordinate ray. For nonsingular candidates, the admissible root may lie on either side of q/2. After the membership and zero-projection tests, the implementation brackets a sign change on an interval satisfying s(ξ), t(ξ) ≥ 0 and applies bisection. Restricting the search to the lower branch can miss valid projections. C.4

Exponential Cone

Let r = (r0 , s0 , t0 ). In addition to the membership and dual-cone cases, r0 ≤ 0 and s0 ≤ 0 give the same projection in both metrics: ΠKexp (r) = ΠM Kexp (r) = (r0 , 0, max{t0 , 0}). Euclidean projection. with ηI (ρ) > 0, where

In the remaining case, a point on the smooth boundary is parameterized as ηI (ρ)(ρ, 1, eρ ) ηI (ρ) :=

ρr0 + s0 + eρ t0 , ρ2 + 1 + e2ρ

(104)

hI (ρ) := ρηI (ρ) − r0 + eρ eρ ηI (ρ) − t0 . The Euclidean projection is obtained by finding a root of hI and returning ηI (ρ)(ρ, 1, eρ ). M -norm projection.

Let M = diag(mr , ms , mt ). For the smooth-boundary case, define mr ρr0 + ms s0 + mt eρ t0 , mr ρ2 + ms + mt e2ρ

hM (ρ) := mr ρηM (ρ) − r0 + mt eρ eρ ηM (ρ) − t0 . ηM (ρ) :=

(105)

The projection is obtained from an admissible root of hM satisfying ηM (ρ) > 0 and equals ηM (ρ)(ρ, 1, eρ ). This is the original-coordinate form of the diagonally rescaled exponential-cone projection in PDCS [30], which gives an explicit admissible interval and the corresponding bisection procedure. A safeguarded Newton step can be used within the same bracket. C.5

Three-Dimensional Power Cone

Let α Kpow := {(x, y, z) : x ≥ 0, y ≥ 0, xα y 1−α ≥ |z|},

Its dual cone is

( α (Kpow )∗ =

(p, q, w) : p, q ≥ 0,

38

 p α  α

q 1−α

1−α

α ∈ (0, 1). ) ≥ |w| .

P REPRINT. U NDER REVIEW.

Let r = (rx , ry , rz ). If rz = 0, the Euclidean and M -norm projections coincide:
α ΠKpow (r) = ΠM (r) = max{rx , 0}, max{ry , 0}, 0 . α Kpow Set R := |rz | and suppose below that R > 0. The membership and zero-projection tests are applied separately for each metric. α α Euclidean projection. Suppose that r ∈ / Kpow and −r ∈ / (Kpow )∗ . The nonzero boundary point is parameterized by ρ = |uz | ∈ (0, R) through  p 1 xI (ρ) = rx + rx2 + 4αρ(R − ρ) , 2 q  1 (106) yI (ρ) = ry + ry2 + 4(1 − α)ρ(R − ρ) , 2 gI (ρ) := α log xI (ρ) + (1 − α) log yI (ρ) − log ρ.

If ρ∗I is the zero of gI specified below, then
α ΠKpow (r) = xI (ρ∗I ), yI (ρ∗I ), sign(rz )ρ∗I .

(107)

α α M -norm projection. For M = diag(mx , my , mz ), suppose that r ∈ / Kpow and −M r ∈ / (Kpow )∗ . The same parameterization gives   r 1 mz xM (ρ) = rx + rx2 + 4α ρ(R − ρ) , 2 mx   r mz 1 (108) yM (ρ) = ry + ry2 + 4(1 − α) ρ(R − ρ) , 2 my

gM (ρ) := α log xM (ρ) + (1 − α) log yM (ρ) − log ρ. For stable numerical evaluation of the positive quadratic roots, define  √ s + s2 + 4c    , s ≥ 0, 2 q(s, c) := 2c   √ , s < 0. 2 s + 4c − s

(109)

Thus, for example, xM (ρ) = q(rx , α(mz /mx )ρ(R − ρ)), with analogous expressions for yM , xI , and yI . The second branch avoids cancellation when the input coordinate is negative and ρ is close to an endpoint. If ρ∗M is the zero of gM specified below, then
ΠM (r) = xM (ρ∗M ), yM (ρ∗M ), sign(rz )ρ∗M . α Kpow

(110)

α Proposition C.1 (Bisection for the power-cone projection). Let M = diag(mx , my , mz ) ≻ 0, R = |rz | > 0, r ∈ / Kpow , α ∗ ∗ and −M r ∈ / (Kpow ) . Then gM has a unique zero ρM ∈ (0, R). Initialize ρL = 0 and ρU = R, using the endpoint limits established below. Each bisection step sets ρ = (ρL + ρU )/2 and updates ( ρL ← ρ, gM (ρ) > 0, (111) ρU ← ρ, gM (ρ) < 0.

For a prescribed scalar tolerance ϵbis > 0, the iterations stop when ρU − ρL ≤ ϵbis max{1, R}, and the midpoint is used in (110). The Euclidean procedure is obtained by setting M = I and replacing gM by gI .

39

P REPRINT. U NDER REVIEW.

Proof. Symmetry in the third coordinate allows a nonzero boundary candidate to be written as uz = sign(rz )ρ, where x, y, ρ > 0 and xα y 1−α = ρ. Its KKT conditions with multiplier λ > 0 are ρ mx (x − rx ) = λα , x

ρ my (y − ry ) = λ(1 − α) , y

(112)

mz (ρ − R) + λ = 0, xα y 1−α = ρ. Eliminating λ = mz (R − ρ) and taking the positive roots of the first two quadratic equations yields xM (ρ) and yM (ρ) in (108). Thus every zero of gM in (0, R) satisfies the KKT conditions and recovers the projection. At the upper endpoint, with log 0 := −∞, lim gM (ρ) = log

ρ↑R

1−α (rx )α + (ry )+ R

! < 0,

(113)

√ α where the strict inequality follows from r ∈ / Kpow . As ρ ↓ 0, xM (ρ), and analogously yM (ρ), is of order 1, ρ, or ρ according as rx is positive, zero, or negative. Hence gM (ρ) → +∞ unless both rx < 0 and ry < 0. In that remaining case,    lim gM (ρ) = log  ρ↓0

mz R 

−mx rx α α

−my ry 1−α

 1−α  > 0.

(114)

α The final inequality is exactly the failure of the zero-projection test −M r ∈ (Kpow )∗ . Thus [0, R] is a sign-changing bracket in the limiting sense. Finally, the weighted projection problem has a strictly convex objective, so its KKT point is unique. Consequently, gM has a unique zero in (0, R), and bisection converges to it. Setting M = I proves the Euclidean statement.

D

Large-Scale Instance Statistics

Table 13 reports the dimensions and densities of the large-scale Lasso-derived QP instances used in Section 5.5; m and n denote the number of samples and features of the underlying design matrix A. Table 13: Instance statistics of the large-scale Lasso-derived QPs. Problem SLS rcv1_test avazu-site.tr avazu-app avazu-site kddb2010_test kdda2010_test kddb2010_train kdda2010_train

m

n

Density

1,748,122 677,399 23,567,843 40,428,967 25,832,830 748,401 510,302 19,264,097 8,407,752

62,729 47,236 1,000,000 1,000,000 1,000,000 1,163,024 20,216,830 1,163,024 20,216,830

6.21 × 10−5 1.55 × 10−3 1.50 × 10−5 1.50 × 10−5 1.50 × 10−5 7.74 × 10−6 1.87 × 10−5 7.97 × 10−6 1.80 × 10−6

Table 14 reports the dimensions of the quasilinear Fisher equilibrium instances used in Section 5.6. Here n and m are the numbers of buyers and goods, e = |E| is the number of nonzero buyer-good valuations, N = e + 4n is the number of stored primal coordinates, and M = m + n is the number of affine equations in (72).

40

P REPRINT. U NDER REVIEW.

Table 14: Instance statistics of the quasilinear Fisher equilibrium benchmark.

E

Buyers n

Goods m

Density ρ

Edges e

Variables N

Equations M

1,000 10,000 100,000 1,000,000 10,000,000

400 4,000 4,000 4,000 4,000

0.20 0.02 0.02 0.02 0.01

79,956 798,310 7,999,211 79,996,397 400,022,237

83,956 838,310 8,399,211 83,996,397 440,022,237

1,400 14,000 104,000 1,004,000 10,004,000

Ablation: Inner Accuracy and Preconditioning

We isolate the effect of the monotone inner tolerance (52) and of Jacobi inner preconditioning on the non-diagonal Maros-Mészáros subset, whose instances invoke the inner solver, at tolerance 10−6 with a 1000 second time limit. The PDHCG baseline sets the inner tolerance to 0.05 times the current KKT residual. Table 15: Ablation of adaptive inner tolerance and inner preconditioning on non-diagonal Maros-Mészáros problems. Solver

Solved

Time SGM10 (s)

Iteration SGM10

76 78 79

11.04 6.15 5.14

18,434 6,495 5,690

PDHCG baseline PDHCG-CQP without Jacobi preconditioning PDHCG-CQP with Jacobi preconditioning

The monotone adaptive tolerance sharply reduces runtime and iteration count relative to the PDHCG baseline, and Jacobi inner preconditioning solves one additional hard instance while lowering the runtime SGM10 from 6.15 to 5.14 seconds.

41

Related documents

• Primal-Dual Inexact Newton-MR for Nonconvex Optimization with Equality Constraints#996513
• A proximal augmented Lagrangian method for nonconvex optimization with equality and inequality constraints#1122944
• Elastic ODYN: Differentiable Optimization for Infeasible Control and Learning in Robotics#1141569
• Online Convex Optimization with Sublinear Noisy Probes#271829
• Halpern-Accelerated Majorized ADMM with Optimal Non-Ergodic Convergence under Degenerate Preconditioning#1143409