arXiv CSopen access

A Curvature-Aware Rank-Adaptive Distributed Augmented-Lagrangian Solver for Large-Scale SDPs

A C URVATURE -AWARE R ANK -A DAPTIVE D ISTRIBUTED AUGMENTED -L AGRANGIAN S OLVER FOR L ARGE -S CALE SDP S

Hongpei Li1 , Huikang Liu2 , Dongdong Ge2 , Yinyu Ye2,3 1 Northwestern University 2 Shanghai Jiao Tong University 3 Stanford University [email protected], [email protected]

arXiv:2607.17933v1 [math.OC] 20 Jul 2026

A BSTRACT We present CARDAL (Curvature-Aware Rank-Adaptive Distributed Augmented Lagrangian), a distributed multi-GPU solver for large-scale semidefinite programs (SDPs) based on a rank-adaptive Burer–Monteiro factorization and an augmented Lagrangian method. At fixed ranks, a matrix-free L-BFGS method with negative-curvature corrections targets an approximate Euclidean second-order stationary point of the factored augmented Lagrangian. A reverse multiplier shift turns a negative dual-slack direction into exact negative curvature after rank expansion, and a small joint rank-lift problem selects a batched low-rank correction. A verified slack lower bound provides an a posteriori approximate KKT certificate. Our analysis establishes generic global-optimality guarantees for heterogeneous products of PSD cones at per-block ranks near the Barvinok–Pataki scale, together with a finite-accuracy counterpart under blockwise cost smoothing. For scalable execution, CARDAL distributes constraint rows, factor columns, and PSD blocks over a C ONSTRAINT × R ANK × C ONE device mesh. The primal residual, gradient, Hessian–vector products, and slack matrix–vector products are evaluated using device-local operations and axis-wise collectives. On the Mittelmann benchmark, CARDAL exhibits stronger robustness than existing low-rank GPU approaches under a uniform accuracy standard. Experiments on large-scale SDP relaxations from robotics, electronic structure, and Max-Cut demonstrate the complementary scaling regimes of the three distribution axes, with observed wall-clock speedups of up to 4× on four H100 GPUs.

1

Introduction

Semidefinite programming (SDP) is a fundamental modeling framework in convex optimization [5, 38], with applications in combinatorial optimization [7, 16], control and robotics [19, 20], quantum information [9, 27], and moment–sum-of-squares relaxations of polynomial optimization problems [22, 23]. We consider the primal–dual pair (P) (D)

min ⟨C, X⟩

s.t.

AX = b,

max ⟨b, y⟩

s.t.

S(y) := C − A∗ y ⪰ 0,

X∈Sn

y∈Rm

X ⪰ 0, (1)

Pm where A : Sn → Rm is defined by (AX)i = ⟨Ai , X⟩ and A∗ y = i=1 yi Ai . For many emerging applications, the difficulty of (1) is determined jointly by the matrix order, the number of affine constraints, and the number and heterogeneity of PSD blocks. Interior-point solvers, such as MOSEK [30], COPT [13], SDPT3 [39], SeDuMi [36], and Clarabel [15], remain the preferred choice for high-accuracy solutions at moderate scale. Their scalability is limited, however, by the cost of solving the Newton systems: as the number of constraints and cone dimensions grow, fill-in and factorization costs can dominate both memory consumption and runtime, even when the original SDP data are sparse. Large-scale non-interior-point approaches include bundle methods [11, 26], low-rank nonlinear programming methods such as SDPLR [6, 8], augmented-Lagrangian and semismooth Newton methods such as SDPNAL+ [37, 42, 46], operator splitting such as SCS [31], and conditional-gradient methods such as CGAL and SketchyCGAL [43–45].

P REPRINT. U NDER REVIEW.

They avoid the dense factorizations of a full interior-point step and can therefore address substantially larger problems, at the price of a different accuracy–runtime tradeoff. An important family is based on the augmented Lagrangian method (ALM). At each ALM iteration, a subproblem for the augmented Lagrangian at penalty parameter ρ > 0 is solved, possibly inexactly: ρ 2 X t+1 ≈ arg min Lρ (X, y) = ⟨C, X⟩ − ⟨y, AX − b⟩ + AX − b 2 . (2) 2 X⪰0 Although this subproblem is convex and has only a conic constraint, solving it globally at large scale still requires substantial spectral work to enforce positive semidefiniteness. Conditional-gradient methods avoid a full PSD projection through rank-one updates, but their global rates are sublinear and the rank of the maintained iterate can grow with the requested accuracy. Low-rank factorization offers a different scaling regime. Write meff := rank(A) for the number of independent affine constraints. For a nonempty compact feasible set, the Barvinok–Pataki bound [2, 32] guarantees the existence of an optimal SDP solution of some rank k satisfying k(k + 1)/2 ≤ meff . The Burer–Monteiro (BM) substitution X = F F ⊤ [8], with F ∈ Rn×k and k ≪ n, therefore replaces an O(n2 ) PSD variable by an O(nk) factor. Applied to the convex ALM subproblem (2), it gives F t+1 ≈ arg min Φρ,yt (F ) := Lρ (F F ⊤ , y t ).

(3)

F ∈Rn×k

This reduction raises two fundamental questions: how to choose the working rank k, and how to solve the resulting problem (3) reliably. Although [12, Proposition 2] shows that the ALM subproblem (2) admits an optimal solution satisfying the Barvinok–Pataki rank bound under the dual Slater condition, this regularity condition need not hold in general. More importantly, (3) remains nonconvex because of the quartic penalty term and may possess spurious local minima. Thus, the existence of a low-rank global solution does not by itself provide an algorithmic mechanism for finding one. A standard response is to adapt the rank until the factorized solver provides an approximate global solution of the current convex ALM subproblem. For a fixed multiplier y, define Zρ (F ; y) := ∇X Lρ (X, y)|X=F F ⊤ . Since ∇F Φρ,y (F ) = 2Zρ (F ; y)F , an exact critical point satisfies Zρ (F ; y)F = 0. If, in addition, Zρ (F ; y) ⪰ 0, then F F ⊤ satisfies the KKT conditions of the convex ALM subproblem and is therefore globally optimal. If Zρ (F ; y) has a negative eigenvalue, a corresponding eigenvector provides an escape direction after rank lifting. This leads to the generic adaptive-rank inner loop  F j ∈ crit F 7→ Lρ (F F ⊤ , y) : F ∈ Rn×kj , (4)
Zρ (F j ; y) ⪰̸ 0 =⇒ F j+1 = [ F j αj vj ], vj ≈ vmin Zρ (F j ; y) . The inner loop repeats until an inner global-accuracy contract of the form Lρ (F F ⊤ , y) − inf Lρ (X, y) ≤ ε X⪰0

(5)

is met. The outer iteration can then invoke standard inexact convex-ALM theory [21, 34]. This idea underlies several prominent low-rank SDP solvers, including SDPLR [6], ManiSDP [40], SDPDAL [41], HALLaR [29], and ALORA [12]. The limitation is structural: solving a transient convex ALM subproblem may require a rank substantially larger than that of the final SDP solution. For example, HALLaR’s analysis permits O(1/ε) inner correction rounds in the worst case [29, Theorem 2.5]; since each non-reset Frank–Wolfe correction may add one column, the factor width can grow with the required accuracy. Consequently, both the O(nk) storage and the cost of the inner solve may become prohibitive. Other methods, such as SDPLR [6] and ALORA [12], instead update the multiplier after an approximate fixed-rank solve and a spectral rank-adaptation step, without first certifying (5). This strategy avoids fitting the rank to every transient convex subproblem, but existing guarantees are local: they require the multiplier and factor initialization to be sufficiently close to a strictly complementary primal–dual solution and therefore do not provide a global mechanism for reaching this regime from an arbitrary initialization. These observations lead to the central question of this work: Under what conditions can an ALM–BM method stop rank growth near the Barvinok–Pataki scale while retaining a verifiable optimality guarantee?

2

P REPRINT. U NDER REVIEW.

Another obstacle to applying this framework to substantially larger instances is systems scalability. Recent generalpurpose low-rank GPU solvers, including cuLoRADS [17, 18], cuHALLaR [1], and ALORA [12], show that matrixfree low-rank iterations map well to GPU hardware. Their implementations primarily target a single device, whereas the memory footprint and computational demands of the largest graph and moment relaxations can exceed singledevice capacity. The named-axis methodology of D-PDLP [25] provides a useful starting point, but a low-rank blockdiagonal SDP has a third structural dimension beyond the row–column layout of a linear program. Distributing only the affine constraints is insufficient: work and storage may instead be dominated by the factor columns, a few large PSD blocks, or thousands of small blocks. Contributions. We address the optimization and systems questions with CARDAL (Curvature-Aware RankAdaptive Distributed Augmented Lagrangian), a rank-adaptive distributed multi-GPU ALM–BM solver. Its two main contributions are as follows. • A curvature-aware rank-adaptive ALM–BM framework with product-cone guarantees. At fixed ranks, we target approximate second-order stationarity using L-BFGS, a matrix-free negative-curvature search, and an exact quartic line search. A reverse multiplier shift converts a negative-slack witness into strict negative curvature after zero-padding. A closed-form step handles one direction or a shared batch, while a joint rank-lift problem couples multiple negative-slack directions across PSD blocks. We formulate the geometry, stationarity transfer, and convergence analysis directly for heterogeneous products of PSD cones. Under the stated fixed-rank stationarity and regularity assumptions, accumulation points are second-order critical. For almost every tuple of block costs, their represented matrices are globally optimal once τ (kc ) = kc (kc + 1)/2 exceeds the affine dimension visible to block c. The resulting exact staircase has finite blockwise rank growth. At finite accuracy, slack lower bounds yield deterministic blockwise approximate KKT certificates. We further extend the cost-smoothed AFAC argument of [10] to independently perturbed cost blocks and transfer the result back to the nominal costs. • A composable three-axis multi-GPU decomposition. We distribute constraint rows, factor columns, and PSD blocks over a C ONSTRAINT ×R ANK ×C ONE device mesh and derive exact distributed formulas for residuals, gradients, Hessian–vector products, and slack matrix–vector products. These primitives support distributed spectral searches, heterogeneous block ranks, cost-aware cone assignment, and batching of small blocks. The three axes address complementary regimes and can be combined for mixed-structure instances. The distributed operator evaluations coincide with the product-cone operators used in the convergence, genericlandscape, smoothed-analysis, and a posteriori certification results. Organization. Section 2 develops the product-cone BM geometry, AL-to-BM stationarity transfer, blockwise landscape, and finite-accuracy theory. Section 3 presents the rank-adaptive algorithm, fixed-rank convergence, exact rankgrowth analysis, and finite-output guarantees. Section 4 gives the multi-GPU decomposition and distributed operator identities. Section 5 collects block operators, scaling, and implementation safeguards. Section 6 reports the numerical experiments, and Section 7 concludes.

2

Theoretical Foundations for Product Cones

CARDAL is designed for semidefinite programs with heterogeneous PSD blocks. We therefore formulate the theory directly on a product cone and treat a single PSD cone as the special case q = 1. This section answers a pointwise question: which stationarity and curvature conditions at a factor imply exact or approximate optimality of the original SDP? Section 3 then gives conditions under which CARDAL produces such factors. 2.1

Product-cone BM formulation and geometry

For a positive integer q, write [q] := {1, . . . , q}. Consider the primal–dual pair min

{Xc ⪰0}qc=1

s.t.

q X ⟨Cc , Xc ⟩,

max

y∈Rm

c=1 q X

b⊤ y, (6)

Ac (Xc ) = b,

s.t.

c=1

3

Sc (y) := Cc − A∗c (y) ⪰ 0

∀c ∈ [q],

P REPRINT. U NDER REVIEW.

where Ac : Snc → Rm and redundant affine rows are allowed. We write A⊕ (H1 , . . . , Hq ) :=

q X

Ac (Hc ),

meff := rank(A⊕ ).

c=1

The formulation includes a conventional single-matrix SDP when q = 1 and a block-diagonal encoding in which the Xc are the diagonal blocks. For a rank profile k = (k1 , . . . , kq ), apply the Burer–Monteiro substitution Xc = Fc Fc⊤ ,

Fc ∈ Rnc ×kc ,

F = (F1 , . . . , Fq ).

Product factor spaces use the direct-sum inner product and norm, ⟨F, U ⟩⊕ :=

q X

∥F ∥2⊕ :=

⟨Fc , Uc ⟩,

c=1

q X

∥Fc ∥2F .

c=1

The corresponding fixed-rank BM problem is min g(F ) := F

q X

⟨Cc , Fc Fc⊤ ⟩

c(F ) :=

s.t.

c=1

q X

Ac (Fc Fc⊤ ) − b = 0,

(BMk )

c=1

with feasible set Mk := {F : c(F ) = 0}. For U = (Uc )c , Dc(F )[U ] =

q X

(7)

Ac (Fc Uc⊤ + Uc Fc⊤ ).

(8)

c=1

Define

q BF ν := A∗c (ν)Fc c=1 ,

G(F ) := BF∗ BF .

(9)

Then Dc(F ) = 2BF∗ . Assumption 1 (Product-cone BM smoothness). For every rank profile k under consideration with Mk ̸= ∅, there is an open neighborhood Uk ⊃ Mk on which Dc(F ) has constant rank. For q = 1, Assumption 1 is the constant-rank alternative of the standard BM smoothness assumption [4, Assumption 1.1(b)]. The constant-rank theorem makes Mk an embedded submanifold, with TF Mk = ker Dc(F ),

NF Mk = range BF .

(10)

PF (Z) := Z − BF G(F )† BF∗ Z = Projker Dc(F ) (Z)

(11)

Moreover, G(F )† and the orthogonal tangent projector

vary smoothly on Uk . Linear independence of the factor-space constraint gradients is not required; if Dc(F ) has full row rank, the compatible multiplier introduced next is unique. For a rank staircase, nonemptiness is required at its initial profile. Zero-padding preserves feasibility at every later profile, but smoothness must hold at each visited profile; neither compactness nor boundedness of the algorithmic iterates is included in Assumption 1. Let CF := (Cc Fc )qc=1 and define the canonical multiplier and block slacks µB (F ) := G(F )† BF∗ CF ,

SB,c (F ) := Cc − A∗c (µB (F )).

(12)

The pseudoinverse fixes a canonical representative when compatible multipliers are not unique. At a feasible factor, the Riemannian gradient and Hessian quadratic form are
q grad g(F ) = 2SB,c (F )Fc c=1 , (13) ⟨U, Hess g(F )[U ]⟩⊕ = 2

q X

⟨Uc , SB,c (F )Uc ⟩,

c=1

4

U ∈ TF Mk .

(14)

P REPRINT. U NDER REVIEW.

Definition 1 (Product-cone BM second-order criticality). A factor F ∈ Mk is a second-order critical point of (BMk ) if q X SB,c (F )Fc = 0 (c ∈ [q]), ⟨Uc , SB,c (F )Uc ⟩ ≥ 0 for every U ∈ ker Dc(F ). (15) c=1

For y ∈ Rm and ρ > 0, define the product factored augmented Lagrangian ρ Φρ,y (F ) := g(F ) − ⟨y, c(F )⟩ + ∥c(F )∥22 , min Φρ,y (F ). F 2

(16)

Although the augmented Lagrangian is quadratic in the matrix tuple (Xc )c , it is quartic in F . Introduce the shifted multiplier yb(F ; y, ρ) := y − ρc(F ). (17) Direct differentiation in the product factor space gives
q ∇Φρ,y (F ) = 2Sc (b y )Fc c=1 , (18)  2  ∗ y )Uc + 2ρAc (Dc(F )[U ]) Fc , (19) ∇ Φρ,y (F )[U ] c = 2Sc (b ⟨U, ∇2 Φρ,y (F )[U ]⟩⊕ = 2

q X ⟨Uc , Sc (b y )Uc ⟩ + ρ∥Dc(F )[U ]∥22 .

(20)

c=1

Proposition 1 (Product-cone SOSP transfer). Suppose Assumption 1 holds at profile k. For any y ∈ Rm and ρ > 0, every feasible full-space second-order stationary point of Φρ,y satisfies (15). For q = 1, Proposition 1 is exactly the single-cone AL-to-BM transfer. The implication is one-way: a manifold SOSP controls curvature only on ker Dc(F ), whereas the augmented Hessian acts on the full factor space. When the constraint gradients are dependent, a compatible ALM multiplier need not equal µB (F ) as a vector, but the two slacks have the same action on F and the same quadratic form on TF Mk . This representation-invariant compatibility is the key step in the proof. The positive penalty term in (20) also shows why exact feasibility matters. It vanishes on tangent directions at a feasible factor, but can mask negative curvature at an infeasible iterate. Finite-accuracy statements therefore retain the feasibility and curvature residuals explicitly. 2.2

Generic landscape and blockwise rank conditions

Let

k(k + 1) . (21) 2 For one cone, the classical BM landscape result states that every second-order critical point is globally optimal for every cost when k ≥ n, and for almost every cost when k < n and rc := rank(Ac ) = dim Im(A∗c ),

τ (k) > rank(A)

τ (k) :=

(22)

under the standard smoothness assumption [4, Proposition 3.1, Corollary 3.2, and Lemma 3.3]. Theorem 1 extends the single-cone conclusions in [4, Theorem 1.4] to heterogeneous products of PSD cones with a separate rank condition for each block. Theorem 1 (Generic product-cone landscape). Suppose Assumption 1 holds at profile k and, for every c ∈ [q], kc ≥ nc or τ (kc ) > rc . (23) Q nc Then, for almost every C = (C1 , . . . , Cq ) ∈ c S , every product-cone BM second-order critical point is globally optimal for (6). If kc ≥ nc for every block, the conclusion holds for every product cost. Condition (23) uses the dimension of the affine directions visible to each block. The coarser choice τ (kc ) > meff P is always sufficient, whereas a condition only on c τ (kc ) is not: columns assigned to one block cannot create null directions with which to test the slack of another block. Theorem 1 reduces to the classical low-rank branch (22) when q = 1.

5

P REPRINT. U NDER REVIEW.

2.3

Finite-accuracy certificates and cost smoothing

At finite accuracy, the feasibility, stationarity, and curvature residuals may be nonzero. We use the route approximate product AL SOSP =⇒ product AFAC pair =⇒ approximate product-SDP optimality. The first implication is deterministic. The second is a blockwise extension of the cost-smoothed single-cone AFAC guarantee in [10, Theorem 6], formalized for product cones in Theorem 2. Definition 2 (Approximate AL second-order stationarity). A product factor F is an (η, ζ)-approximate second-order stationary point of (16) if
2 ∥∇Φρ,y (F )∥⊕ ≤ η, λEuc (24) min ∇ Φρ,y (F ) ≥ −ζ. The minimum eigenvalue is taken over the full direct-sum factor space. Definition 3 (Approximate product-SDP optimality). A pair ((Xc )c , λ) is (ε0 , ε1 , ε2 )-approximately optimal for (6) if X Ac (Xc ) − b ≤ ε0 , Xc ⪰ 0 (c ∈ [q]), c

2

(25)

!1/2 X

∥Sc (λ)Xc ∥2F

≤ ε1 ,

Sc (λ) ⪰ −ε2 Inc

(c ∈ [q]).

c

Proposition 2 (Objective interpretation of approximate KKT conditions). Let ((Xc )c , λ) satisfy (25), and let (Xc⋆ )c be any feasible tuple. Then s X X X X ⋆ ⟨Cc , Xc ⟩ − ⟨Cc , Xc ⟩ ≤ ε0 ∥λ∥2 + nc ε1 + ε2 tr(Xc⋆ ). (26) c

c

c

c

Thus bounded multipliers and a bounded comparison solution turn the approximate KKT statements below into explicit one-sided objective bounds. Whenever Ac is applied to a nonsymmetric matrix, it denotes the extension induced by the symmetric matrices defining that operator. Hence X Ac (Uc Fc⊤ ). Dc(F )[U ] = 2 c

Definition 4 (Product approximate feasible approximate criticality). For ϵ = (ϵ0 , ϵ1 , ϵ2 ) and γ > 0, a pair (F, λ) is a product (ϵ, γ)-AFAC pair if !1/2 X 2 ∥c(F )∥2 ≤ ϵ0 , ∥Sc (λ)Fc ∥F ≤ ϵ1 , (27) c

X ⟨Uc , Sc (λ)Uc ⟩ ≥ −ϵ2

X

∥U ∥⊕ = 1,

whenever

c

Ac (Uc Fc⊤ )

c

≤ γ.

(28)

2

Proposition 3 (AL residuals imply product AFAC). Suppose F is an (η, ζ)-approximate product AL SOSP and set λ = yb(F ; y, ρ). Then, for every γ > 0, (F, λ) is a product (ϵ, γ)-AFAC pair with ϵ0 = ∥c(F )∥2 ,

ϵ1 =

η , 2

ϵ2 =

ζ + 2ργ 2 . 2

(29)

For the smoothed statements, fix a linearly independent representation of the shared affine constraints. All residuals, multiplier bounds, operator norms, and the dimensions rc below refer to this representation. Removing affine redundancies leaves the primal feasible set unchanged but changes multiplier coordinates. Let ∥Ac ∥F →2 denote the operator norm from the matrix Frobenius norm to the Euclidean norm. Theorem 2 (Cost-smoothed product-cone finite-accuracy guarantee). Fix ϵ ∈ R3+ , γ, RF , Rλ > 0, and a rank profile k before sampling the block costs. Independently draw Cc ∼ Unif{D ∈ Snc : ∥D − C̄c ∥F ≤ σc }. Let J := {c : kc ≤ nc ,

ac := ∥Ac ∥F →2 > 0},

6

δc :=

ϵ1 ac , γ

κc := Rλ ac .

P REPRINT. U NDER REVIEW.

Suppose, for every c ∈ J , τ (kc ) > rc , 0 < δc ≤ κc , P Then, with probability at least 1 − c∈J pc , where  pc := 4e

3κc δc

rc 

4n3c δc σc

δc <

σc . 4n3c

(30)

τ (kc ) ,

(31)

every product (ϵ, γ)-AFAC pair satisfying ∥F ∥⊕ ≤ RF and ∥λ∥2 ≤ Rλ obeys X

Xc := Fc Fc⊤ ⪰ 0,

≤ ϵ0 ,

Ac (Xc ) − b

c

2

(32)

!1/2 X

∥Sc (λ)Xc ∥2F

≤ RF ϵ1 ,

Sc (λ) ⪰ −ϵ2 Inc

(c ∈ [q]).

c

Blocks with kc > nc or ac = 0 satisfy the final inequality deterministically and do not contribute a failure probability. The familiar near-Barvinok–Pataki scale follows blockwise. If, for fixed χ, tc > 0,  1/χ  1+1/χ 1 σc τ (kc ) ≥ (1 + χ)rc + χtc , δc ≤ , 3κc 4n3c then  pc ≤ 4e

σc 12κc n3c

σc < 1, 12κc n3c

(33)

 tc .

p For fixed χ and tc , the condition is kc ∼ 2(1 + χ)rc , rather than a full-width factor. Taking q = 1 recovers the single-cone cost-smoothed AFAC guarantee. Remark 1 (Scope of the smoothed guarantee). Theorem 2 applies when the sampled block costs are used throughout the run. Its event is uniform over all bounded AFAC pairs at the prescribed rank profile, which is essential because the output depends on the sampled costs. It does not describe a run performed only at the nominal costs; Section 3.3 transfers the sampled-cost output back to those costs. We use the AFAC-to-optimality argument of [10], but do not claim its polynomial iteration bound for the implemented ALM–L-BFGS–NC procedure. A data-dependent profile or tolerance envelope requires a finite union bound or a separate covering argument.

3

CARDAL Algorithm and Convergence Guarantees

We now address the algorithmic question: under what conditions does CARDAL produce the factors characterized in Section 2, and when does blockwise rank growth stop? The outer method applies ALM directly to the nonconvex product BM problem (BMk ). With the rank profile fixed, its inner routine targets an approximate full-space SOSP of (16); the rank-adaptive wrapper then tests each dual-slack block and enlarges only the blocks with detected negative curvature. 3.1

Rank-adaptive CARDAL

We organize CARDAL into three nested routines. The rank-adaptive wrapper in Algorithm 1 is the only routine that changes the rank profile k = (k1 , . . . , kq ). At each profile, it calls the fixed-rank ALM in Algorithm 2, whose subproblems are solved by the L-BFGS–NC method in Algorithm 3. Because all blocks enter the shared residual c(F ), they use one multiplier and one penalty and are optimized jointly in the direct-sum factor space. Throughout Algorithms 2 and 3, F denotes the tuple (Fc )qc=1 , rather than a single common-size matrix. Its entries may be assembled into a single optimization vector, but the objective, operators, gradients, and search directions retain their blockwise definitions. Only Algorithm 1 changes the widths kc of these blocks. Here fixed-rank means that every entry of k remains unchanged during the ALM call; the block ranks need not be equal. For fixed (y, ρ, k), the inner solver targets ∥∇Φρ,y (F )∥⊕ ≤ ω,

2 λEuc min ∇ Φρ,y (F ) ≥ −ζ.

7

(34)

P REPRINT. U NDER REVIEW.

Rank-adaptive wrapper. Algorithm 1 first runs the fixed-rank ALM and then tests every returned slack block. For block c, the algorithm retains at most rinc,c directions below −τdual . If no block violates the test, the wrapper terminates; otherwise, one joint rank-lift problem couples all retained directions through the shared affine residual. Algorithm 1 CARDAL rank-adaptive wrapper 1: Input: product factor F = (Fc )c , profile k, multiplier ybase , penalty ρ > 0, slack tolerance τdual , and increment caps (rinc,c )c 2: loop 3: (F, y, ρ) ← F IXED R ANK ALM(F, ybase , ρ) 4: for c = 1, . . . , q do 5: Ec ← N EGATIVE E IGENPAIRS(Sc (y), τdual , rinc,c ) 6: end for 7: if Ec = ∅ for every c then 8: return (F, y) 9: else 10: ybase ← y + ρc(F ) 11: (Wc )c ← J OINT R ANK L IFT(y, ρ, (Ec , Ac )c ) 12: for all c with Wc ̸= ∅ do 13: Fc ← [ Fc Wc ] and kc ← kc + cols(Wc ) 14: end for 15: end if 16: end loop The following single-cone lemma isolates the reverse-shift mechanism and its one-dimensional quartic line search. The product update then applies the same identity simultaneously to the selected blocks. To ensure that these rank-lift constructions are well posed, we assume that the primal feasible set is nonempty and compact. This compactness condition also appears in Assumption 2 and rules out nonzero PSD recession directions in the nullspace of the affine map. Lemma 1 (Reverse-shift rank lift and exact line search). Let F ∈ Rn×k be produced by a subproblem with base multiplier ybase and penalty ρ, and set yout = ybase − ρc(F ). ⊤ For any unit vector v, set qv := v S(yout )v and suppose qv < 0. Define F̄ = [ F 0 ],

Uv = [ 0 v ],

av = A(vv ⊤ ).

Then Dc(F̄ )[Uv ] = 0, Moreover, the enlarged-rank objective satisfies

Uv , ∇2 Φρ,ybase (F̄ )[Uv ] = 2qv < 0.

ρ Φρ,ybase (F̄ + tUv ) = Φρ,ybase (F̄ ) + qv t2 + ∥av ∥22 t4 . 2 If av ̸= 0, its exact nonnegative minimizer is s −qv t⋆ = . ρ∥av ∥22

(35) (36)

(37)

Proof. The added column is orthogonal to the zero-padded factor in the sense that F̄ Uv⊤ + Uv F̄ ⊤ = 0. Hence Dc(F̄ )[Uv ] = 0. Moreover, c(F̄ ) = c(F ), so yb(F̄ ; ybase , ρ) = ybase − ρc(F ) = yout . This proves (35) using (20). For any t ∈ R, c(F̄ + tUv ) = c(F ) + t2 av ,

g(F̄ + tUv ) = g(F ) + t2 v ⊤ Cv.

Expanding the augmented Lagrangian and using yout = ybase − ρc(F ) gives (36). Differentiating its even quartic polynomial gives (37). Compactness ensures av ̸= 0 for nonzero v: otherwise X + svv ⊤ would be feasible for every s ≥ 0 whenever X is feasible. Hence (37) is well defined and eliminates the rank-growth scale hyperparameter. The ybase assignment in Algorithm 1 recovers the base multiplier for which this line search is exact.

8

P REPRINT. U NDER REVIEW.

Joint rank lift. Let J be the set of violated blocks. For each c ∈ J , collect the retained orthonormal directions in Vc = [ vc,1 · · · vc,ℓc ] and define Hc := Vc⊤ Sc (yout )Vc ,

⊤ qc,j := vc,j Sc (yout )vc,j < 0,

⊤ ac,j := Ac (vc,j vc,j ).

1/2

For Zc ∈ Sℓ+c , set Fc (Zc ) = [ Fc Vc Zc ] on selected blocks and leave the remaining factors unchanged. Lemma 1 gives X Φρ,ybase (F (Z)) = Φρ,ybase (F ) + ⟨Hc , Zc ⟩ c∈J

ρ X Ac (Vc Zc Vc⊤ ) + 2 c∈J

(38)

2

. 2

Thus all selected cones and directions are coupled through the convex quadratic SDP   2 X  X ρ (Zc⋆ )c∈J ∈ arg min ⟨Hc , Zc ⟩ + . Ac (Vc Zc Vc⊤ ) Zc ⪰0, c∈J   2 c∈J

c∈J

(39)

2

The PSD variables mix directions only within their own cone, while the shared residual couples the cones. The same recession argument makes the objective coercive on the selected PSD subspaces. The correction Wc = Vc (Zc⋆ )1/2 is compressed to the positive eigenspace of Zc⋆ before columns are allocated. Restricting each matrix to Zc = Diag(zc,1 , . . . , zc,ℓc ) gives the joint nonnegative QP   2   ℓc ℓ c  X X ρ XX qc,j zc,j + min . zc,j ac,j z≥0   2 c∈J j=1  c∈J j=1

(40)

2

2

Restricting further to Zc = t Iℓc gives the shared closed-form amplitude v P Pc u u − c∈J ℓj=1 qc,j t⋆ = u t P 2. Pℓc ρ c∈J j=1 ac,j

(41)

2

For one selected direction in one cone, this reduces to (37). Fixed-rank ALM. With k fixed, the method alternates between an approximate subproblem solve and the multiplier update y + = y − ρc(F + ). The penalty is increased only when the primal residual fails to contract:  ρ, ∥c(F + )∥2 ≤ θ∥c(F )∥2 , + ρ = 0 < θ < 1, γρ > 1. (42) γρ ρ, otherwise, Algorithm 2 Fixed-rank augmented-Lagrangian method q 1: Input: product factor F = (Fc )c=1 , profile k, multiplier y, and initial penalty ρ > 0 2: loop 3: δ ← ∥c(F )∥2 and choose inner tolerances (ω, ζ) 4: F ← L-BFGS–NC(F, y, ρ, ω, ζ) 5: r ← c(F ) and y ← y − ρr 6: if ∥r∥2 ≤ ϵ0 then 7: return (F, y, ρ) 8: end if 9: if ∥r∥2 > θδ then 10: ρ ← γρ ρ 11: end if 12: end loop The fixed-rank routine neither computes a slack eigenpair nor changes k. The convergence analysis assumes that the current profile can represent a feasible point. All rank updates are determined by the outer dual-slack test in Algorithm 1.

9

P REPRINT. U NDER REVIEW.

L-BFGS–NC subproblem solver. For fixed (y, ρ, k), Algorithm 3 first reduces the gradient using safeguarded LBFGS steps. Once the gradient tolerance is met, a matrix-free spectral test is applied to the true Hessian. If negative curvature is detected, the method takes a curvature step and resets the L-BFGS memory; otherwise it returns when the second-order tolerance is met. To limit the cost of Hessian–vector products, the practical implementation invokes the negative-curvature search only intermittently after sufficient first-order progress, while the analysis uses the secondorder termination condition above. Algorithm 3 L-BFGS–NC for one augmented-Lagrangian subproblem q 1: Input: product factor F = (Fc )c=1 , profile k, multiplier y, penalty ρ, tolerances (ω, ζ), and maximum step τmax 2: loop 3: while ∥∇Φρ,y (F )∥⊕ > ω do 4: Compute a safeguarded L-BFGS direction D = (Dc )qc=1 5: τ⋆ ← arg min0≤τ ≤τmax Φρ,y (F + τ D) 6: F ← F + τ⋆ D and update the L-BFGS memory 7: end while 2 8: if λEuc min (∇ Φρ,y (F )) ≥ −ζ then 9: return F 10: end if 11: Let V = (Vc )qc=1 be the returned unit negative-curvature direction 12: τ⋆ ← arg min|τ |≤τmax Φρ,y (F + τ V ) 13: F ← F + τ⋆ V and clear the L-BFGS memory 14: end loop Along either search direction, the augmented-Lagrangian objective is a quartic polynomial in the step length, so the exact line search requires only the interval endpoints and the admissible real roots of its cubic derivative. The Hessian test in Algorithm 3 targets second-order stationarity at the current profile, whereas the blockwise slack tests in Algorithm 1 determine which ranks should increase. 3.2

Fixed-rank ALM convergence

Consider a tail on which the rank profile k is fixed. For the asymptotic analysis, suppress the finite residual stopping test. Given (Fℓ , yℓ , ρℓ ), let Fℓ+1 be the accepted output of the inner solve with base multiplier yℓ , and set yℓ+1 = yℓ − ρℓ c(Fℓ+1 ),

δℓ := ∥c(Fℓ )∥2 .

(43)

Assumption 2 (Fixed-rank ALM conditions). Fix a rank profile k and suppose that: (A1) Mk is nonempty, the SDP feasible set is compact, and the accepted factors and relevant inner level sets lie in a compact set; (A2) Assumption 1 holds at k; (A3) every Fℓ+1 satisfies (34) for (yℓ , ρℓ ), with ωℓ → 0 and ζℓ → 0; (A4) ρ0 > 0, the multipliers {yℓ } are bounded, and the uncapped residual-shrinkage rule (42) is used. Lemma 2 (Vanishing primal residual). Under Assumption 2, δℓ −→ 0.

(44)

Proof. If ρℓ is eventually constant, successful residual tests eventually give δℓ+1 ≤ θδℓ . Otherwise ρℓ → ∞ along the penalty-increase indices, and boundedness of the multipliers in (43) gives 2 supj ∥yj ∥2 −→ 0 ρℓ along those indices; every intervening successful test contracts the residual by at least θ. Hence the full sequence converges to zero. δℓ+1 ≤

Remark 2 (Finite penalties in the local regime). The diverging-penalty branch above is a proof device, not a prediction. Standard local ALM theory permits a sufficiently large finite penalty to remain constant when an exact local minimizer branch satisfies the usual regularity and second-order sufficient conditions [3, Proposition 2.7]. That result explains practical penalty stabilization but does not by itself prove local convergence of the inexact L-BFGS–NC inner solver.

10

P REPRINT. U NDER REVIEW.

Theorem 3 (Fixed-rank product-cone ALM convergence). Under Assumption 2, every accumulation point of the accepted factor sequence is feasible and satisfies (15). If, in addition, the blockwise rank conditions (23) hold and the product cost lies outside the exceptional set of Theorem 1, then every accumulation point is globally optimal for (6),
⊤ g(Fℓ ) → p⋆ , dist (Fℓ,c Fℓ,c )c , X ⋆ → 0, (45) where X ⋆ is the product-SDP solution set. If kc ≥ nc for every block, the conclusion holds for every product cost. The single-cone result is recovered when q = 1. Because each factor is invariant under Fc 7→ Fc Qc , the theorem is stated in terms of accumulation points and represented matrices; convergence of an entire factor sequence would require additional quotient-manifold or finite-length assumptions. 3.3

Blockwise rank growth and finite outputs

We separate three conclusions. An idealized exact staircase explains finite blockwise rank growth, an a posteriori slack test certifies any finite output meeting the stated residual bounds, and cost smoothing gives a low-rank finite-accuracy statement. None of these conclusions bounds the numerical work at a given rank. Define the first certifying rank of block c by kcgen := min{k ≥ kc,0 : k ≥ nc or τ (k) > rc }.

(46)

An exact product staircase solves each visited profile to a feasible full-space AL SOSP, tests every slack block with zero tolerance, and enlarges precisely the blocks with a negative slack eigenvalue. Corollary 1 (Finite exact blockwise rank growth). Suppose Assumption 1 holds at every visited profile up to kgen , and exclude the finite union of the exceptional cost sets from Theorem 1 over those profiles. Then no block in an exact product staircase is enlarged beyond kcgen . If block c receives rinc,c columns per update, with the final update truncated at kcgen , the total number of accepted block updates is at most  q  gen X kc − kc,0 . (47) rinc,c c=1 The terminal tuple satisfies the primal–dual KKT conditions and is globally optimal. This exact statement concerns staircase stages; it does not assert finite time to reach an exact SOSP. More generally, any infinite run with a bounded integer rank profile has a fixed-rank tail after its final accepted update, so Theorem 3 applies if Assumption 2 holds there. Reaching a rank cap alone is not a certificate: every block at its rank cap must still pass the slack test. Proposition 4 (Finite-output product approximate KKT certificate). Let Xc = Fc Fc⊤ and let yout = ybase − ρc(F ) be the multiplier returned after a fixed-rank ALM subproblem. Suppose ∥c(F )∥2 ≤ δfeas ,

∥∇Φρ,ybase (F )∥⊕ ≤ ω,

Sc (yout ) ⪰ −τdual Inc (c ∈ [q]),

∥F ∥⊕ ≤ RF .

(48)

Then ((Xc )c , yout ) is  δfeas ,

 RF ω , τdual -approximately optimal 2

(49)

in the sense of Definition 3. This deterministic certificate does not invoke a generic landscape theorem: the blockwise slack tests directly supply approximate dual feasibility. Together with Proposition 2, it also yields a one-sided objective bound whenever the multiplier and a comparison solution are bounded. For the smoothed output statement, choose deterministic envelopes before sampling: δ̄, ω̄, ζ̄, ρ̄, γ, RF , Ry > 0 and set ω̄ ζ̄ , ϵ2 := + 2ρ̄γ 2 . (50) 2 2 Corollary 2 (Cost-smoothed product-cone ALM output). Suppose the independently sampled block costs of Theorem 2 are used throughout the run and, at a rank profile chosen before sampling, the terminal output satisfies
2 ∥c(F )∥2 ≤ δ̄, ∥∇Φρ,y (F )∥⊕ ≤ ω̄, λEuc ρ ≤ ρ̄, (51) min ∇ Φρ,y (F ) ≥ −ζ̄, ϵ0 := δ̄,

ϵ1 :=

11

P REPRINT. U NDER REVIEW.

and, with yb = y − ρc(F ),

∥F ∥⊕ ≤ RF , ∥b y ∥2 ≤ Ry . (52) Under Pthe conditions of Theorem 2, with Rλ = Ry , the sampled-cost output satisfies (32) with probability at least 1 − c∈J pc . Let σmax := maxc σc . For the nominal costs (C̄c )c , the same output is approximately optimal with tolerances
ϵ0 , RF ϵ1 + σmax RF2 , ϵ2 + σmax . (53) If τdual ≥ ϵ2 , no block covered by the theorem can produce a sampled-slack direction below −τdual on the same event. The additive σmax terms are the accuracy floor incurred when transferring the sampled-cost result to the nominal problem. For a terminal profile chosen from a finite predeclared set, apply a union bound to Theorem 2 over the eligible profiles. No additional union over iterates is needed because the event for each prescribed rank profile is uniform over bounded AFAC pairs.

4

Distributed Multi-GPU Parallelization

The difficulty of scaling a semidefinite program is not captured by its ambient matrix dimension alone. Large application instances exhibit at least three distinct forms of structural complexity. Electronic-structure and orbital-relaxation models can contain very large families of affine constraints, making the applications of A and A∗ the dominant cost. Some difficult Max-Cut relaxations require a wide working factor before a rank-adaptive method passes its terminal dual-slack test; factor storage, slack–factor products, and Hessian–vector products then become the bottleneck. Finally, sparse sum-of-squares and moment relaxations naturally produce many PSD cones, often with highly nonuniform dimensions and sparsity patterns. These regimes are not mutually exclusive: a single instance may have many constraints, wide factors, and many heterogeneous cones. This structural diversity motivates a three-axis distribution rather than a one-dimensional data split. We partition the affine constraints, the factor columns, and the PSD cones over a C ONSTRAINT × R ANK × C ONE device mesh. The resulting layout is an exact decomposition of the factorized SDP, not an approximation. Moreover, all global quantities required by the ALM–BM method can be evaluated through two complementary collective patterns: a forward reduction over the Rank and Cone axes, and an adjoint reduction over the Constraint axis. 4.1

Structural regimes and a three-axis decomposition

Consider the block-diagonal form of the primal SDP, minq

{Xc ⪰0}c=1

s.t.

q X ⟨Cc , Xc ⟩, c=1

(54)

q X ⟨Ac,i , Xc ⟩ = bi ,

i = 1, . . . , m,

c=1

where Ac,i is the restriction of the ith constraint matrix to cone c. For clarity, we omit the optional nonnegative component; it can be represented by 1 × 1 PSD cones and distributed in the same way. Let [m] =

PR G

Ip ,

[q] =

p=1

PC G

Λℓ

(55)

ℓ=1

be partitions of the constraints and cones. For each cone, partition the columns of its factor as [kc ] =

P K G

Kc,s ,

  Fc = Fc[1] · · · Fc[PK ] ,

Fc[s] := Fc [:, Kc,s ].

(56)

s=1

Empty Kc,s are allowed when kc < PK . Column partitioning is particularly convenient because it introduces no cross terms: PK X Xc = Fc Fc⊤ = Fc[s] Fc[s]⊤ . (57) s=1

The subscripts R, K, and C denote the Constraint (row), Rank, and Cone mesh axes, respectively; R is therefore not the factor rank.

12

P REPRINT. U NDER REVIEW.

For i ∈ Ip , define the constraint–cone block operator X   ⟨Ac,i , Xc ⟩. Apℓ (XΛℓ ) i :=

(58)

c∈Λℓ

Substituting (57) into (54) gives the fully partitioned BM problem min [s]

{Fc }c∈[q], s∈[PK ]

s.t.

PK X D PC X E X Cc , Fc[s] Fc[s]⊤ , ℓ=1 s=1 c∈Λℓ PC X PK X

(59) 



Apℓ {Fc[s] Fc[s]⊤ }c∈Λℓ = bIp ,

p = 1, . . . , PR .

ℓ=1 s=1

Equation (59) exposes the three independent sums that drive the distributed design. The Constraint axis partitions the output rows of A; the Rank axis partitions the sum forming each Fc Fc⊤ ; and the Cone axis partitions the block sum in each constraint, which may couple multiple cones. The Ndev = PR PK PC devices are indexed by (p, s, ℓ). Device (p, s, ℓ) stores the operator shard {Ac,i : i ∈ Ip , c ∈ [s] Λℓ } and the factor slabs {Fc : c ∈ Λℓ }. Operator shards are replicated along the Rank axis, factors along the Constraint axis, and each objective block Cc along both the Constraint and Rank axes within its cone group. Row vectors such as the multiplier and primal residual are replicated along the Rank and Cone axes. We write ARR , ARK , and ARC for summation by AllReduce along the corresponding mesh axis; the result remains replicated along the reduced axis. 4.2

Distributed evaluation of the ALM–BM operators

The primal residual, gradient, Hessian action used by the curvature test, and matrix–vector products used by the dual-slack test are all compositions of a local forward product and a local adjoint–factor product. Making these two primitives explicit both shortens the derivation and clarifies where communication is unavoidable. For factor-shaped families U and V with the same distribution as F , device (p, s, ℓ) forms the local forward partial E XD   Πpsℓ (U, V ) i := Ac,i , Uc[s] Vc[s]⊤ , i ∈ Ip . (60) c∈Λℓ

For a row vector zp ∈ R|Ip | and a factor-shaped family U , the corresponding local adjoint–factor partial is X Ψc,ps (zp ; U ) := (zp )i Ac,i Uc[s] , c ∈ Λℓ .

(61)

i∈Ip

The first object is row-shaped and is reduced over the Rank and Cone axes; the second is factor-shaped and is reduced over the Constraint axis. In particular, ARR (Ψc,ps (zp ; U )) = (A∗ z)c Uc[s] .

(62)

This ordering communicates only a thin nc × |Kc,s | product rather than an assembled nc × nc slack matrix. Primal residual.

Using (57), the residual slice owned by constraint group p is rp = ARK (ARC (Πpsℓ (F, F ))) − bIp .

(63)

There is no reduction over the Constraint axis: every component of rp belongs to exactly one constraint group. The shifted multiplier ybp = yp − ρrp is therefore formed locally. Augmented-Lagrangian gradient. For c ∈ Λℓ , the local factor slab of the gradient is  [s] ∇Φρ,y (F ) c = 2Cc Fc[s] − 2ARR (Ψc,ps (b yp ; F )) = 2Sc (b y )Fc[s] .

(64)

Thus the residual and gradient realize the same forward/adjoint duality as the distributed Ax/A⊤ y pattern in D-PDLP, with the additional Cone reduction induced by block-diagonal semidefinite structure [25].

13

P REPRINT. U NDER REVIEW.

Curvature test. on

Let V be a spectral-search vector distributed like F . The nonlocal part of the Hessian action depends

ξ(V ) := A(F V ⊤ + V F ⊤ ). Its constraint-group slice is obtained by the same forward reduction as the primal residual: ξp (V ) = ARK (ARC (Πpsℓ (F, V ) + Πpsℓ (V, F ))) .

(65) (66)

Define the combined local thin partial Θc,ps (V ) := −Ψc,ps (b yp ; V ) + ρ Ψc,ps (ξp (V ); F ).

(67)

The two adjoint contributions can therefore be accumulated before communication, giving the local Hessian slab with one Constraint-axis reduction:  [s] c ∈ Λℓ . (68) Hρ,y (F )[V ] c = 2Cc Vc[s] + 2ARR (Θc,ps (V )) , Consequently, every Hessian–vector product requires one Rank–Cone forward reduction and one Constraint-axis thinproduct reduction; no dense Hessian is formed. Dual-slack test. When the fixed-rank ALM returns, its final shifted multiplier yb is the returned multiplier yout in Algorithm 2. The outer rank-adaptation test then examines the smallest eigenvalue of each relevant slack block Sc (yout ) = Cc − (A∗ yout )c . Its matrix–vector products are evaluated without assembling the slack:   X Sc (yout )v = Cc v − ARR  (yout,p )i Ac,i v  . (69) i∈Ip

Each cone group tests its local cones, assigning each cone to a single Rank coordinate s = πK (c) ∈ [PK ] to avoid redundant work. Different cone tests may be assigned to different Rank coordinates. A negative Rayleigh quotient provides a valid witness for increasing that block’s rank. The resulting extremal-eigenvalue estimates and rank-growth flags require only scalar gathers or reductions among the owning devices. This formulation is consistent with the factor-shaped communication in (64); it does not assume that a full dual-slack matrix has already been assembled. For each retained subspace Vc , the operator evaluations Ac (Vc Zc Vc⊤ ) required by the joint rank-lift problem reuse the forward-reduction pattern in (63). The coefficients of this small problem are assembled by scalar reductions. Each resulting correction is compressed to its nonzero factor columns before those columns are assigned to Rank shards. Exactness of the distributed operator actions. The partitions in (55) and (56) are disjoint and exhaustive. Linearity of A and A∗ therefore shows that (63), (64), (68), and (69) coincide exactly with their single-device counterparts in exact arithmetic. In floating-point arithmetic they are algebraically identical up to the reassociation induced by collective reductions. Thus distribution does not change the operator being evaluated. The remaining inner-solver operations require only small-payload scalar collectives. Factor-space quantities such as gradient norms, L-BFGS curvature-pair products, and spectral-search normalization reduce over the Rank and Cone axes, whereas residual-space quantities reduce over the Constraint axis. The coefficients of the quartic line search combine both types and are reduced over the axes on which their nonreplicated partial sums reside. Each reduction is performed on one replica whenever its input is replicated, so that the same contribution is not counted PR , PK , or PC times. Their payload is small compared with the thin operator reductions. The product-cone decomposition in this section is algebraic: it changes neither the objective nor the derivatives. Consequently, its assembled operators are exactly those used by the theory and algorithm in Sections 2 and 3. 4.3

Workload-aware partitioning

Equal numbers of constraints, factor columns, or cones rarely imply equal work: operator nonzeros may be clustered and cone dimensions may vary widely. Following the nonzero-aware data-layout principle of D-PDLP [25], CARDAL first reorders constraint rows to improve column locality and then partitions contiguous row ranges by cumulative operator nonzeros. PSD cones remain indivisible. Large cones and batches of compatible small cones are assigned to Cone coordinates by a longest-processing-time heuristic with weights that reflect each cone’s dimension and rank cap. At initialization, factor columns are divided as evenly as possible over the Rank coordinates. These operations change only data ordering and ownership; the distributed operators remain the exact rearrangements established above. The mesh dimensions (PR , PK , PC ) are runtime configuration parameters satisfying PR PK PC = Ndev . Increasing one dimension reduces the corresponding local constraint, factor, or cone work, but also introduces collectives along that axis. The appropriate shape therefore depends on the dominant structure of the instance and on whether its local work is large enough to amortize collective latency.

14

P REPRINT. U NDER REVIEW.

5

Practical Enhancements

This section describes the numerical transformations, implementation safeguards, and structured operator organization used by CARDAL. They do not modify the ALM–BM control flow of Section 3, the exact distributed identities of Section 4, or the stationarity conditions in the analysis. 5.1

Block-structured cone operations

Let F = blkdiag(F1 , . . . , Fq ) and write Ac,i for the restriction of the ith constraint matrix to cone c. The inner solver accesses the block SDP through three algebraic primitives: q X
A(F F ⊤ ) i = ⟨Ac,i , Fc Fc⊤ ⟩, (70) c=1 ∗

(A z)c Uc =

X

zi Ac,i Uc ,

(71)

i

Sc (z)Uc = Cc Uc − (A∗ z)c Uc . (72) Here Uc may be a factor slab, a Hessian–vector direction, or a slack Lanczos vector. Thus (72) covers the gradient, curvature, and dual-slack actions derived in Section 4 without requiring a dense slack matrix. For each cone, define the union sparsity pattern [

Sc := supp(Cc ) ∪

supp(Ac,i ).

(73)

i:Ac,i ̸=0

Only entries of Fc Fc⊤ indexed by Sc are required in (70); likewise, the adjoint and slack actions are supported on Sc . The fixed data and union patterns are stored sparsely, and the adjoint is contracted with Uc before the distributed reduction, as in Section 4. Consequently, each pass scales with the active nonzeros times the relevant factor or vector width, rather than with the dense cone sizes n2c . Neither the primal matrix F F ⊤ , the Hessian, nor a dense global slack is formed or communicated. 5.2

Batching small cones

Sparse SOS, moment, and chordally decomposed relaxations often contain a few large cones together with many small cones. Although the operators in (70)–(72) are already sparse, applying them cone by cone leaves too little work in each small block to amortize fixed operator overhead. The direct-sum structure permits an exact aggregation. For a batch B of cones, define


AB,i := blkdiag Ac,i c∈B , SB (z) := blkdiag Sc (z) c∈B , FB := blkdiag Fc c∈B . Then ⟨AB,i , FB FB⊤ ⟩ =

X

(74)

⟨Ac,i , Fc Fc⊤ ⟩, (75)

c∈B

SB (z)FB = blkdiag Sc (z)Fc c∈B . Here blkdiag denotes a direct sum; it does not imply explicit storage of the off-diagonal zero blocks. Batching therefore changes only the granularity at which the same block operators are evaluated; it does not change the residual, gradient, or spectral tests. Cones with similar dimensions, ranks, and sparsity densities are grouped together, while large cones remain separate. This distinction is particularly important when q is large but most nc are small. 5.3

Scaling and preconditioning

The constraint rows and cone blocks can have widely different magnitudes, which degrades both the ALM penalty update and the conditioning of the factorized subproblems. CARDAL combines ℓ∞ Ruiz equilibration with Pock– Chambolle scaling [33], subject to the requirement that variable scaling preserve the PSD cones. A generic entrywise scaling of a vectorized symmetric matrix is not, in general, a congruence and need not preserve positive semidefiniteness. We therefore restrict the variable scaling within each cone to a positive diagonal congruence. Let ri > 0 be a constraint-row scale, let Dc ≻ 0 be diagonal, and let ηb , ηC > 0 be global scaling factors for the right-hand side and objective. The primal variables are related by ec Dc−1 , Xc = ηb−1 Dc−1 X

15

−1/2

Fc = ηb

Dc−1 Fec .

(76)

P REPRINT. U NDER REVIEW.

The corresponding scaled data are ec,i = r−1 D−1 Ac,i D−1 , A i

c

c

ebi = ηb r−1 bi , i

ec = ηC Dc−1 Cc Dc−1 . C

(77)

If ye is the multiplier of the scaled problem, then yei yi = , Sec (e y ) = ηC Dc−1 Sc (y)Dc−1 . (78) ri ηC Because Dc is positive diagonal, these maps preserve PSD order, rank, primal and dual feasibility, and complementarity. They also give X X ec , X ec ⟩ = ηb ηC ⟨C ⟨Cc , Xc ⟩. c

c

The Ruiz and Pock–Chambolle iterations determine ri and the diagonal entries of Dc ; the default implementation uses this per-coordinate diagonal mode rather than one scalar per cone. Scalar nonnegative variables use the analogous one-dimensional map. 5.4

Solver defaults

P Let nlp denote the number of optional scalar nonnegative components not included in c nc ; algebraically, they may be treated as 1 × 1 PSD cones. The settings in Table 1 define the experimental protocol unless an experiment states otherwise. We use a comprehensive KKT stopping test following the DIMACS error framework of the Mittelmann benchmark [28]. Let X = blkdiag(Fc Fc⊤ )c , let p = ⟨C, X⟩ and d = b⊤ y, and let S(y) = blkdiag(Sc (y))c ; optional bc (y) denote the smallest slack eigenvalue returned scalar nonnegative components are included as 1 × 1 blocks. Let λ b b by the block spectral routine and set λ⊕ (y) := minc λc (y), including the scalar blocks. In the original minimization convention, the solver-reported relative measures are b⊕ (y)]+ [−λ |p − d| ∥A(X) − b∥2 (79) , , dres := , gap := pres := 1 + ∥b∥1 1 + ∥C∥1 1 + |p| + |d| where [t]+ := max{t, 0} and ∥ · ∥1 is the entrywise data norm. Thus dres tracks the most negative detected slack direction rather than the Frobenius distance to the product cone. The gap compares the primal minimization value with the candidate dual value, which is a lower bound when S(y) ⪰ 0. All quantities in (79), all objectives, and every reported slack tolerance are evaluated after applying (76)–(78) back to the original, unscaled problem; internal scaled residuals are not reported. The spectral routines are implemented matrix-free. The inner Lanczos search estimates the smallest eigenpair of the fixed-rank ALM Hessian, while the outer search estimates a capped set of the most negative Ritz pairs of each slack block for rank adaptation. A shifted power iteration supplies the terminal slack estimate used in the reported dual residual. None of these operations forms a dense Hessian or slack matrix. Table 1: Default numerical settings and implementation policies used by CARDAL. Quantity Termination Initial BM rank Per-cone rank cap Inner quasi-Newton solve Curvature search Rank-growth search Rank-lift update Terminal slack estimate ALM penalty Inner-iteration cap Scaling

Default Relative primal residual, dual residual, and objective gap at most 10−4 Raw value ⌈2 ln m⌉, subject to the cone dimension and the per-cone cap below √  ( 8mc + 1 − 1)/2 , where mc = |{i : Ac,i ̸= 0}|; used as a heuristic safeguard L-BFGS history size 5; exact quartic line search on [0, 1] for L-BFGS steps and on [−1, 1] for negative-curvature steps At most 15 matrix-free Hessian–Lanczos steps At most 30 Lanczos steps per cone; at most the prescribed rinc,c mostnegative directions are retained Small joint rank-lift problem over the retained slack subspaces; resulting corrections are compressed before rank sharding Shifted power pP iteration per cone, used in the reported5 dual residual ρ0 = 2/ c nc + nlp ; multiplier 3.33; cap 5 × 10 30 000 Ten Ruiz sweeps, Pock–Chambolle scaling with α = 1, per-coordinate positive diagonal PSD congruences, and bound/objective rescaling enabled

16

P REPRINT. U NDER REVIEW.

6

Numerical Experiments

The evaluation addresses two complementary questions. We first assess the single-GPU solver on the heterogeneous Mittelmann sparse-SDP benchmark, using both solver-reported termination and the benchmark’s uniform DIMACS accuracy test. We then examine how the distributed implementation responds to three structural regimes that motivate the device mesh of Section 4: vehicle-landing moment relaxations contain many PSD blocks and stress the Cone axis; electronic-structure SOS relaxations contain millions of moment-matching constraints and stress the Constraint axis; and large Max-Cut relaxations produce wide Burer–Monteiro factors and stress the Rank axis. These families emphasize, but do not statistically isolate, the corresponding axes: problem size, sparsity, rank trajectory, and ALM iteration count can change together. The experiments therefore characterize observed end-to-end scaling and identify when computation along each axis amortizes the associated collectives. Experimental protocol. All CARDAL runs use FP64 arithmetic and a wall-clock limit of 3600 s. Unless stated otherwise, the stopping criteria and numerical parameters are those of Section 5.4; MOSEK and cuLoRADS use their default parameter configurations. In the Mittelmann study, baseline results are taken from the public leaderboard [28], and CARDAL is run on an NVIDIA H100 GPU matching the hardware reported for cuLoRADS [18]. All CARDAL experiments use one node with two 48-core Intel Xeon Platinum 8468 CPUs and four NVIDIA H100 SXM GPUs, each with 80 GB of memory; the distributed studies use at most four of these GPUs, interconnected by NVLink/NVSwitch. In the (PR , PK , PC ) convention of Section 4, the Cone-, Constraint-, and Rank-axis studies use meshes (1, 1, p), (p, 1, 1), and (1, p, 1), respectively. For p GPUs, we report the observed wall-clock speedup Sp :=

T1 , Tp

(80)

computed from the underlying timings before display rounding. The symbol f means that the target specified for that table was not reached within 3600 s; for a native-termination table, this means that the corresponding solver did not report successful termination under its own criterion. The symbol does not distinguish a timeout from other unsuccessful termination modes. For runs using the default stopping rule, native success is determined by the 10−4 tolerances on the three solver-reported quantities in (79). Table-specific throughput checkpoints introduced below are reported only as scaling measurements. Columns for comparison solvers provide contextual reference points and are not used to compute Sp . 6.1

Mittelmann sparse-SDP benchmark

We begin with the standard sparse-SDP benchmark of Mittelmann [28], which contains 75 instances spanning combinatorial optimization, control, and relaxation problems. The leaderboard reports three complementary statistics. Solved counts the instances on which a solver terminates according to its own stopping rule; because these rules differ, this column is not a uniform accuracy comparison. DIMACS counts the instances satisfying the benchmark’s common tolerance on the six scaled measures of primal infeasibility, dual infeasibility, and duality gap. SGM10 is the shifted geometric mean of solution times with a 10 s shift and is the leaderboard’s aggregate efficiency metric; lower values are better. Table 2: Mittelmann sparse-SDP benchmark (75 instances; 3600 s limit). Solved uses each solver’s native stopping rule, DIMACS applies the benchmark’s uniform accuracy test, and lower values of SGM10 are better. CARDAL is shaded. Solver

solved/75

DIMACS/75

SGM10

GPU solvers CARDAL cuLoRADS

65 68

65 56

113.03 96.00

CPU solvers COPT MOSEK CSDP SDPT3 SeDuMi SDPA

75 66 65 67 52 60

70 66 60 54 51 42

37.69 122.66 158.27 156.08 634.05 185.85

17

P REPRINT. U NDER REVIEW.

On SGM10, CARDAL ranks third among the listed solvers, behind COPT and cuLoRADS and ahead of MOSEK, SDPT3, CSDP, SDPA, and SeDuMi. Its DIMACS count is 65/75, the highest among the noncommercial solvers in the table and nine more than cuLoRADS. The Solved and DIMACS columns answer different questions and should therefore be read together. For example, cuLoRADS reports 68 instances solved by its internal criterion, of which 56 satisfy the uniform DIMACS test. For CARDAL, both counts are 65. This agreement is consistent with its comprehensive native KKT test, but the table does not make solver-specific stopping rules directly comparable. We therefore use DIMACS for cross-solver accuracy comparisons and SGM10 for aggregate timing. 6.2

Cone-axis scaling on sparse moment relaxations

We study Cone-axis parallelism using semidefinite relaxations of the vehicle-landing optimal-control problem of [19], which we abbreviate as VL. Its sparse Lasserre relaxation organizes consecutive trajectory states and controls into overlapping cliques. Each stage contributes one 190 × 190 moment block and approximately ten 19 × 19 localizing blocks, while consensus constraints couple adjacent stages. The resulting family has 11N − 1 cones at horizon N , so increasing N exposes more indivisible blocks to the Cone-axis assignment. We generate VLN for N ∈ {5, 10, 20, 30, 50, 100, 150}; the original formulation uses N = 50. Table 3 reports the time to each solver’s native successful termination, so its baseline columns are contextual rather than same-criterion comparisons. Table 3: Cone-axis scaling on vehicle-landing moment relaxations. Entries are times in seconds to each solver’s native successful termination; speedups use the one-GPU time; f denotes failure within 3600 s; the fastest CARDAL time in each row is bold.

N cones 5 10 20 30 50 100 150

n

m

1 GPU

2 GPUs

4 GPUs

time

time speedup

time speedup

54 1 881 89 942 9.16 8.99 109 3 781 177 906 72.07 63.93 219 7 581 353 823 92.17 76.65 329 11 381 529 706 145.30 106.78 549 18 981 872 903 211.96 156.55 1099 37 981 1 743 628 537.82 342.99 1649 56 981 2 614 431 1797.68 968.48

1.02× 17.04 1.13× 45.87 1.20× 69.11 1.36× 78.76 1.35× 131.96 1.57× 224.14 1.86× 686.03

0.54× 1.57× 1.33× 1.84× 1.61× 2.40× 2.62×

MOSEK cuLoRADS time

time

214.11 432.87 800.81 1294.26 2400.62 f f

f f f f f f f

Table 3 exhibits the expected crossover. VL5 is too small to benefit from four GPUs, whereas four GPUs are the fastest CARDAL configuration from VL10 onward. Relative to one GPU, the four-GPU speedup ranges from 1.33× on VL20 to 2.62× on VL150; on the two largest instances, the runtime falls from 537.82 to 224.14 s and from 1797.68 to 686.03 s, respectively. Thus Cone-axis parallelism becomes more useful as the number of assignable blocks grows, while the smallest instance remains communication-bound. MOSEK reports native successful termination for instances up to VL50; neither baseline reaches its own native target on VL100 or VL150. Because the baseline stopping criteria differ from CARDAL’s native criterion, these columns are contextual timing references rather than same-criterion speedups. 6.3

Constraint-axis scaling on electronic-structure SOS relaxations

We use the spin-free level-2 fermionic SOS relaxation of Low et al. [27, Appendix G]. The exact optimum of this convex relaxation lower bounds the ground-state energy. Each N -orbital instance has three PSD cones of orders  2  2N , N, N , n = 2N (N + 1), m = Θ(N 4 ), (81) so the computational burden is dominated by the moment-matching constraints rather than by the number of cones. The benchmark suite from [27] contains Fe2 S2 , Fe4 S4 , FeMoCo54 , FeMoCo76 , and the P450CPD1X active space. FeMoCo76 is the largest chemistry instance in that study, where cuLoRADS was used for the SDP computation. The default solver tolerance is not an energy-accuracy statement. Moreover, these instances exhibit a long, highpenalty tail during which conditioning effects in the outer ALM obscure the scaling of the distributed operators. We therefore measure time to the relaxed throughput checkpoint pres < 10−3 ,

gap < 10−2 .

18

(82)

P REPRINT. U NDER REVIEW. Table 4: Constraint-axis scaling on chemistry SOS-SDPs. Entries are times in seconds to (82); speedups use the oneGPU time; f denotes failure within 3600 s; the fastest CARDAL time in each row is bold.

n

1 GPU

2 GPUs

4 GPUs

time

time speedup

time speedup

Problem

N

m

Fe2 S2 Fe4 S4 FeMoCo54 P450CPD1X FeMoCo76

20 840 160 K 2.63 2.48 36 2 664 1.68 M 11.87 5.62 54 5 940 8.5 M 34.72 17.15 58 6 844 11.3 M 121.45 74.82 76 11 704 33 M 318.91 131.24

1.06× 2.37 2.11× 4.93 2.02× 8.84 1.62× 44.68 2.43× 78.36

1.11× 2.41× 3.93× 2.72× 4.07×

MOSEK cuLoRADS time

time

f f f f f

f f f f f

This checkpoint is deliberately looser than the default stopping test and is used only to compare Constraint-axis wallclock scaling, not to assess final energy accuracy. Under the same checkpoint and time limit, neither MOSEK nor cuLoRADS reaches the target in Table 4. The smallest instance shows little benefit from additional devices. For the four larger problems, the observed four-GPU speedup ranges from 2.41× to 4.07×, with the largest values on the FeMoCo instances. P450CPD1X has a comparable largest block order but a smaller speedup; this difference may reflect both workload balance and a different ALM or rank trajectory, so it should not be attributed to row imbalance alone. Overall, the larger values of m provide more local operator work with which to amortize the Constraint-axis thin-product reduction. 6.4

Rank-axis scaling on Max-Cut relaxations

Finally, we study the Goemans–Williamson Max-Cut relaxation [14] on graphs from the Stanford SNAP collection [24]: 1 max ⟨L(W ), X⟩ s.t. diag(X) = 1, X ⪰ 0, (83) X∈SN 4 where L(W ) is the graph Laplacian. This is a single-cone SDP with n = m = N ; the main dimension available for distribution is the width of the adaptive factor F ∈ RN ×k . The Rank axis partitions the columns of F , while the corresponding constraint images are combined by the forward reduction of Section 4. To match the minimization convention used throughout the paper, the solver is given 1 C = − L(W ), 4

A(X) = diag(X),

b = 1.

(84)

Hence the reported primal value ⟨C, X⟩ is the negative of the relaxation’s maximization value, and the residuals and primal–dual gap use the definitions in (79) without a sign convention change. We construct nested real-world graph instances using h-cores [35]: thresholds h ∈ {8, 6, 5, 4, 3, 2} are applied to the SNAP com-DBLP graph, and the 2-core of com-Amazon supplies a similarly sized instance from a different domain. Table 5 reports the surviving vertex and edge counts. For this scaling study, time is measured to the same throughput checkpoint as in (82), pres < 10−3 , gap < 10−2 . Table 5: Rank-axis scaling on Max-Cut SDP relaxations of SNAP graphs. Entries are times in seconds to (82); speedups use unrounded timings; f denotes failure within 3600 s; the fastest CARDAL time in each row is bold. 1 GPU Instance DBLP 8-core DBLP 6-core DBLP 5-core DBLP 4-core DBLP 3-core DBLP 2-core Amazon 2-core

N

|E|

33 875 285 867 68 070 459 741 98 942 584 570 143 129 730 284 204 031 884 337 271 646 1 004 432 305 892 896 901

2 GPUs

4 GPUs

MOSEK cuLoRADS

time time speedup time speedup

time

time

2.89× 1.79× 1.68× 1.56× 2.51× 1.87× 2.54×

f f f f f f f

1.32 1.97 2.76 3.73 4.69 7.79 16.63

1.22 1.30 4.53 2.10 3.74 7.11 7.36

0.42 0.73 2.70 1.35 1.49 3.80 2.90

19

0.84 1.12 1.64 0.69 1.01 2.32 2.14

1.45× 1.16× 2.77× 3.06× 3.69× 3.07× 3.44×

P REPRINT. U NDER REVIEW.

The two smallest DBLP instances are latency dominated: two GPUs are faster than four, and the one-GPU runtimes are only 1–2 s. From DBLP 5-core onward, four GPUs are consistently fastest, with observed speedups between 2.77× and 3.69×. These ratios should be read as wall-clock observations rather than precise estimates of parallel efficiency, especially for the shortest runs and for values slightly above the device count. Collective reassociation and changes in the adaptive-rank or ALM trajectory can also make the multi-GPU execution non-work-preserving. The crossover is nevertheless consistent with the expected mechanism: increasing graph and factor work improves the amortization of the Rank-axis reduction. Across all three families, a material benefit from additional devices generally appears only after local operator work dominates collective latency. Small instances often favor one or two GPUs, whereas larger instances with many cones, many constraint rows, or wide factors benefit from the axis aligned with their dominant structural cost. Thus neither the device count nor the mesh shape is uniformly best across problem families.

7

Conclusion

We presented CARDAL, a distributed low-rank solver for large-scale, block-structured SDPs. At fixed ranks, the method uses matrix-free L-BFGS–NC iterations to target an approximate Euclidean second-order stationary point of the factored augmented Lagrangian. The outer staircase then tests the shifted dual slack: a negative Rayleigh quotient drives blockwise rank growth, while satisfaction of the slack tolerance supports termination. The reverse multiplier shift makes a negative-slack witness an exact negative-curvature direction in a newly appended zero column. The resulting update admits a closed-form shared amplitude and a joint rank-lift problem for batches of directions. Our analysis establishes three distinct conclusions. Under the stated fixed-rank ALM and product-manifold assumptions, accumulation points are second-order critical for the equality-constrained BM problem. For almost every product cost, a per-block rank condition then yields global SDP optimality and finite exact blockwise rank growth. At finite accuracy, blockwise slack lower bounds provide an approximate KKT certificate that is deterministic and verifiable a posteriori. Under independent full-dimensional perturbations of the cost blocks, a blockwise tube argument also gives a conditional low-rank finite-accuracy result and transfers its sampled-cost bounds back to the nominal costs. None of these statements supplies a polynomial iteration bound for the implemented ALM–L-BFGS–NC procedure. For scalable execution, CARDAL distributes constraint rows, factor columns, and PSD blocks over a Constraint × Rank × Cone device mesh. The distributed residual, gradient, Hessian action, and slack matrix–vector products are exact rearrangements of their single-device counterparts in exact arithmetic. Workload-aware partitioning, batching of small cones, and PSD-preserving scaling complement this decomposition. On the Mittelmann benchmark, CARDAL passes the uniform DIMACS test on 65 of 75 instances. The three multi-GPU studies show an empirical crossover: additional devices help once the local cone, constraint, or factor work is large enough to amortize collective latency, whereas small instances can be faster on fewer GPUs. The reported chemistry and Max-Cut checkpoint timings measure throughput rather than final solution accuracy. Limitations. The generic product-cone result excludes a measure-zero set of block costs, while the finite-accuracy low-rank result requires independent smoothing of the cost blocks, prespecified rank profiles and tolerance envelopes, and per-block rank conditions. It therefore does not certify an unsmoothed run on an arbitrary fixed cost. Performance also depends on ALM conditioning, workload balance, and the adaptive rank and iteration trajectories, which can differ across device meshes. Removing the smoothing requirement and improving topology-aware balancing and high-penalty preconditioning are natural directions for future work.

References [1] Jacob M. Aguirre, Diego Cifuentes, Vincent Guigues, Renato D. C. Monteiro, Victor Hugo Nascimento, and Arnesh Sujanani. cuHALLaR: A GPU-accelerated low-rank augmented Lagrangian method for large-scale semidefinite programming. arXiv preprint arXiv:2505.13719, 2025. [2] Alexander I. Barvinok. Problems of distance geometry and convex properties of quadratic maps. Discrete & Computational Geometry, 13(2):189–202, 1995. [3] Dimitri P. Bertsekas. Constrained Optimization and Lagrange Multiplier Methods. Computer Science and Applied Mathematics. Academic Press, New York, 1982. ISBN 0-12-093480-9. [4] Nicolas Boumal, Vladislav Voroninski, and Afonso S. Bandeira. Deterministic guarantees for Burer–Monteiro factorizations of smooth semidefinite programs. Communications on Pure and Applied Mathematics, 73(3):581–608, 2020. [5] Stephen Boyd and Lieven Vandenberghe. Convex Optimization. Cambridge University Press, 2004. [6] Samuel Burer and Changhui Choi. Computational enhancements in low-rank semidefinite programming. Optimisation Methods and Software, 21(3):493–512, 2006. [7] Samuel Burer and Renato D. C. Monteiro. A projected gradient algorithm for solving the Max-Cut SDP relaxation. Optimization Methods and Software, 15(3-4):175–200, 2001.

20

P REPRINT. U NDER REVIEW.

[8] Samuel Burer and Renato D. C. Monteiro. A nonlinear programming algorithm for solving semidefinite programs via lowrank factorization. Mathematical Programming, 95(2):329–357, 2003. [9] Andrew M Childs, Andrew J Landahl, and Pablo A Parrilo. Quantum algorithms for the ordered search problem via semidefinite programming. Physical Review A: Atomic, Molecular, and Optical Physics, 75(3):032335, 2007. [10] Diego Cifuentes and Ankur Moitra. Polynomial time guarantees for the Burer–Monteiro method. Advances in Neural Information Processing Systems, 35:23923–23935, 2022. [11] Lijun Ding and Benjamin Grimmer. Revisiting spectral bundle methods: Primal-dual (sub)linear convergence rates. SIAM Journal on Optimization, 33(2):1305–1332, 2023. [12] Lijun Ding, Haihao Lu, and Jinwen Yang. New understandings and computation on augmented Lagrangian methods for low-rank semidefinite programming. arXiv preprint arXiv:2505.15775, 2025. [13] Dongdong Ge, Qi Huangfu, Zizhuo Wang, Jian Wu, and Yinyu Ye. Cardinal optimizer (COPT) user guide. arXiv preprint arXiv:2208.14314, 2022. [14] Michel X. Goemans and David P. Williamson. Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming. Journal of the ACM, 42(6):1115–1145, 1995. [15] Paul J Goulart and Yuwen Chen. Clarabel: An interior-point solver for conic programs with quadratic objectives. Mathematical Programming Computation, pages 1–83, 2026. [16] Qiaoming Han, Yinyu Ye, and Jiawei Zhang. An improved rounding method and semidefinite programming relaxation for graph partition. Mathematical Programming, 92(3):509–535, 2002. [17] Qiushi Han, Chenxi Li, Zhenwei Lin, Caihua Chen, Qi Deng, Dongdong Ge, Huikang Liu, and Yinyu Ye. A low-rank ADMM splitting approach for semidefinite programming. arXiv preprint arXiv:2403.09133, 2024. [18] Qiushi Han, Zhenwei Lin, Hanwen Liu, Caihua Chen, Qi Deng, Dongdong Ge, and Yinyu Ye. Accelerating low-rank factorization-based semidefinite programming algorithms on GPU. arXiv preprint arXiv:2407.15049, 2024. [19] Shucheng Kang, Xiaoyang Xu, Jay Sarva, Ling Liang, and Heng Yang. Fast and certifiable trajectory optimization. In International Workshop on the Algorithmic Foundations of Robotics, pages 43–65. Springer, 2024. [20] Shucheng Kang, Guorui Liu, and Heng Yang. Global contact-rich planning with sparsity-rich semidefinite relaxations. arXiv preprint arXiv:2502.02829, 2025. [21] Guanghui Lan and Renato D. C. Monteiro. Iteration-complexity of first-order augmented Lagrangian methods for convex programming. Mathematical Programming, 155(1):511–547, 2016. [22] Jean B Lasserre. Global optimization with polynomials and the problem of moments. SIAM Journal on Optimization, 11(3): 796–817, 2001. [23] Jean B Lasserre. Convergent SDP-relaxations in polynomial optimization with sparsity. SIAM Journal on Optimization, 17 (3):822–843, 2006. [24] Jure Leskovec and Andrej Krevl. SNAP Datasets: Stanford large network dataset collection. http://snap.stanford. edu/data, 2014. [25] 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. [26] Feng-Yi Liao, Lijun Ding, and Yang Zheng. An overview and comparison of spectral bundle methods for primal and dual semidefinite programs. Computational Optimization and Applications, 93(2):523–566, 2026. [27] Guang Hao Low, Robbie King, Dominic W Berry, Qiushi Han, A Eugene DePrince III, Alec F White, Ryan Babbush, Rolando D Somma, and Nicholas C Rubin. Fast quantum simulation of electronic structure by spectral amplification. Physical Review X, 15(4):041016, 2025. [28] Hans D. Mittelmann. Benchmarks for optimization software: Sparse SDP. https://plato.asu.edu/ftp/sparse_sdp. html, 2026. Accessed 2026. [29] Renato D. C. Monteiro, Arnesh Sujanani, and Diego Cifuentes. A low-rank augmented Lagrangian method for large-scale semidefinite programming based on a hybrid convex–nonconvex approach. Mathematical Programming, pages 1–57, 2026. [30] MOSEK ApS. MOSEK Optimization Suite. MOSEK ApS, 2019. [31] 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(3):1042–1068, 2016. [32] Gábor Pataki. On the rank of extreme matrices in semidefinite programs and the multiplicity of optimal eigenvalues. Mathematics of Operations Research, 23(2):339–358, 1998. [33] Thomas Pock and Antonin Chambolle. Diagonal preconditioning for first-order primal-dual algorithms in convex optimization. In 2011 International Conference on Computer Vision (ICCV), pages 1762–1769, 2011. [34] R. Tyrrell Rockafellar. Augmented Lagrangians and applications of the proximal point algorithm in convex programming. Mathematics of Operations Research, 1(2):97–116, 1976. [35] Stephen B. Seidman. Network structure and minimum degree. Social Networks, 5(3):269–287, 1983. [36] Jos F. Sturm. Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones. Optimization Methods and Software, 11(1-4):625–653, 1999. [37] Defeng Sun, Kim-Chuan Toh, Yancheng Yuan, and Xin-Yuan Zhao. SDPNAL+: A MATLAB software for semidefinite programming with bound constraints (version 1.0). Optimization Methods and Software, 35(1):87–115, 2020. [38] Michael J. Todd. Semidefinite optimization. Acta Numerica, 10:515–560, 2001. [39] Kim-Chuan Toh, Michael J Todd, and Reha H Tütüncü. On the implementation and usage of SDPT3–a MATLAB software package for semidefinite-quadratic-linear programming, version 4.0. In Handbook on semidefinite, conic and polynomial optimization, pages 715–754. Springer, 2012.

21

P REPRINT. U NDER REVIEW.

[40] Jie Wang and Liangbing Hu. Solving low-rank semidefinite programs via manifold optimization. Journal of Scientific Computing, 104(1):33, 2025. [41] Yifei Wang, Kangkang Deng, Haoyang Liu, and Zaiwen Wen. A decomposition augmented Lagrangian method for low-rank semidefinite programming. SIAM Journal on Optimization, 33(3):1361–1390, 2023. [42] Liuqin Yang, Defeng Sun, and Kim-Chuan Toh. SDPNAL+: a majorized semismooth Newton–CG augmented Lagrangian method for semidefinite programming with nonnegative constraints. Mathematical Programming Computation, 7(3):331– 366, 2015. [43] Alp Yurtsever, Olivier Fercoq, Francesco Locatello, and Volkan Cevher. A conditional gradient framework for composite convex minimization with applications to semidefinite programming. In International Conference on Machine Learning, pages 5727–5736. PMLR, 2018. [44] Alp Yurtsever, Olivier Fercoq, and Volkan Cevher. A conditional-gradient-based augmented Lagrangian framework. In International Conference on Machine Learning, pages 7272–7281. PMLR, 2019. [45] Alp Yurtsever, Joel A Tropp, Olivier Fercoq, Madeleine Udell, and Volkan Cevher. Scalable semidefinite programming. SIAM Journal on Mathematics of Data Science, 3(1):171–200, 2021. [46] Xin-Yuan Zhao, Defeng Sun, and Kim-Chuan Toh. A Newton–CG augmented Lagrangian method for semidefinite programming. SIAM Journal on Optimization, 20(4):1737–1765, 2010.

A

Additional Proofs

A.1

Product-cone stationarity and landscape

Proof of Proposition 1. Let F ∈ Mk be a feasible full-space second-order stationary point. Feasibility gives yb = y, and (18) gives Sc (y)Fc = 0, c ∈ [q]. (85) Let µB = µB (F ) and δ := y − µB . The compatibility relation in (85) implies CF = BF y, and therefore µB = G(F )† BF∗ CF = G(F )† G(F )y. Because G(F )† G(F ) is the orthogonal projector onto range G(F ),
q A∗c (δ)Fc c=1 = 0.

δ ∈ ker G(F ) = ker BF , It follows immediately that SB,c (F )Fc = 0 for every block.

For U ∈ TF Mk , Assumption 1 provides a smooth feasible curve Γ(t) = (Γc (t))c through F with velocity U . Twice differentiating ⟨δ, c(Γ(t))⟩ ≡ 0 and using BF δ = 0 yields q X

⟨Uc , A∗c (δ)Uc ⟩ = 0.

c=1

Consequently, X X ⟨Uc , Sc (y)Uc ⟩ = ⟨Uc , SB,c (F )Uc ⟩. c

(86)

c

The penalty term in (20) vanishes because Dc(F )[U ] = 0. Full-space second-order stationarity and (86) therefore give the curvature inequality in (15). Together with block stationarity, this proves the proposition. Proof of Theorem 1. For each block with kc < nc , let Vc,kc := {H ∈ Snc : rank(H) ≤ nc − kc }. This determinantal variety has codimension τ (kc ) in Snc . Since dim Im(A∗c ) = rc ,
dim Vc,kc + Im(A∗c ) ≤ τ (nc ) − τ (kc ) + rc < τ (nc ). It is therefore a measure-zero subset of Snc . Suppose a first-order critical factor has a full-column-rank block Fc with kc < nc , and set λ := µB (F ). Stationarity Sc (λ)Fc = 0 implies rank Sc (λ) ≤ nc − kc , and hence Cc = Sc (λ) + A∗c (λ) ∈ Vc,kc + Im(A∗c ). Q The union of the corresponding cylinders in c Snc is measure zero. Outside that union, every block with kc < nc is column-rank deficient.

22

P REPRINT. U NDER REVIEW.

Now let F be second-order critical. If Fc is column-rank deficient, choose 0 ̸= zc ∈ ker Fc . For arbitrary uc ∈ Rnc , let U be zero outside block c and set Uc = uc zc⊤ . Then Dc(F )[U ] = 0, and (15) gives 0 ≤ ⟨Uc , SB,c (F )Uc ⟩ = ∥zc ∥22 u⊤ c SB,c (F )uc . Thus SB,c (F ) ⪰ 0. Otherwise Fc has full column rank. The generic argument excludes this case when kc < nc , and kc > nc cannot have full column rank. Hence kc = nc , Fc is invertible, and SB,c (F )Fc = 0 gives SB,c (F ) = 0. Every slack block is therefore positive semidefinite. Primal feasibility and complementarity complete the product-SDP KKT conditions. A.2

Finite-accuracy certificates and smoothing P Proof of Proposition 2. Let r := c Ac (Xc ) − b. Since Cc = A∗c (λ) + Sc (λ), X X ⟨Cc , Xc ⟩ = ⟨λ, b + r⟩ + tr(Sc (λ)Xc ). c

c

Cauchy–Schwarz across the block traces gives X

tr(Sc (λ)Xc ) ≤

s X

c

!1/2 nc

c

X

∥Sc (λ)Xc ∥2F

.

c

Hence X ⟨Cc , Xc ⟩ ≤ b⊤ λ + ε0 ∥λ∥2 +

s

X

nc ε1 .

c

c

For any feasible comparison tuple (Xc⋆ )c , X X X X b⊤ λ = ⟨Cc , Xc⋆ ⟩ − ⟨Sc (λ), Xc⋆ ⟩ ≤ ⟨Cc , Xc⋆ ⟩ + ε2 tr(Xc⋆ ). c

c

c

c

Combining these inequalities proves (26). Proof of Proposition 3. Set λ = yb(F ; y, ρ). Equation (18) gives !1/2 X η 2 ∥Sc (λ)Fc ∥F ≤ . 2 c For an admissible unit direction in (28), ∥Dc(F )[U ]∥2 = 2

X

Ac (Uc Fc⊤ )

c

≤ 2γ. 2

The approximate Hessian bound and (20) then imply X 2 ⟨Uc , Sc (λ)Uc ⟩ ≥ −ζ − 4ργ 2 . c

This is (28) with the tolerances in (29). Proof of Theorem 2. The primal residual bound is part of (27). Moreover, !1/2 !1/2 X X 2 2 2 ∥Sc (λ)Xc ∥F ≤ ∥Sc (λ)Fc ∥F ∥Fc ∥op c

c

≤ RF ϵ1 . It remains to prove approximate dual feasibility. If kc > nc , choose a unit zc ∈ ker Fc . For any unit uc ∈ Rnc , the product direction supported on block c with Uc = uc zc⊤ satisfies Ac (Uc Fc⊤ ) = 0. Equation (28) then gives u⊤ c Sc (λ)uc ≥ −ϵ2 . The same conclusion holds when ac = 0, now using any unit zc ∈ Rkc . Hence these blocks are deterministically approximately dual feasible.

23

P REPRINT. U NDER REVIEW.

Fix c ∈ J . If σmin (Fc ) ≤ γ/ac , choose a corresponding unit right singular vector zc and again take Uc = uc zc⊤ . Then ∥Ac (Uc Fc⊤ )∥2 ≤ ac ∥Uc Fc⊤ ∥F = ac ∥Fc zc ∥2 ≤ γ, so (28) implies Sc (λ) ⪰ −ϵ2 I. Any block violating approximate dual feasibility must therefore satisfy γ σmin (Fc ) > . ac

(87)

Let Vc,kc be the determinantal variety used in the proof of Theorem 1. The min–max characterization of singular values and (87) give
∥Sc (λ)Fc ∥F ϵ1 ac < = δc . dist Sc (λ), Vc,kc ≤ σmin (Fc ) γ Thus a bad block can occur only if Cc ∈ tubeδc (Vc,kc ) + A∗c (BRλ ), (88) where BRλ is the multiplier ball. The set A∗c (BRλ ) lies in an rc -dimensional subspace and in the Frobenius ball of radius κc . Since δc ≤ κc , it has a δc -net with at most (3κc /δc )rc points. For each net point, (88) is contained in a 2δc -tube around a translate of Vc,kc . That variety has codimension τ (kc ) and is defined by minors of degree nc −kc +1. Applying the tube-volume estimate used in [10, Theorem 6], together with δc < σc /(4n3c ), bounds the probability of (88) by  r  τ (kc ) 3κc c 4n3c δc 4e = pc . δc σc A union bound over c ∈ J proves (32). Substituting (33) into (31) gives the simplified bound stated after the theorem, by the same algebra as [10, Corollary 1]. A.3

CARDAL convergence and finite outputs

Proof of Theorem 3. Lemma 2 gives ∥c(Fℓ )∥2 → 0. Let Fℓj +1 → F⋆ and, after taking a further subsequence, let yℓj +1 → y⋆ . The residual and gradient bounds imply c(F⋆ ) = 0,

Sc (y⋆ )F⋆,c = 0

(c ∈ [q]).

The product tangent projector is PF = I − BF G(F )† BF∗ . Assumption 1 makes PF smooth near Mk . Given U⋆ ∈ ker Dc(F⋆ ), set Uj = PFℓj +1 U⋆ . Then Uj → U⋆ and Dc(Fℓj +1 )[Uj ] = 0. Since yℓj +1 is the shifted multiplier of the subproblem indexed by ℓj , the inner Hessian bound and (20) give X 2 ⟨Uj,c , Sc (yℓj +1 )Uj,c ⟩ ≥ −ζℓj ∥Uj ∥2⊕ . c

Passing to the limit and applying the compatible-multiplier argument from the proof of Proposition 1 yields (15). Under the additional rank and generic-cost hypotheses, Theorem 1 makes every accumulation point globally optimal. If objective convergence failed, compactness and continuity would produce an accumulation point whose objective differs from p⋆ . If convergence in distance to X ⋆ failed, the same argument would produce a represented matrix tuple outside the closed solution set. Both contradict global optimality of every accumulation point. Proof of Corollary 1. There are finitely many block profiles between k0 and kgen , so the union of their measure-zero exceptional cost sets remains measure zero. Consider a stage at which block c has reached kcgen . Exact feasibility and first-order stationarity give Sc (y)Fc = 0. If Fc has full column rank, the generic block-deficiency argument rules out kcgen < nc , while kcgen > nc cannot have full column rank. Thus kcgen = nc , Fc is invertible, and Sc (y) = 0. Otherwise choose a nonzero zc ∈ ker Fc . With U supported only on block c and Uc = uc zc⊤ , we have Dc(F )[U ] = 0. Full-space second-order stationarity and (20) give 0 ≤ 2∥zc ∥22 u⊤ for every uc , c Sc (y)uc so again Sc (y) ⪰ 0. The exact slack test cannot enlarge this block again. Every accepted update increases at least one block that has not reached its certifying rank. Summing the maximum number of updates of each block gives (47). At termination, every slack block is positive semidefinite, while exact feasibility and stationarity give complementarity. The primal–dual KKT conditions prove global optimality.

24

P REPRINT. U NDER REVIEW.

Proof of Proposition 4. Equation (18) gives !1/2 X

∥Sc (yout )Fc ∥2F

c

≤

ω . 2

Consequently, !1/2 X

∥Sc (yout )Xc ∥2F

!1/2 ≤

X

c

∥Sc (yout )Fc ∥2F ∥Fc ∥2op

c

RF ω ≤ . 2 Primal positive semidefiniteness, the primal-residual bound, and the blockwise slack bounds are assumed directly, which proves (49). Proof of Corollary 2. Proposition 3 shows that (F, yb) is a product (ϵ, γ)-AFAC pair with the fixed tolerances in (50). Theorem 2 gives the sampled-cost conclusion. Write Cc = C̄c + Wc . Since ∥Wc ∥op ≤ ∥Wc ∥F ≤ σc , y ) ⪰ −(ϵ2 + σmax )I. SC̄c (b P The direct-sum triangle inequality, ∥Xc ∥F ≤ ∥Fc ∥2F , and c ∥Fc ∥2F ≤ RF2 also give !1/2 X y )Xc ∥2F ≤ RF ϵ1 + σmax RF2 . ∥SC̄c (b c

The affine constraints are unchanged, proving (53). Finally, the sampled bound SCc (b y ) ⪰ −ϵ2 I excludes a Rayleigh quotient below −τdual whenever τdual ≥ ϵ2 .

25