Paper deep dive

Reinforcement Learning to Accelerate Primal-Dual Hybrid Gradient for Linear Programming

Jinhwan Sul, Alex Oshin, Evangelos A. Theodorou

Year: 2026Venue: arXiv preprintArea: math.OCType: PreprintEmbeddings: 146

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%

Last extracted: 10/3/2026, 3:31:06 AM

Summary

The paper introduces GALLOP, a reinforcement learning framework that accelerates Primal-Dual Hybrid Gradient (PDHG) methods for large-scale linear programming. GALLOP uses a dimension-agnostic policy to jointly learn continuous algorithm parameters (step size, primal weight) and discrete restart decisions within a generalized accelerated PDHG (GA-PDHG) update scheme. Evaluated on six LP families, GALLOP achieves significant speedups (up to 16x) and iteration reductions (1.9x-5.6x) over MPAX, with policies generalizing to problems up to 400x larger than training instances without retraining.

Entities (8)

Linear Programming · problemclass · 99%GALLOP · algorithm · 98%PDHG · algorithm · 95%Reinforcement Learning · method · 95%GA-PDHG · algorithm · 92%MPAX · algorithm · 90%PPO · algorithm · 90%Transport LPs · dataset · 85%

Relation Signals (5)

PDHG → solves → Linear Programming

confidence 95% · Primal-dual hybrid gradient (PDHG) methods solve large-scale linear programs (LPs)

GALLOP → uses → Reinforcement Learning

confidence 95% · We introduce GALLOP, which uses reinforcement learning to jointly learn continuous algorithm parameters

GALLOP → implements → GA-PDHG

confidence 92% · We augment PDHG with a generalized acceleration step and a restart mechanism, and refer to the resulting update scheme as generalized accelerated PDHG (GA-PDHG)... We call our proposed algorithm GALLOP

GALLOP → outperforms → MPAX

confidence 90% · achieves up to a 16.0× speedup in algorithm wall-clock time over MPAX.

GALLOP → generalizesto → Transport LPs

confidence 85% · the learned policies generalize without retraining to within-family LPs... including Transport LPs with 10.24 million variables.

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Primal-dual hybrid gradient (PDHG) methods solve large-scale linear programs (LPs) using GPU-friendly matrix-vector products and projections, but their practical performance depends on coordinating algorithm parameters, acceleration, and restarts. We introduce GALLOP, which uses reinforcement learning to jointly learn continuous algorithm parameters and discrete restart decisions without differentiating through the solver. Its generalized accelerated PDHG update combines separate primal and dual extrapolation, history corrections, and restart anchoring with independently adjustable coefficients. We train a dimension-agnostic feedback policy using a groupwise proximal policy optimization objective that clips likelihood ratios separately for different control groups and excludes inactive acceleration controls on restart transitions. We evaluate GALLOP on six LP families and a public item-placement benchmark. On the main evaluation settings across the six families, GALLOP reduces iteration counts by factors of $1.9$-$5.6$ and achieves up to a $16.0\times$ speedup in algorithm wall-clock time over MPAX. With one policy trained per family, the learned policies generalize without retraining to within-family LPs $3\times$-$400\times$ larger than the largest training instances, including Transport LPs with $10.24$ million variables.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

142,516 characters extracted from source content.

Expand or collapse full text

Reinforcement Learning to Accelerate Primal-Dual Hybrid Gradient for Linear Programming Jinhwan Sul Affiliation: Georgia Institute of Technology, USA Email: jsul7@gatech.edu Alex Oshin Affiliation: Georgia Institute of Technology, USA Email: alexoshin@gatech.edu Evangelos A. Theodorou Affiliation: Georgia Institute of Technology, USA Email: evangelos.theodorou@gatech.edu Abstract Primal–dual hybrid gradient (PDHG) methods solve large-scale linear programs (LPs) using GPU-friendly matrix–vector products and projections, but their practical performance depends on coordinating algorithm parameters, acceleration, and restarts. We introduce GALLOP, which uses reinforcement learning to jointly learn continuous algorithm parameters and discrete restart decisions without differentiating through the solver. Its generalized accelerated PDHG update combines separate primal and dual extrapolation, history corrections, and restart anchoring with independently adjustable coefficients. We train a dimension-agnostic feedback policy using a groupwise proximal policy optimization objective that clips likelihood ratios separately for different control groups and excludes inactive acceleration controls on restart transitions. We evaluate GALLOP on six LP families and a public item-placement benchmark. On the main evaluation settings across the six families, GALLOP reduces iteration counts by factors of 1.91.9–5.65.6 and achieves up to a 16.0×16.0× speedup in algorithm wall-clock time over MPAX. With one policy trained per family, the learned policies generalize without retraining to within-family LPs 3×3×–400×400× larger than the largest training instances, including Transport LPs with 10.2410.24 million variables. 1 Introduction Linear programming (LP) is a fundamental problem class in constrained optimization (Dantzig, 1963). The simplex and interior-point methods are widely used to solve LPs (Dantzig, 1963; Wright, 1997). Recently, due to the scalability limitations of these traditional methods (Applegate et al., 2026), primal–dual hybrid gradient (PDHG) method (Chambolle & Pock, 2011) based LP algorithms have been explored for large-scale problems (Applegate et al., 2021; Lu & Yang, 2025; Lu et al., 2025). PDHG iteration consists of matrix–vector products and projections which are inexpensive and GPU-friendly. However, PDHG can require a large number of iterations to attain high-accuracy solutions (Applegate et al., 2023). PDHG-based LP solvers (Lu et al., 2025; Zhou et al., 2025) improve their practical convergence through heuristic adjustment of algorithmic parameters, such as the step size and primal weight, and acceleration techniques, such as restarted Halpern acceleration (Halpern, 1967; Lu & Yang, 2024) and Anderson acceleration (Anderson, 1965). Many of these heuristics and acceleration schemes are motivated by theoretical convergence analyses (Lu & Yang, 2024; Applegate et al., 2023; Zhou et al., 2025). However, existing global convergence guarantees do not cover the full combination of heuristic adaptations used in practical solvers, whose effectiveness has been demonstrated empirically (Lu et al., 2023; Lu et al., 2025). Beyond convergence guarantees, a separate design challenge is to determine how these adaptations should be selected and coordinated to reduce solution time for a given family of LPs. This leaves an opportunity to further accelerate PDHG-based LP solvers by learning algorithm parameters and heuristic rules from LP data. Learning to optimize (L2O) uses problem data to accelerate optimization (Gregor & LeCun, 2010). We focus on L2O methods that use a learned policy throughout the solve until convergence, without handing the remaining iterations to an external solver. Venkataraman & Amos (2021) and Martin & Belgioioso (2026) use algorithm unrolling for fixed-point optimization. However, their neural network designs are tied to the problem dimension and therefore do not support size generalization. PDHG-Net (Li et al., 2024) unrolls the basic PDHG iteration without incorporating discrete restart decisions or the additional acceleration schemes used by practical PDHG-based solvers. It does not run the learned updates until convergence, but instead uses the network output to warm-start PDLP. RLQP (Ichnowski et al., 2021) uses reinforcement learning (RL) to learn the penalty parameter in OSQP (Stellato et al., 2020). While effective, RLQP applies model-free RL to tune continuous parameters of an iterative algorithm whose updates can be differentiated. In this setting, deep unfolding enables gradients of the optimization objective to be propagated directly through the solver, avoiding the additional value-function estimation and exploration required by RL. For example, DeepQP (Saravanos et al., 2025) unfolds OSQP and learns its algorithm parameters end-to-end, providing a more direct gradient-based alternative. Figure 1: Conceptual comparison of (a) vanilla PDHG, (b) PDHG with heuristic acceleration and restarts, and (c) GALLOP with learned solver parameters and restarts. To solve large LPs faster, we learn a policy that jointly controls continuous algorithm parameters and discrete restart decisions in a PDHG-based LP algorithm. We use RL to learn these controls jointly because ordinary backpropagation through the solver does not directly provide gradients for hard restart decisions. The policy uses a fixed number of solver-state features and outputs a fixed number of controls, independent of the problem size. This dimension-agnostic design allows us to train on smaller LPs and apply the learned policy to larger LPs without retraining. We augment PDHG with a generalized acceleration step and a restart mechanism, and refer to the resulting update scheme as generalized accelerated PDHG (GA-PDHG). This step combines separate primal and dual extrapolation along the current PDHG update, history corrections from recent iterate differences, and a correction toward the restart anchor. We learn the coefficients of these terms together with the step size, primal weight, and restart decision. These controls provide a broad action space for RL to explore combinations of acceleration and restarting. We call our proposed algorithm GALLOP, short for Generalized Acceleration and Reinforcement Learning for Linear Optimization with PDHG. In experiments on six LP families, we demonstrate GALLOP is 1.9×1.9×–5.6×5.6× faster in terms of number of iterations relative to MPAX’s restarted reflected Halpern solver (Lu et al., 2024). Algorithm wall-clock speedups reach 16.0×16.0× on Random LPs with 300,000300,000 variables, 10.1×10.1× on Auction LPs with 358,400358,400 variables, and 3.7×3.7× on Transport LPs with 10.2410.24 million variables. Our trained policies generalize to problems 3×3×–400×400× larger than the largest training size in a zero shot setting, using the family-specific size measures. In particular, the Transport policy trained on LPs with at most 25,60025,600 variables, scales to LPs with 10.2410.24 million variables, a 400×400× size generalization without retraining. Ablations on Transport show that all five learned control groups contribute to reducing iteration counts, including the generalized-acceleration coefficients. Our main contributions are as follows. • We propose GA-PDHG, a generalized accelerated PDHG update scheme that combines independently controlled primal and dual extrapolation, history corrections, and a restart-anchor correction. This parameterization encompasses established acceleration schemes and their combinations, providing a broad search space for RL to discover acceleration strategies. • We introduce GALLOP, an RL-based L2O framework for large-scale LPs. To accommodate the conditional discrete-continuous action space of GA-PDHG, we develop a groupwise PPO objective that clips likelihood ratios separately for each control group and masks inactive acceleration controls on restart transitions. • We demonstrate acceleration and size generalization across six LP families, with iteration-count reductions by factors of 1.91.9–5.65.6 and up to a 16.0×16.0× speedup in algorithm wall-clock time over MPAX. The learned policies generalize to LPs 3×3×–400×400× larger than the largest training instances under family-specific size measures, including Transport LPs with 10.2410.24 million variables, without retraining. 2 Primal–Dual Hybrid Gradient for Large-Scale Linear Programming PDLP (Applegate et al., 2021) is a first-order LP solver based on the primal–dual hybrid gradient (PDHG) method (Chambolle & Pock, 2011). Its updates use matrix–vector products and projections, avoiding matrix factorization. Basic PDHG can require many iterations to reach high accuracy. When PDHG is applied to LP, acceleration techniques and heuristics are introduced to reduce the number of iterations. This section presents the LP formulation and PDHG updates and reviews these enhancements in different PDLP variants. 2.1 Problem Formulation: Linear Programming We consider the LP minx∈Xc⊤xsubject toPx=q,Gx≥h,X=x:l≤x≤u, _x∈ X c x to Px=q, Gx≥ h, X=\x:l≤ x≤ u\, (1) where x∈ℝnx ^n is the primal variable and c∈ℝnc ^n is the objective coefficient vector. The LP has mem_e equality constraints specified by P∈ℝme×nP ^m_e× n and q∈ℝmeq ^m_e, and mim_i inequality constraints specified by G∈ℝmi×nG ^m_i× n and h∈ℝmih ^m_i. The primal variable x is bounded with X where l∈(ℝ∪−∞)nl∈ (R∪\-∞\ )^n, u∈(ℝ∪+∞)nu∈ (R∪\+∞\ )^n, and l≤ul≤ u componentwise. Define A=[P⊤G⊤]⊤∈ℝm×nA=[P \;G ] ^m× n and b=[q⊤h⊤]⊤∈ℝmb=[q \;h ] ^m, where m=me+mim=m_e+m_i. The corresponding dual vector is y=[(yeq)⊤(yin)⊤]⊤∈ℝmy=[ (y^eq ) \; (y^in ) ] ^m, with equality multipliers yeq∈ℝmey^eq ^m_e and inequality multipliers yin∈ℝ+miy^in _+^m_i. Linear programs (LPs) are commonly solved with the simplex (Dantzig, 1963) or interior-point methods (Wright, 1997). Standard implementations of these methods include matrix factorizations, whose computational and memory costs can become prohibitive at large scales. 2.2 Primal–Dual Hybrid Gradient for Linear Programming PDHG uses proximal gradients of the Lagrangian to update the primal and dual variables without matrix factorization. To derive these first-order updates, we express Eq. (1) as the equivalent convex–concave min-max problem minx∈X⁡maxy∈YL⁡(x,y):=c⊤x+(yeq)⊤(q−Px)+(yin)⊤(h−Gx)=c⊤x−y⊤Ax+b⊤y, _x∈ X\ _y∈ Y L(x,y):=c x+ (y^eq ) (q-Px )+ (y^in ) (h-Gx )=c x-y Ax+b y, (2) where Y:=ℝme×ℝ+miY:=R^m_e×R_+^m_i leaves the equality multipliers unrestricted and the inequality multipliers nonnegative. For fixed y∈Yy∈ Y, minimizing the Lagrangian over X gives the dual objective infx∈XL⁡(x,y)=b⊤y+infx∈X(c−A⊤y)⊤x. _x∈ XL (x,y )=b y+ _x∈ X (c-A y ) x. (3) For the dual objective to be finite, the reduced cost c−A⊤yc-A y must satisfy the sign conditions imposed by X. We denote the set of such reduced costs by Λ . Throughout the paper, z:=(x⊤,y⊤)⊤∈ℝn+mz:=(x ,y ) ^n+m denotes a stacked primal–dual vector. Subscripts and accents on z apply to both blocks x and y. At iteration t, z^t z_t is the PDHG input and zt∈X×Yz_t∈ X× Y is its output. At iteration t, projected descent and ascent with primal extrapolation give xt x_t =projX⁡[x^t−ηω(c−A⊤y^t)], =proj_X [ x_t- ηω (c-A y_t ) ], (4a) yt y_t =projY⁡[y^t+ωη(b−A⁡(2xt−x^t))]. =proj_Y [ y_t+ωη (b-A (2x_t- x_t ) ) ]. (4b) Here projproj denotes Euclidean projection, η>0η>0 is the step size, and ω>0ω>0 is the primal weight, giving primal and dual step sizes η/ωη/ω and ηωηω, respectively. Writing ℋH for this PDHG map, we have zt=ℋ⁡(z^t,ω,η)z_t=H ( z_t;ω,η ). With iteration-dependent parameters ηt,ωt>0 _t, _t>0, zt+1=ℋ⁡(z^t+1,ωt,ηt)z_t+1=H ( z_t+1; _t, _t ). Thus the current pair (z^t,zt) ( z_t,z_t ) was computed with (ωt−1,ηt−1) ( _t-1, _t-1 ), with (ω0,η0) ( _0, _0 ) reserved for initialization. For any ω,η>0ω,η>0, z is a fixed point of the PDHG map, i.e. ℋ⁡(z,ω,η)=zH (z;ω,η )=z, if and only if z satisfies the KKT conditions of Eq. (1). PDHG enhancements for LP. PDHG-based LP solvers choose a constant step size or adaptively choose it through trial-update acceptance tests (Chambolle & Pock, 2011; Applegate et al., 2021). The primal weight is fixed or adjusted at restarts or fixed intervals to balance primal–dual displacements, typically through exponential moving averaging (EMA) in log space or proportional–integral–derivative (PID) control (Lu et al., 2025). Restarting (Hinder & Lubin, 2020) resets the averaging or acceleration state from a selected iterate after sufficient progress, a loss of progress following an initial decrease, or an excessively long epoch. The progress is measured by the normalized duality gap, KKT error, or fixed-point residual (Applegate et al., 2023; Lu & Yang, 2024). Acceleration uses step-size-weighted averaging (Applegate et al., 2021), Halpern anchoring with a diminishing anchor weight and optional reflection, or Anderson history corrections computed by least-squares residual minimization (Anderson, 1965; Zhou et al., 2025). With suitable parameter choices and restart rules, restarted averaging and Halpern acceleration exploit LP sharpness to achieve accelerated linear convergence (Applegate et al., 2023; Lu & Yang, 2024). See Appendix C for the details on these enhancements across representative solvers. 3 Reinforcement Learning for PDHG with Acceleration The enhancement rules reviewed in PDHG-based LP solvers are hand-designed feedback controls based on solver state. Their criteria address individual decisions, such as balancing primal–dual displacements, testing restart thresholds, or minimizing a local residual model. These heuristics motivate exploring whether a more expressive policy that coordinates these decisions using a broader set of solver-state features can further accelerate convergence. We therefore learn such a policy from solver trajectories using a common reward that induces faster convergence. 3.1 PDHG with Generalized Acceleration To enable a learned policy to explore both established acceleration rules and their combinations, we introduce a family of accelerated PDHG updates. Specifically, we combine extrapolation, history corrections, and anchoring with independently adjustable coefficients to formulate generalized accelerated PDHG (GA-PDHG). For a fixed, finite memory length k≥0k≥ 0, define the output-history matrix t:=[zt−zt−1⋯zt−k+1−zt−k]∈ℝ(n+m)×kK_t:= bmatrixz_t-z_t-1&·s&z_t-k+1-z_t-k bmatrix (n+m )× k. Unavailable differences and those discarded at a restart are replaced by zero. The vector z⋆z_ is the fixed anchor within the current epoch, and s counts the non-restart iterations within the current epoch. The coefficients αtx,αty≥−1 _t^x, _t^y≥-1 control primal and dual relaxation or extrapolation, t∈ℝk β_t ^k controls the history correction, and ζt∈ℝ _t controls the anchor correction. Let IdI_d denote the d-dimensional identity matrix and et∈0,1e_t∈\0,1\ the restart decision. When et=0e_t=0 (no restart), the GA-PDHG update is z^t+1 z_t+1 =zt+(αtxIn00αtyIm)⁡(zt−z^t)+tt+ζt(z⋆−zt), =z_t+ pmatrix _t^xI_n&0\\ 0& _t^yI_m pmatrix (z_t- z_t )+K_t β_t+ _t (z_ -z_t ), (5a) zt+1 z_t+1 =ℋ⁡(z^t+1,ωt,ηt). =H ( z_t+1; _t, _t ). (5b) When et=1e_t=1 (restart), the restart branch sets z^t+1=zt z_t+1=z_t. Note that GA-PDHG always ends with a PDHG mapping and takes its output as the iterate, allowing the acceleration coefficients to explore combinations beyond convex combinations and generate proposals outside X×YX× Y while ensuring zt+1∈X×Yz_t+1∈ X× Y. After a restart, we set z⋆←zt+1z_ ← z_t+1, t+1←0K_t+1← 0, and s←0s← 0. Table 1 shows how coefficient choices in Eq. (5a) correspond to five acceleration and relaxation schemes for fixed-point problems. GA-PDHG thus offers a unified search space that encompasses these schemes as special cases. For general coefficient choices, Proposition 3.1 gives sufficient conditions under which the output of a convergent GA-PDHG trajectory is primal–dual optimal in the limit. Note that we do not claim global convergence for arbitrary learned control sequences. See Appendix B for details. Table 1: Acceleration and relaxation schemes for fixed-point problems via GA-PDHG coefficients. Method αtx=αty=αt _t^x\!=\! _t^y\!=\! _t t β_t ζt _t Halpern (Halpern, 1967) 00 0 1s+2 1s+2 Reflected Halpern (Lu & Yang, 2024) s+1s+2 s+1s+2 0 1s+2 1s+2 Anderson (Anderson, 1965) 00 tAA β_t^A 00 Krasnosel’skiĭ–Mann (KM) (Krasnosel’skii, 1955; Mann, 1953) λ−1λ-1 0 00 Fast KM (Boţ & Nguyen, 2023) ν2(s+ν)−1 ν2 (s+ν )-1 s+ν ss+ν 00 tAA β_t^A is obtained by solving a least squares problem. λ∈(0,2)λ\!∈\!(0,2). For Fast KM, k=1k=1, ν>2ν>2 and s≥1s≥ 1. Proposition 3.1 (Optimality of convergent GA-PDHG). Suppose t β_t is bounded, ηt,ωt _t, _t are uniformly positive and bounded, and there exists δ>0δ>0 such that 1+αtx,1+αty≥δ1+ _t^x,1+ _t^y≥δ whenever et=0e_t=0. If both z^t\ z_t\ and zt\z_t\ converge to finite limits and ‖ζt(z⋆−zt)‖2⟶0 _t (z_ -z_t ) _2 0 along the trajectory, then the limit of zt\z_t\ is a primal–dual optimal solution of Eq. (1). Selecting the coefficients of GA-PDHG can be viewed as an optimal control problem in which the GA-PDHG iteration defines a discrete-time dynamical system and the objective is to accelerate convergence. Figure 2 illustrates a feedback control loop in which the solver inputs are selected based on the current solver state. At iteration t, the solver state StS_t contains the LP data, input–output pair (zt,z^t)(z_t, z_t), anchor z⋆z_ , and history tK_t, the previous control variables, convergence records, and iteration counters. From this state, we construct an observation ot=Obs⁡(St)o_t=Obs (S_t ) whose dimension is independent of the problem size. A policy network inputs oto_t and outputs an action ata_t comprising the restart decision ete_t and GA-PDHG coefficients αtx,αty,t,ζt _t^x, _t^y, β_t, _t, and PDHG parameters ωt,ηt _t, _t. Using this action, the solver performs the GA-PDHG update and updates its diagnostics and counters to obtain St+1S_t+1. The next observation ot+1=Obs⁡(St+1)o_t+1=Obs (S_t+1 ) is then constructed and passed to the policy for the next decision. Figure 2: Feedback control loop for GA-PDHG. Red symbols denote control variables selected by the policy network from solver observations. 3.2 Reinforcement Learning Formulation We use reinforcement learning (RL) to train the feedback policy πθ _θ for GA-PDHG and refer to the resulting solver as GALLOP. Each LP solve corresponds to an episode, during which the policy selects solver controls at each iteration to maximize the expected cumulative reward. To support size generalization, we design the observation and control spaces to have dimensions independent of n and m, allowing the same policy to be evaluated on LPs of different sizes. At each decision step, the policy receives the observation ot=Obs⁡(St)o_t=Obs(S_t), which summarizes the solver state. Log-normalized primal and dual residuals, the duality gap, their changes, and fixed-point residual progress describe current error and recent improvement. Iteration count, restart age, progress from the anchor, primal–dual displacement balance, relative primal weight, and its one-step change characterize the current epoch. A normalized history Gram matrix, inner products with the PDHG direction, and anchor magnitude and alignment describe the scales and relative directions of the available corrections. Given this observation, the policy selects controls in a conditional hybrid action space, in which the discrete restart decision determines which continuous controls affect the solver transition. The policies for ete_t, ηt _t, and ωt _t are trained using both restart and non-restart transitions, whereas the policies for the acceleration coefficients (αtx,αty,t,ζt) ( _t^x, _t^y, β_t, _t ) are trained only on non-restart transitions. We measure the quality of iterate z∈X×Yz∈ X× Y in original LP coordinates using the primal residual, dual residual, and duality-gap measure. Let c¯:=c−A⊤y c:=c-A y be the reduced cost and λ:=projΛ⁡(c¯)λ:=proj_ ( c ) its projection onto Λ . The primal and dual residuals and the duality-gap measure are defined as rp(z):=projY⁡(b−Ax),rd(z):=c¯−λ,rg(z):=c⊤x−(b⊤y+minx∈X⁡λ⊤x).r_p (z ):=proj_Y (b-Ax ), r_d (z ):= c-λ, r_g (z ):=c x- (b y+ _x∈ Xλ x ). (6) We denote the logarithms of their normalized magnitudes by p⁡(z)p (z ), d⁡(z)d (z ), and g⁡(z)g (z ) and define the merit function as Mrel(z):=max⁡p⁡(z),d⁡(z),g⁡(z).M^rel (z ):= \p (z ),d (z ),g (z ) \. (7) For the target tolerance ϵtol _tol, an episode ends successfully when the convergence criterion Mrel(zt+1)≤log⁡ϵtolM^rel (z_t+1 )≤ _tol is met. If the episode does not converge, it is truncated after TmaxT_ controlled iterations. With [a]+:=max⁡a,0[a]^+:= \a,0\, we use the following reward function to encourage reaching tolerance in fewer iterations, rt=Mrel(zt)−Mrel(zt+1)−χ[Mrel(zt+1)−logϵtol]+−1Tmax+[Mrel(zt+1)≤logϵtol]. r_t\!=\!M^rel\! (z_t )\!-\!M^rel\! (z_t+1 )\!-\!χ [M^rel\! (z_t+1 )\!-\! _tol ]^\!+\!\!\!-\! 1T_ \!+\!1\! [M^rel\! (z_t+1 )\!≤\! _tol ]\!. (8) Here, χ≥0χ≥ 0 weights the excess-error penalty [Mrel(zt+1)−log⁡ϵtol]+ [M^rel\! (z_t+1 )\!-\! _tol ]^\!+ and [⋅]1[·] is the indicator function. The difference in merit Mrel(zt)−Mrel(zt+1)M^rel\! (z_t )\!-\!M^rel\! (z_t+1 ) rewards a reduction in KKT error from one iteration to the next, even when the error remains above the target tolerance. The excess-error penalty encourages reducing errors earlier along the trajectory. Summed over an episode without discounting, the last two terms contribute 1−T/Tmax1-T/T_ for success after T iterations and −1-1 for truncation without success, favoring successful solves that reach tolerance in fewer iterations. 3.3 Policy Learning with Proximal Policy Optimization We design a groupwise proximal policy optimization (PPO) objective that reflects both conditional activation and the distinct control roles in GA-PDHG. We assign actions to five groups, =e,β,(α,ζ),η,ωG=\e,β,(α,ζ),η,ω\, with a Bernoulli policy for restart and Gaussian policies for the continuous controls. For trajectories collected under fixed policy parameters θold _old, let atja_t^j be the recorded action of group j∈j , and let πθj(atj|ot)π^j_θ(a_t^j|o_t) denote its probability mass or density under policy parameters θ. On a transition where group j is active, its likelihood ratio is ρtj:=πθj(atj|ot)πθoldj(atj|ot) _t^j:= π^j_θ\;(a_t^j|o_t)π^j_ _old(a_t^j|o_t). On restart transitions, we marginalize the unused acceleration controls in groups β and (α,ζ)(α,ζ) out of the likelihood and set ρtj≡1 _t^j≡ 1 for these inactive groups, so they contribute no policy gradient. A shared state-value critic Vϕ(ot)V^φ (o_t ) provides value estimates for generalized advantage estimation (GAE) (Schulman et al., 2016). Let Ψ~t _t denote the resulting advantage normalized over the collected transitions. All groups use this common advantage to learn from the return of their joint execution. The objective of our groupwise PPO is ⁡(θ)=∑j∈^t[min⁡(ρtjΨ~t,clip⁡(ρtj,1−εclip,1+εclip)Ψ~t)].J (θ )= _j E_t [ ( _t^j _t,\, ( _t^j,1- ^clip,1+ ^clip ) _t ) ]. (9) Here, ^t E_t averages over the collected transitions. εclip>0 ^clip>0 is the clipping parameter, which is applied to each group’s likelihood ratio ρtj _t^j. For the standard PPO (Schulman et al., 2017), clipping is applied to the joint action likelihood ratio ∏j∈ρtj _j _t^j. Under the joint clipping, changes in one group’s likelihood can suppress the gradient for the other group’s policy even when the the other group’s likelihood ratio is close to one. In particular, restart probability is often below 10%10\%, hence a small logit change has little effect on the non-restart likelihood. On the other hand β is sampled from k-dimensional Gaussian, which accumulates coordinate contributions in its log-likelihood ratio. Table 2: Clipping ratio (%). Problem class β e Transport 3.79 0.35 Auction 8.99 0.94 Random LPs 2.01 0.13 Table 2 shows that clipping of β happens more frequently than restart e, motivating separate clipping decisions. Under our groupwise PPO objective in Eq. (9), each group’s clipping depends only on its own ratio, so a change in one group’s ratio cannot by itself clip another group’s ratio. Similar discussions can be found in Fan et al. (2019). At inference, we fix θ, threshold the restart probability at 0.50.5, and transform the Gaussian means into continuous controls. GALLOP then executes the feedback loop in Figure 2 until the KKT tolerance is reached or the iteration limit is exhausted. The critic is used only during training. See Appendix F for details. Off-policy actor–critic methods, such as SAC (Haarnoja et al., 2018), can also be extended to hybrid action spaces (Delalleau et al., 2019). However, these methods require the hybrid structure to be represented in both the policy and the action-value critic, including the treatment of inactive continuous controls. Because environment interaction is inexpensive in our setting, the sample-efficiency gains from off-policy data reuse offer limited benefit. 4 Experiments 4.1 Experiment Setting Table 3: Density nnz⁡(A)/(mn)nnz (A )/(mn) and m/nm/n at each family’s largest training size. Family Density m/nm/n Random LP 1.28%1.28\% 0.50.5 Transport 0.625%0.625\% 0.01250.0125 Auction 0.313%0.313\% 0.250.25 Min-cost flow 0.0417%0.0417\% 0.1670.167 Robust inventory 0.326%0.326\% 0.3410.341 L∞L_∞ regression 100%100\% 26.726.7 We evaluate our method on six different LP families in Table 3; Random LP, Transport, Auction, Min-cost flow, Robust inventory management, and L∞L_∞ regression. We train one policy per family and evaluate on unseen LPs with varying sizes within the same family. Appendix D gives the problem details and size lists. We compare GALLOP with MPAX’s r2r^2HPDHG solver (Lu et al., 2024) and GA-PDHG (HT), a hand-tuned baseline that shares our solver implementation. GA-PDHG (HT) uses the reflected Halpern coefficients in Table 1, and ηt=0.998/∥A∥2 _t=0.998/ A _2, with restarting rules and PID controller weight updates from cuPDLPx (Lu et al., 2025). Because MPAX checks convergence every 100 iterations by default whereas our implementation checks every iteration, we use saved trajectories to identify the first iterate satisfying Mrel(zt)≤log⁡ϵtolM^rel (z_t )≤ _tol, with ϵtol=10−6 _tol=10^-6. We report iterations and matrix–vector products to this common tolerance, including products from MPAX’s rejected adaptive steps. Algorithm wall clock measures GPU execution of the main solve loop after compilation and excludes JIT compilation and preconditioning. 4.2 Convergence and Size Generalization Figure 3 plots the exponentiated merit exp⁡(Mrel(zt)) (M^rel(z_t)) against iteration count. GALLOP reaches the tolerance 10−610^-6 in fewer iterations than both baselines on all six plotted instances. Table 4 reports medians over four unseen LPs per family and size. (a) Random LP (b) Transport (c) Auction (d) Min-cost flow (e) Robust inventory (f) L∞L_∞ regression Figure 3: Convergence on one unseen LP per family at the largest training size. Labels mark the first iteration meeting tolerance 10−610^-6 (dashed line). The first row for each family corresponds to the largest training size (1×1×), and the following rows evaluate the same policy on larger LPs without retraining. At 1×1×, GALLOP requires fewer iterations and matrix–vector products for all six families. GALLOP generalizes to 15×15× for Random LP, 400×400× for Transport, 14×14× for Auction, 12×12× for Min-cost flow, 3×3× for Robust inventory, and 7×7× for L∞L_∞ regression, with gains in all three metrics at these sizes. Size multipliers are measured by the number of variables, except for Robust inventory (item count) and L∞L_∞ regression (constraint count at fixed n=9n=9). Table 4: Test results (four-LP medians). The “Size ratio” column reports the ratio of test size to each family’s largest training size. Size is measured in n except for robust inventory (item count N) and L∞L_∞ regression (m at fixed n=9n=9). Ours denotes GALLOP, HT denotes GA-PDHG (HT). Iterations Matrix–vector products Algorithm wall clock (ms) Family (m,n)(m,n) Size ratio Ours HT MPAX Ours HT MPAX Ours HT MPAX Random LP (10,000, 20,000) 1×1× 284 440 664 570 883 1,384 110 140 408 (30,000, 60,000) 3×3× 280 421 663 561 844 1,384 151 199 1,175 (150,000, 300,000) 15×15× 290 428 702 582 859 1,468 450 626 7,202 Transport (320, 25,600) 1×1× 164 633 548 329 1,268 1,150 39 107 87 (1,600, 640,000) 25×25× 153 533 536 308 1,068 1,128 102 286 308 (6,400, 10,240,000) 400×400× 200 570 524 403 1,142 1,105 1,998 5,222 7,334 Auction (6,400, 25,600) 1×1× 190 928 980 381 1,859 2,025 70 256 198 (32,000, 128,000) 5×5× 190 923 994 382 1,848 2,055 84 316 479 (89,600, 358,400) 14×14× 192 926 1,066 387 1,853 2,198 134 521 1,352 Min-cost flow (4,800, 28,800) 1×1× 164 302 440 329 606 928 49 68 76 (24,000, 144,000) 5×5× 234 356 444 469 715 939 80 94 109 (57,600, 345,600) 12×12× 242 370 482 485 741 1,010 115 142 147 Robust inventory (1,500, 4,401) 1×1× 1,905 4,266 3,875 3,812 8,535 7,872 507 915 432 (2,980, 8,801) 2×2× 2,311 7,768 5,348 4,624 15,537 10,854 570 1,494 661 (4,460, 13,201) 3×3× 2,874 25,798 6,121 5,749 51,599 12,391 858 5,978 895 L∞L_∞ regression (240, 9) 1×1× 150 586 834 303 1,175 1,734 24 54 82 (1,200, 9) 5×5× 236 673 760 475 1,348 1,583 41 78 69 (1,680, 9) 7×7× 358 752 782 718 1,505 1,626 62 93 78 Figure 4: Size generalization (four-LP medians). Shading marks the largest training size. Figure 5: Transport control ablations (four-LP medians, without retraining). Control group Replaced by Iterations — (GALLOP as trained) 165.0 Step size ηt _t 0.998/∥A∥20.998/ A _2 457.0 Primal weight ωt _t PID rule Eq. (65) 224.5 αtx,αty,ζtα^x_t,α^y_t, _t Halpern 209.0 History t β_t t= β_t=0 190.5 History t β_t Anderson 908.0 Restart flag ete_t three conditions Eq. (66) 181.0 All groups GA-PDHG (HT) 633.0 Figure 5 highlights the corresponding algorithm wall clock savings for random LP and Transport. At 400×400×, Transport takes 1,998 ms for GALLOP, compared with 5,222 ms for GA-PDHG (HT) and 7,334 ms for MPAX. Tables 10 and 11 provide results for all test sizes. 4.3 Ablation study We evaluate the contribution of each learned control group on Transport at the largest training size, m=320m=320 and n=25,600n=25,600. At evaluation time, we replace one control group with a hand-tuned rule and retain the other learned controls without retraining. We also test Anderson acceleration in place of the learned history coefficients. Table 5 reports the median number of iterations to reach ϵtol=10−6 _tol=10^-6 on all four LPs. Among substitutions with hand-tuned rules, replacing the learned step size has the largest effect, increasing the iterations from 165 to 457. Replacing any of the other four groups increases the median by 10%10\%–36%36\%. These results support the effectiveness of GA-PDHG and retaining all five learned control groups in the trained policy. Since the remaining controls are kept fixed rather than retrained, this experiment measures reliance of the learned joint policy on each control group, rather than the best performance achievable after removing that group and retraining. 4.4 Comparison on Public dataset Table 5: IP-S. Mean iteration counts. Solver Iterations GALLOP 186.7 MPAX 1,089.0 PDHG-Net + PDLP 422 PDLP 550.4 Table 6: Comparison results on IP-L. Zero-shot evaluation. Solver Iterations GALLOP 272.2 MPAX OOM PDLP 851.2 OOM: out of memory on RTX 4090. We evaluate on the public item-placement benchmark from Li et al. (2024). IP-S has (m,n)=(15,525,31,350) (m,n )= (15,525,31,350 ) and nnz⁡(A)=5,291,250nnz (A )=5,291,250, while IP-L has (m,n)=(91,575,266,450) (m,n )= (91,575,266,450 ) and nnz⁡(A)=94,826,250nnz (A )=94,826,250. We train GALLOP on 128 IP-S instances and evaluate the same policy on ten unseen instances of IP-S and IP-L. IP-L evaluation tests zero-shot size generalization without retraining. Tables 6 and 6 report mean iteration counts for each solvers. GALLOP and MPAX use tolerance 10−610^-6. PDLP uses the released PDHG-Net settings with absolute tolerance 10−610^-6 and relative tolerance 10−410^-4. This stopping criterion is looser than GALLOP’s. PDLP is the OR-Tools 9.8.3296 version run on our machine. The IP-S result of PDHG-Net is quoted from Table 2 in Li et al. (2024), which predicts a initialization for PDLP, hence the label PDHG-Net + PDLP. On IP-S, GALLOP requires 5.83×5.83× and 2.95×2.95× less number of iterations relative to MPAX and PDLP, respectively, and uses fewer iterations than the reported PDHG-Net + PDLP result. On IP-L, the reduction relative to PDLP is a factor of 3.13×3.13×. Its iteration count increases only 1.46×1.46× despite 8.50×8.50× more variables and 17.92×17.92× more nonzeros, demonstrating size generalization on this benchmark. 5 Related Work RL for discrete-continuous hybrid action spaces. P-DQN (Xiong et al., 2018) and MP-DQN (Bester et al., 2019) use action-value critics, unlike GALLOP’s state-value critic. H-PPO (Fan et al., 2019) also uses a state-value critic and separately clips actor objectives; DISC (Han & Sung, 2019) instead clips coordinatewise importance ratios. GALLOP specializes actor-wise clipping to five functional solver-control groups, marginalizing acceleration controls on restart transitions while retaining step-size and primal-weight updates. L2O for LP. Prior work learns simplex initial bases (Fan et al., 2023), imitates pivot experts (Liu et al., 2024), or emulates interior-point methods (Qian et al., 2024). These approaches target initialization or expert-derived decisions rather than jointly learning PDHG acceleration and restarting. More directly related, PDHG-Net (Li et al., 2024) unrolls PDHG to warm-start PDLP. GALLOP instead retains learned feedback throughout the solve, without handing subsequent iterations to an external solver. RL in L2O. RL learns simplex pivot-rule selection (Suriyanarayana et al., 2022) and presolve routines (Kuang et al., 2025), addressing discrete solver decisions rather than hybrid iterative control. RLQP (Ichnowski et al., 2021) adapts continuous penalty parameters of the alternating direction method of multipliers (ADMM), whereas TFPnP (Wei et al., 2020) combines continuous parameter tuning with discrete termination for inverse imaging. GALLOP instead jointly controls PDHG step sizes, primal weights, generalized acceleration, and restarts through a dimension-agnostic policy. Its discrete decisions reset acceleration state rather than terminate optimization. 6 Conclusion and Limitations We introduced GALLOP, which uses reinforcement learning to jointly control PDHG parameters, acceleration coefficients, and restart decisions. Its GA-PDHG update combines primal and dual extrapolation, history corrections, and an anchor correction in a common action space. A dimension-agnostic policy enables deployment on larger LPs without retraining. Across six LP families, GALLOP reduces iteration counts relative to MPAX and achieves algorithm wall-clock speedups of up to 16.0×16.0×. Policies trained on smaller LPs generalize to sizes up to 400×400× larger, with additional evidence of size generalization on the public item-placement benchmark. Limitations. Available GPU memory limited the maximum LP sizes evaluated in this study. Experiments with more GPU memory would allow us to test whether the observed size generalization extends to substantially larger LPs. References Amos (2023) Brandon Amos. Tutorial on amortized optimization. Foundations and Trends in Machine Learning, 16(5):592–732, 2023. ISSN 1935-8237. doi: 10.1561/2200000102. Anderson (1965) Donald G Anderson. Iterative procedures for nonlinear integral equations. Journal of the ACM (JACM), 12(4):547–560, 1965. Applegate et al. (2021) 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. Applegate et al. (2023) David Applegate, Oliver Hinder, Haihao Lu, and Miles Lubin. Faster first-order primal-dual methods for linear programming using restarts and sharpness. Mathematical Programming, 201(1):133–184, 2023. doi: 10.1007/s10107-022-01901-9. Applegate et al. (2026) David Applegate, Mateo Díaz, Oliver Hinder, Haihao Lu, Miles Lubin, Brendan O’Donoghue, and Warren Schudy. PDLP: A practical first-order method for large-scale linear programming. Mathematical Programming Computation, 2026. doi: 10.1007/s12532-026-00309-2. Bester et al. (2019) Craig J. Bester, Steven D. James, and George D. Konidaris. Multi-pass Q-networks for deep reinforcement learning with parameterised action spaces. arXiv preprint arXiv:1905.04388, 2019. URL https://arxiv.org/abs/1905.04388. Boţ & Nguyen (2023) Radu Ioan Boţ and Dang Khoa Nguyen. Fast Krasnosel’skiĭ–Mann algorithm with a convergence rate of the fixed point iteration of o⁡(1/k)o(1/k). SIAM Journal on Numerical Analysis, 61(6):2813–2843, 2023. doi: 10.1137/22M1504305. Chambolle & Pock (2011) 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(1):120–145, 2011. Chen et al. (2022) Tianlong Chen, Xiaohan Chen, Wuyang Chen, Howard Heaton, Jialin Liu, Zhangyang Wang, and Wotao Yin. Learning to optimize: A primer and a benchmark. Journal of Machine Learning Research, 23(189):1–59, 2022. Dantzig (1963) George B. Dantzig. Linear Programming and Extensions. Princeton University Press, 1963. Delalleau et al. (2019) Olivier Delalleau, Maxim Peter, Eloi Alonso, and Adrien Logut. Discrete and continuous action representation for practical RL in video games. arXiv preprint arXiv:1912.11077, 2019. URL https://arxiv.org/abs/1912.11077. Fan et al. (2023) Zhenan Fan, Xinglu Wang, Oleksandr Yakovenko, Abdullah Ali Sivas, Owen Ren, Yong Zhang, and Zirui Zhou. Smart initial basis selection for linear programs. In Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, p. 9650–9664. PMLR, 2023. URL https://proceedings.mlr.press/v202/fan23d.html. Fan et al. (2019) Zhou Fan, Rui Su, Weinan Zhang, and Yong Yu. Hybrid actor-critic reinforcement learning in parameterized action space. In Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, p. 2279–2285. International Joint Conferences on Artificial Intelligence Organization, 2019. doi: 10.24963/ijcai.2019/316. URL https://w.ijcai.org/proceedings/2019/316. Gregor & LeCun (2010) Karol Gregor and Yann LeCun. Learning fast approximations of sparse coding. In Proceedings of the 27th international conference on machine learning, p. 399–406, 2010. Haarnoja et al. (2018) Tuomas Haarnoja, Aurick Zhou, Pieter Abbeel, and Sergey Levine. Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor. In Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, p. 1861–1870. PMLR, 2018. URL https://proceedings.mlr.press/v80/haarnoja18b.html. Halpern (1967) Benjamin Halpern. Fixed points of nonexpanding maps. Bulletin of the American Mathematical Society, 73(6):957–961, 1967. Han & Sung (2019) Seungyul Han and Youngchul Sung. Dimension-wise importance sampling weight clipping for sample-efficient reinforcement learning. In Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, p. 2586–2595. PMLR, 2019. URL https://proceedings.mlr.press/v97/han19b.html. Hershey et al. (2014) John R. Hershey, Jonathan Le Roux, and Felix Weninger. Deep unfolding: Model-based inspiration of novel deep architectures, 2014. Hinder & Lubin (2020) Oliver Hinder and Miles Lubin. A generic adaptive restart scheme with applications to saddle point algorithms. arXiv preprint arXiv:2006.08484, 2020. Ichnowski et al. (2021) Jeffrey Ichnowski, Paras Jain, Bartolomeo Stellato, Goran Banjac, Michael Luo, Francesco Borrelli, Joseph Gonzalez, Ion Stoica, and Ken Goldberg. Accelerating quadratic optimization with reinforcement learning. In Advances in Neural Information Processing Systems, volume 34, p. 21043–21055. Curran Associates, Inc., 2021. Krasnosel’skii (1955) Mark Aleksandrovich Krasnosel’skii. Two remarks on the method of successive approximations. Uspekhi matematicheskikh nauk, 10(1):123–127, 1955. Kuang et al. (2025) Yufei Kuang, Xijun Li, Jie Wang, Fangzhou Zhu, Meng Lu, Zhihai Wang, Jia Zeng, Houqiang Li, Yongdong Zhang, and Feng Wu. Accelerate presolve in large-scale linear programming via reinforcement learning. IEEE Transactions on Pattern Analysis and Machine Intelligence, 47(8):6660–6672, 2025. doi: 10.1109/TPAMI.2025.3562286. Li et al. (2024) Bingheng Li, Linxin Yang, Yupeng Chen, Senmiao Wang, Haitao Mao, Qian Chen, Yao Ma, Akang Wang, Tian Ding, Jiliang Tang, and Ruoyu Sun. PDHG-unrolled learning-to-optimize method for large-scale linear programming. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, p. 29164–29180. PMLR, 2024. Liu et al. (2024) Tianhao Liu, Shanwen Pu, Dongdong Ge, and Yinyu Ye. Learning to pivot as a smart expert. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, p. 8073–8081, 2024. doi: 10.1609/aaai.v38i8.28646. Lu & Yang (2024) Haihao Lu and Jinwen Yang. Restarted Halpern PDHG for linear programming. arXiv preprint arXiv:2407.16144, 2024. Lu & Yang (2025) Haihao Lu and Jinwen Yang. cuPDLP.jl: A GPU implementation of restarted primal-dual hybrid gradient for linear programming in Julia. Operations Research, 73(6):3440–3452, 2025. doi: 10.1287/opre.2024.1069. Lu et al. (2023) 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. URL https://arxiv.org/abs/2312.14832. Lu et al. (2024) Haihao Lu, Zedong Peng, and Jinwen Yang. MPAX: mathematical programming in JAX. arXiv preprint arXiv:2412.09734, 2024. Lu et al. (2025) 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. Mann (1953) W Robert Mann. Mean value methods in iteration. Proceedings of the American Mathematical Society, 4(3):506–510, 1953. Martin & Belgioioso (2026) Andrea Martin and Giuseppe Belgioioso. Learning to accelerate Krasnosel’skii–Mann fixed-point iterations with guarantees. arXiv preprint arXiv:2601.07665, 2026. URL https://arxiv.org/abs/2601.07665. Metz et al. (2019) Luke Metz, Niru Maheswaranathan, Jeremy Nixon, Daniel Freeman, and Jascha Sohl-Dickstein. Understanding and correcting pathologies in the training of learned optimizers. In Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, p. 4556–4565. PMLR, 2019. Monga et al. (2021) Vishal Monga, Yuelong Li, and Yonina C. Eldar. Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing. IEEE Signal Processing Magazine, 38(2):18–44, 2021. ISSN 1053-5888, 1558-0792. doi: 10.1109/MSP.2020.3016905. Oshin et al. (2026a) Alex Oshin, Rahul Vodeb Ghosh, Augustinos D. Saravanos, and Evangelos A. Theodorou. Deep FlexQP: Accelerated nonlinear programming via deep unfolding. In The Fourteenth International Conference on Learning Representations, Rio de Janeiro, Brazil, 2026a. doi: 10.48550/arXiv.2512.01565. Oshin et al. (2026b) Alex Oshin, Rahul Vodeb Ghosh, and Evangelos A. Theodorou. Scalable deep unfolding of conic optimizers. In Proceedings of the 10th Conference on Robot Learning. PMLR, 2026b. doi: 10.48550/ARXIV.2606.13825. Pock & Chambolle (2011) Thomas Pock and Antonin Chambolle. Diagonal preconditioning for first order primal-dual algorithms in convex optimization. In 2011 International Conference on Computer Vision, p. 1762–1769. IEEE, 2011. Qian et al. (2024) Chendi Qian, Didier Chételat, and Christopher Morris. Exploring the power of graph neural networks in solving linear optimization problems. In Proceedings of the 27th International Conference on Artificial Intelligence and Statistics, volume 238 of Proceedings of Machine Learning Research, p. 1432–1440. PMLR, 2024. URL https://proceedings.mlr.press/v238/qian24a.html. Ruiz (2001) Daniel Ruiz. A scaling algorithm to equilibrate both rows and columns norms in matrices. Institut National Polytechnique de Toulouse, Toulouse,…, 2001. Sambharya & Stellato (2024) Rajiv Sambharya and Bartolomeo Stellato. Learning algorithm hyperparameters for fast parametric convex optimization, 2024. Sambharya & Stellato (2025) Rajiv Sambharya and Bartolomeo Stellato. Data-driven performance guarantees for classical and learned optimizers. Journal of Machine Learning Research, 26(171):1–49, 2025. Saravanos et al. (2025) Augustinos Saravanos, Hunter Kuperman, Alex Oshin, Arshiya Taj Abdul, Vincent Pacelli, and Evangelos Theodorou. Deep distributed optimization for large-scale quadratic programming. In International Conference on Learning Representations, volume 2025, p. 93580–93619, 2025. Schulman et al. (2016) John Schulman, Philipp Moritz, Sergey Levine, Michael I. Jordan, and Pieter Abbeel. High-dimensional continuous control using generalized advantage estimation. In International Conference on Learning Representations (ICLR), 2016. Schulman et al. (2017) John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347, 2017. Stellato et al. (2020) 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. Sucker et al. (2025) Michael Sucker, Jalal Fadili, and Peter Ochs. Learning-to-optimize with PAC-bayesian guarantees: Theoretical considerations and practical implementation. Journal of Machine Learning Research, 26(211):1–53, 2025. Suriyanarayana et al. (2022) Varun Suriyanarayana, Onur Tavaslıoğlu, Ankit B. Patel, and Andrew J. Schaefer. Reinforcement learning of simplex pivot rules: A proof of concept. Optimization Letters, 16(8):2513–2525, 2022. doi: 10.1007/s11590-022-01880-y. Venkataraman & Amos (2021) Shobha Venkataraman and Brandon Amos. Neural fixed-point acceleration for convex optimization. arXiv preprint arXiv:2107.10254, 2021. Wei et al. (2020) Kaixuan Wei, Angelica Aviles-Rivero, Jingwei Liang, Ying Fu, Carola-Bibiane Schönlieb, and Hua Huang. Tuning-free plug-and-play proximal algorithm for inverse imaging problems. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, p. 10158–10169. PMLR, 2020. URL https://proceedings.mlr.press/v119/wei20b.html. Wright (1997) Stephen J. Wright. Primal-Dual Interior-Point Methods. Society for Industrial and Applied Mathematics, 1997. Xiong et al. (2018) Jiechao Xiong, Qing Wang, Zhuoran Yang, Peng Sun, Lei Han, Yang Zheng, Haobo Fu, Tong Zhang, Ji Liu, and Han Liu. Parametrized deep Q-networks learning: Reinforcement learning with discrete-continuous hybrid action space. arXiv preprint arXiv:1810.06394, 2018. URL https://arxiv.org/abs/1810.06394. Zhou et al. (2025) Yingxin Zhou, Stefano Cipolla, and Phan Tu Vuong. Anderson Accelerated Primal-Dual Hybrid Gradient for solving LP. arXiv preprint arXiv:2508.08062, 2025. Appendix Appendix A Extended Related Work Learning to Optimize and Deep unfolding. Learning to optimize (L2O) leverages data-driven machine learning in order to develop and improve optimization methods. Key survey papers include (Chen et al., 2022) and (Amos, 2023). Deep unfolding is a model-based L2O approach that interprets each iteration of an optimizer as a layer of a neural network (Hershey et al., 2014). Learnable parameters are introduced and the optimizer is unrolled for a fixed number of iterations and trained using backpropagation-through-time, much like a recurrent neural network (Monga et al., 2021). The idea has seen great success in the signal and image processing communities after the seminal work by Gregor & LeCun (2010) on learned iterative shrinkage and soft-thresholding (LISTA). These approaches focus on learning iteration-varying parameters that can accelerate the optimizer, amortize the cost of tuning, and improve generalization on unseen instances (Sambharya & Stellato, 2024; Amos, 2023). Deep unfolding has been recently extended to general classes of convex optimization problems, such as linear programming (Li et al., 2024), quadratic programming (Saravanos et al., 2025; Oshin et al., 2026a), and conic programming (Oshin et al., 2026b). Another line of recent work has studied guarantees of learned optimizers, including PAC-Bayes generalization bounds (Sucker et al., 2025) and worst-case performance bounds (Sambharya & Stellato, 2025). However, deep unfolding often struggles to scale to high-dimensional problems in machine learning. This is primarily due to the intractable memory requirement of storing an entire trajectory of optimizer iterations in memory in order to backpropagate through them via supervised learning. These long-horizon unrolls are also numerically unstable, as they suffer from the vanishing/exploding gradient problem (Metz et al., 2019). Moreover, deep unfolding struggles to learn discrete decisions, such as when to restart in PDLP. Finally, supervised L2O requires optimal solutions, which are expensive to compute, especially for large-scale problems. These points make reinforcement learning a more suitable choice for the large-scale LPs that we target in this work. Appendix B Proofs from the Paper B.1 Fixed Points of PDHG and LP Optimality We first establish that solving the LP is equivalent to finding a fixed point of the PDHG operator ℋH. Fix ω,η>0ω,η>0 and write ℋ=ℋ⁡(⋅,ω,η)H=H (·;ω,η ). Let Z⋆Z denote the set of primal–dual optimal solutions of Eq. (1). For α∈[0,1]α∈[0,1], define ℛα(z):=(1+α)ℋ(z)−αzR_α (z ):= (1+α )H (z )-α z. Define the primal–dual KKT mapping ℱ:X×Y→2ℝn+mF:X× Y→ 2^R^n+m by ℱ⁡(z)=ℱ⁡(x,y):=(c−A⊤y+NX(x)Ax−b+NY(y)),F(z)=F(x,y):= pmatrixc-A y+N_X(x)\\ Ax-b+N_Y(y) pmatrix, (10) where NXN_X and NYN_Y denote the normal-cone mappings of X and Y, respectively. Lemma B.1 (KKT characterization). A primal–dual point z⋆=(x⋆,y⋆)z =(x ,y ) is optimal for Eq. (1) if and only if 0∈ℱ⁡(z⋆).0 (z ). (11) Equivalently, 0 0 ∈c−A⊤y⋆+NX(x⋆), ∈ c-A y +N_X(x ), (12) 0 0 ∈Ax⋆−b+NY(y⋆), ∈ Ax -b+N_Y(y ), (13) where Eq. (13) is equivalent to Px⋆ Px =q, =q, (14) Gx⋆ Gx ≥h, ≥ h, (15) yin,⋆ y^in, ≥0, ≥ 0, (16) yin,⋆⊙(Gx⋆−h) y^in, (Gx -h) =0. =0. (17) Proof. The saddle-point formulation in Eq. (2) is minx∈X⁡maxy∈Y⁡(c⊤x−y⊤Ax+b⊤y). _x∈ X _y∈ Y (c x-y Ax+b y ). (18) Its first-order optimality conditions are 0∈c−A⊤y⋆+NX(x⋆),0∈Ax⋆−b+NY(y⋆),0∈ c-A y +N_X(x ), 0∈ Ax -b+N_Y(y ), (19) which are necessary and sufficient because the problem is convex–concave and polyhedral. For the equality component, yeq∈ℝmey^eq ^m_e is unrestricted, and therefore the corresponding normal cone is zero. Hence Px⋆=qPx =q. For the inequality component, 0∈Gx⋆−h+Nℝ+mi(yin,⋆),0∈ Gx -h+N_R^m_i_+(y^in, ), (20) which is equivalent to primal feasibility, dual feasibility, and complementary slackness. ∎ Lemma B.2 (Equivalence between LP solutions and PDHG fixed points). The primal–dual optimal set coincides with the fixed-point set of the PDHG operator: Z⋆=Fix⁡(ℋ):=z∈X×Y:ℋ⁡(z)=z.Z =Fix(H):=\z∈ X× Y:H(z)=z\. (21) Furthermore, Fix⁡(ℛα)=Fix⁡(ℋ)=Z⋆.Fix(R_α)=Fix(H)=Z . (22) Proof. Let z=ℋ⁡(z^)=(x,y).z=H( z)=(x,y). (23) Eq. (4) gives x x =projX⁡[x^−ηω(c−A⊤y^)], =proj_X [ x- ηω(c-A y) ], (24) y y =projY⁡[y^+ωη(b−A⁡(2x−x^))]. =proj_Y [ y+ωη (b-A(2x- x) ) ]. (25) Using the projection characterization v=projC(a)⟺a−v∈NC(v),v=proj_C(a) a-v∈ N_C(v), (26) Eq. (24) and Eq. (25) imply ωη(x^−x)+A⊤y^−c ωη( x-x)+A y-c ∈NX(x), ∈ N_X(x), (27) 1ωη(y^−y)+b−A⁡(2x−x^) 1ωη( y-y)+b-A(2x- x) ∈NY(y). ∈ N_Y(y). (28) Suppose first that z^=ℋ⁡(z^) z=H( z). Then x=x^x= x and y=y^y= y, so A⊤y^−c∈NX(x^),b−Ax^∈NY(y^).A y-c∈ N_X( x), b-A x∈ N_Y( y). (29) Equivalently, 0∈c−A⊤y^+NX(x^),0∈Ax^−b+NY(y^).0∈ c-A y+N_X( x), 0∈ A x-b+N_Y( y). (30) By Lemma B.1, z^∈Z⋆ z∈ Z . Conversely, let z^∈Z⋆ z∈ Z . Then the KKT conditions imply A⊤y^−c∈NX(x^),b−Ax^∈NY(y^).A y-c∈ N_X( x), b-A x∈ N_Y( y). (31) The projection characterization gives x^=projX⁡[x^−ηω(c−A⊤y^)] x=proj_X [ x- ηω(c-A y) ] (32) and y^=projY⁡[y^+ωη(b−Ax^)]. y=proj_Y [ y+ωη(b-A x) ]. (33) Since the first equation gives x=x^x= x, the second equation agrees with Eq. (25). Hence ℋ⁡(z^)=z^H( z)= z. Finally, z^−ℛα(z^)=(1+α)(z^−ℋ⁡(z^)), z-R_α( z)=(1+α) ( z-H( z) ), (34) and therefore Fix⁡(ℛα)=Fix⁡(ℋ)Fix(R_α)=Fix(H). ∎ Lemma B.3 (Metric and nonexpansiveness). For fixed ω,η>0ω,η>0, define Q⁡(ω,η):=(ωηInA⊤A1ωηIm).Q (ω,η ):= pmatrix ωηI_n&A \\ A& 1ωηI_m pmatrix. (35) Then Q≻0Q 0 if and only if η∥A∥2<1η A _2<1. Under this condition, ℋ=ℋ⁡(⋅,ω,η)H=H (·;ω,η ) is firmly nonexpansive in the norm ∥v∥Q:=v⊤Qv v _Q:= v Qv, ‖ℋ⁡(z^)−ℋ⁡(z^′)‖Q2≤(ℋ⁡(z^)−ℋ⁡(z^′))⊤Q(z^−z^′). ( z )-H ( z ) _Q^2≤ (H ( z )-H ( z ) ) Q ( z- z ). (36) Consequently, it is nonexpansive in this norm, and the fixed-point residual Mfp(z^,ω,η):=‖ℋ⁡(z^,ω,η)−z^‖QM^fp ( z;ω,η ):= ( z;ω,η )- z _Q (37) vanishes if and only if z^∈Z⋆ z∈ Z . Proof. The Schur complement of the primal diagonal block in Eq. (35) gives Q≻0⟺Im−η2A⊤≻0⟺η∥A∥2<1.Q 0 I_m-η^2A 0 η A _2<1. (38) Fix parameters satisfying this condition. For z=ℋ⁡(z^,ω,η)z=H ( z;ω,η ), adding the two projection inequalities gives, for every z′∈X×Yz ∈ X× Y, (z^−z)⊤Q(z′−z) ( z-z ) Q (z -z ) ≤(c−A⊤y)⊤(x′−x) ≤ (c-A y ) (x -x ) (39) +(Ax−b)⊤(y′−y). + (Ax-b ) (y -y ). Apply this inequality to outputs z=ℋ⁡(z^)z=H ( z ) and z′=ℋ⁡(z^′)z =H ( z ), interchange the primed and unprimed variables, and add. The coupling terms cancel, yielding ∥z−z′∥Q2≤(z−z′)⊤Q(z^−z^′). z-z _Q^2≤ (z-z ) Q ( z- z ). (40) This proves firm nonexpansiveness (Lu & Yang, 2024). Applying Cauchy–Schwarz in the Q-inner product gives ∥z−z′∥Q≤∥z^−z^′∥Q z-z _Q≤ z- z _Q. Positive definiteness gives Mfp(z^,ω,η)=0M^fp ( z;ω,η )=0 if and only if ℋ⁡(z^)=z^H ( z )= z. Lemma B.2 then gives z^∈Z⋆ z∈ Z . ∎ B.2 Proof of Proposition 3.1 Proof of Proposition 3.1. By convergence, write z^t→z^¯ z_t→ z and zt→z¯z_t→ z. Then ‖zt−zt−1‖2⟶0,‖z^t+1−z^t‖2⟶0. z_t-z_t-1 _2 0, z_t+1- z_t _2 0. (41) Since t β_t is bounded, choose Bβ<∞B_β<∞ such that ∥t∥2≤Bβ β_t _2≤ B_β for all t. The fixed, finite history length gives, for sufficiently large t, ‖tt‖2 _t β_t _2 ≤Bβ(∑j=0k−1‖zt−j−zt−j−1‖22)1/2⟶0. ≤ B_β ( _j=0^k-1 z_t-j-z_t-j-1 _2^2 )^1/2 0. (42) For k=0k=0, the history correction and the sum are zero. At a non-restart step, rearranging Eq. (5a) gives ((1+αtx)(xt−x^t)(1+αty)(yt−y^t))=z^t+1−z^t−tt−ζt(z⋆−zt). pmatrix (1+ _t^x ) (x_t- x_t )\\ (1+ _t^y ) (y_t- y_t ) pmatrix= z_t+1- z_t-K_t β_t- _t (z_ -z_t ). (43) Set Cδ:=max⁡1,δ−1C_δ:= \1,δ^-1\. The lower bound 1+αtx,1+αty≥δ1+ _t^x,1+ _t^y≥δ on nonrestart steps, the triangle inequality, and the vanishing anchor correction assumed in Proposition 3.1 give ‖zt−z^t‖2≤Cδ(‖z^t+1−z^t‖2+‖tt‖2+‖ζt(z⋆−zt)‖2)⟶0. z_t- z_t _2≤ C_δ ( z_t+1- z_t _2+ _t β_t _2+ _t (z_ -z_t ) _2 ) 0. (44) At a restart, z^t+1=zt z_t+1=z_t, so the same bound holds because Cδ≥1C_δ≥ 1. Thus z^¯=z¯ z= z. Since ηt _t and ωt _t are uniformly positive and bounded, there exist constants satisfying 0<η¯≤ηt≤η¯<∞,0<ω¯≤ωt≤ω¯<∞.0< η≤ _t≤ η<∞, 0< ω≤ _t≤ ω<∞. (45) Hence a subsequence tj→∞t_j→∞ satisfies (ωtj,ηtj)⟶(ω∞,η∞)∈[ω¯,ω¯]×[η¯,η¯]. ( _t_j, _t_j ) ( _∞, _∞ )∈ [ ω, ω ]× [ η, η ]. (46) Euclidean projections onto the nonempty closed convex sets X and Y are continuous, so ℋH is jointly continuous in its input and positive parameters. Taking limits in Eq. (5b) along this subsequence gives z¯=limj→∞ztj+1=limj→∞ℋ⁡(z^tj+1,ωtj,ηtj)=ℋ⁡(z¯,ω∞,η∞). z= _j→∞z_t_j+1= _j→∞H ( z_t_j+1; _t_j, _t_j )=H ( z; _∞, _∞ ). (47) Thus z¯ z is a fixed point of ℋ⁡(⋅,ω∞,η∞)H (·; _∞, _∞ ). By Lemma B.2, z¯∈Z⋆ z∈ Z , completing the proof. ∎ B.3 Satisfying the assumptions in implementation We relate the deterministic inference rule in Appendix F.1 to the assumptions of Proposition 3.1. The trained network parameters are fixed during each solve. The implementation enforces the coefficient bounds and makes the anchor correction vanish whenever the output sequence converges. Convergence of both z^t\ z_t\ and zt\z_t\ remains a premise of the proposition. Bounded history coefficients and positive relaxation. Each continuous actor has a final tanh hidden layer followed by a linear output layer. For any scalar mean output, write μ⁡(o)=wμ⊤hμ(o)+bμ,hμ(o)∈[−1,1]64,|μ⁡(o)|≤∥wμ∥1+|bμ|.μ (o )=w_μ h_μ (o )+b_μ, h_μ (o )∈[-1,1]^64, |μ (o )|≤ w_μ _1+|b_μ|. (48) The right-hand side is a finite constant for the trained network, independent of the observation and iteration. At inference, t=tβ β_t= μ_t^β, so the fixed number k of history coefficients is uniformly bounded. The exponential transforms give αtx,αty>0 _t^x, _t^y>0, satisfying the proposition’s relaxation lower bound with δ=1δ=1, and a bounded positive anchor multiplier ζ~t=exp⁡(μtζ)≤Cζ ζ_t= ( _t^ζ )≤ C_ζ for some Cζ<∞C_ζ<∞. Step-size and primal-weight bounds. Let aA>0a_A>0 denote the fixed spectral-norm denominator, including the numerical safeguard, for the matrix supplied to PDHG. Let Lω>0L_ω>0 be the configured primal-weight range in Eq. (93). The sigmoid and tanh transforms give 0.1aA≤ηt≤5aA,ω0e−Lω≤ωt≤ω0eLω. 0.1a_A≤ _t≤ 5a_A, _0e^-L_ω≤ _t≤ _0e^L_ω. (49) Since aAa_A and ω0 _0 are positive and finite for each LP, both parameters are uniformly positive and bounded. The weight is set relative to ω0 _0 at every iteration, so these bounds do not accumulate over successive updates. Proposition 3.1 does not require ηt∥A∥2<1 _t A _2<1. Vanishing normalized anchor correction. The implemented coefficient is ζt=ζ~t‖zt−zt−1‖2‖z⋆−zt‖2,z⋆≠zt,0,z⋆=zt. _t= cases ζ_t z_t-z_t-1 _2 z_ -z_t _2,&z_ ≠ z_t,\\[6.0pt] 0,&z_ =z_t. cases (50) The zero case is handled explicitly before division. For either case, if zt→z¯z_t→ z, then ‖ζt(z⋆−zt)‖2≤Cζ‖zt−zt−1‖2⟶0. _t (z_ -z_t ) _2≤ C_ζ z_t-z_t-1 _2 0. (51) Thus the anchor correction vanishes as required by Proposition 3.1, independently of the restart pattern. Only the anchor correction must vanish; the effective coefficient ζt _t need not itself be bounded. Alternative Halpern-style parameterization. Using the epoch counter s from Section 3.1, an alternative parameterization is ζt=min⁡1,ζ~ts+1,0≤ζt≤min⁡1,Cζs+1. _t= \1, ζ_ts+1 \, 0≤ _t≤ \1, C_ζs+1 \. (52) The associated relaxation coefficients αtx=exp⁡(utαx)(1−ζt),αty=exp⁡(utαy)(1−ζt) _t^x= (u_t^α^x ) (1- _t ), _t^y= (u_t^α^y ) (1- _t ) (53) are nonnegative, satisfying the proposition’s relaxation lower bound with δ=1δ=1. Suppose zt→z¯z_t→ z. If only finitely many restarts occur, then z⋆z_ is eventually fixed and s→∞s→∞, so ‖ζt(z⋆−zt)‖2≤Cζs+1‖z⋆−zt‖2⟶0. _t (z_ -z_t ) _2≤ C_ζs+1 z_ -z_t _2 0. (54) If infinitely many restarts occur, the iteration at which the current anchor was set tends to infinity. Since each new anchor is an output iterate, z⋆→z¯z_ → z. Using 0≤ζt≤10≤ _t≤ 1 gives ‖ζt(z⋆−zt)‖2≤‖z⋆−zt‖2⟶0. _t (z_ -z_t ) _2≤ z_ -z_t _2 0. (55) Thus the anchor correction vanishes as required by Proposition 3.1 for either restart pattern. During training, Gaussian exploration has unbounded support, so the uniform bounds on t β_t and ζ~t ζ_t established here do not extend automatically to an infinite sampled trajectory. Proposition 3.1 applies to any trajectory satisfying its hypotheses, whether its controls are deterministic or random. B.4 Coefficient mappings for classical acceleration Fix ω,η>0ω,η>0 and write ℋ=ℋ⁡(⋅,ω,η)H=H (·;ω,η ). Within an uninterrupted run, write us:=z^tu_s:= z_t and vs:=zt=ℋ⁡(us)v_s:=z_t=H (u_s ) at local iteration s. The following mappings concern the acceleration recurrence with each method’s standard initialization, independently of the rollout initialization and restart branch in Algorithm 1. Halpern and reflected Halpern. Set u0=z⋆u_0=z_ and t=0 β_t=0. Substituting the first two rows of Table 1 into Eq. (5a) gives, respectively, us+1 u_s+1 =s+1s+2ℋ(us)+1s+2u0, = s+1s+2H (u_s )+ 1s+2u_0, (56) us+1 u_s+1 =s+1s+2(2ℋ(us)−us)+1s+2u0. = s+1s+2 (2H (u_s )-u_s )+ 1s+2u_0. These are Halpern iteration on ℋH and on its reflection 2ℋ−I2H-I, respectively (Halpern, 1967; Lu & Yang, 2024). Krasnosel’skiĭ–Mann. For constant 0<λ<20<λ<2, setting αtx=αty=λ−1 _t^x= _t^y=λ-1, t=0 β_t=0, and ζt=0 _t=0 gives us+1=(1−λ)us+λℋ(us).u_s+1= (1-λ )u_s+ (u_s ). (57) This is the KM relaxation of ℋH, with overrelaxation when λ>1λ>1 (Boţ & Nguyen, 2023). Fast KM. For ν>2ν>2, k≥1k≥ 1, and s≥1s≥ 1, the final row of Table 1 gives us+1 u_s+1 =vs+(ν2(s+ν)−1)(vs−us)+s+ν(vs−vs−1) =v_s+ ( ν2 (s+ν )-1 ) (v_s-u_s )+ ss+ν (v_s-v_s-1 ) (58) =(1−ν2(s+ν))us+ν2(s+ν)ℋ(us)+s+ν(ℋ⁡(us)−ℋ⁡(us−1)). = (1- ν2 (s+ν ) )u_s+ ν2 (s+ν )H (u_s )+ ss+ν (H (u_s )-H (u_s-1 ) ). This is equation (2.8) of Boţ & Nguyen (2023), with their damping parameter α replaced by ν and their step parameter set to one. Initialize u0,u1u_0,u_1 and store v0=ℋ⁡(u0)v_0=H (u_0 ) and v1=ℋ⁡(u1)v_1=H (u_1 ). The general Fast KM recurrence with a nonunit step also contains an input-difference term, so the table specifies the unit-step variant. Anderson coefficients. The vector tAA∈ℝk β_t^A ^k specifies how much of each past output difference to add to the current PDHG output ztz_t. Anderson acceleration chooses these coefficients by using past changes in the fixed-point residual to approximately cancel the current residual. For the PDHG input–output pairs (z^i,zi) ( z_i,z_i ), where zi=ℋ⁡(z^i)z_i=H ( z_i ), define dt:=zt−z^t,t:=[dt−dt−1⋯dt−k+1−dt−k].d_t:=z_t- z_t, _t:= bmatrixd_t-d_t-1&·s&d_t-k+1-d_t-k bmatrix. (59) The residual dtd_t vanishes exactly when z^t z_t is a fixed point of ℋH. The columns of tF_t record residual changes over the same history as the output differences in tK_t. For trial coefficients b∈ℝkb ^k, the residual model dt+tbd_t+F_tb combines the current residual with these past changes. Undamped type-I Anderson acceleration chooses b so that tbF_tb cancels dtd_t as closely as possible in the weighted norm ∥v∥Q2=v⊤Qv v _Q^2=v Qv. Here, Q≻0Q 0 is the metric in Eq. (35), with η∥A∥2<1η A _2<1. Using the available history columns and fixing all unavailable coefficients to zero, this gives tAA∈argminb∈ℝk⁡∥dt+tb∥Q2. β_t^A∈ *arg\,min_b ^k d_t+F_tb _Q^2. (60) The optimization variables are the history coefficients, and the objective is the squared norm of this residual model, not the exact residual at the next iterate. Apply the fitted coefficient of each residual difference to the matching output difference to form the next PDHG input, z^t+1=zt+ttAA. z_t+1=z_t+K_t β_t^A. (61) This is the undamped Anderson update in difference form (Anderson, 1965), recovered from Eq. (5a) by setting αtx=αty=ζt=0 _t^x= _t^y= _t=0. For a full window and b=tAAb= β_t^A, both the residual model and the update use the affine weights 1+b11+b_1, bj+1−bjb_j+1-b_j for 1≤j<k1≤ j<k, and −bk-b_k on the current and past residuals or outputs, respectively. These weights sum to one. Table 1 describes the algebraic GA-PDHG family. The exponential parameterization in Appendix F.1 restricts the learned extrapolation coefficients to positive values, so it does not implement the negative coefficients of underrelaxed KM or the displayed Fast KM specialization. Appendix C Details on PDHG-based LP solver enhancements This appendix gives the enhancement rules summarized in Section 2.2. Table 7 compares these enhancements across representative solvers. The four subsections follow the table columns. Table 7: Representative PDHG-based LP solver papers and their enhancement techniques. Paper / solver Step size Primal weight Restarting criteria Acceleration method PDLP (Applegate et al., 2021) Adaptive EMA Normalized duality gap Averaging raPDHG (Applegate et al., 2023) Constant Constant Normalized duality gap Averaging cuPDLP-C (Lu et al., 2023) Adaptive EMA KKT error Averaging cuPDLP.jl (Lu & Yang, 2025) Adaptive EMA KKT error Averaging rHPDHG / r2HPDHG (Lu & Yang, 2024) Adaptive EMA Fixed-point residual Halpern / Reflected Halpern MPAX (LP) (Lu et al., 2024) Adaptive EMA Fixed-point residual Reflected Halpern cuPDLPx (Lu et al., 2025) Constant PID Fixed-point residual Reflected Halpern A-PDHG (Zhou et al., 2025) Constant †EMA N/A Anderson EMA denotes exponential averaging in log space. PID denotes PID controller. N/A denotes no epoch restart. Paired acceleration entries follow the order of the variants. †A-PDHG uses periodic EMA and resets the Anderson history at weight updates. C.1 Step size raPDHG, cuPDLPx, and A-PDHG use a constant step size ηt≡η _t≡η. The condition η∥A∥2<1η A _2<1 ensures nonexpansiveness of the PDHG map in the induced norm (Appendix B.1). Choosing a constant step size from this condition can be conservative (Applegate et al., 2021), because it must hold for all possible updates, even when the current update permits a larger step. PDLP, cuPDLP.jl/C, MPAX, and rHPDHG/r2HPDHG instead choose the step size ηt _t adaptively from trial updates. Concretely, for a trial update zt+1=ℋ⁡(z^t+1,ωt,ηt)z_t+1=H ( z_t+1; _t, _t ), the step size acceptance test verifies ηt≤∥zt+1−z^t+1∥ωt22|(yt+1−y^t+1)⊤A(xt+1−x^t+1)|,where∥z∥ω2:=ω∥x∥22+ω−1∥y∥22. _t≤ z_t+1- z_t+1 _ _t^22 | (y_t+1- y_t+1 ) A (x_t+1- x_t+1 ) |, z _ω^2:=ω x _2^2+ω^-1 y _2^2. (62) If the test fails, it reduces ηt _t and recomputes zt+1z_t+1 until the condition holds. Checking only the current trial can permit steps beyond the sufficient spectral-norm bound. Rejected trials incur extra matrix–vector products. Although the test is motivated by the nonexpansiveness proof of the PDHG map, to the best of our knowledge, no convergence guarantee has been established for the adaptive step-size rule used in practical PDLP-style solvers (Applegate et al., 2021; Lu & Yang, 2024). C.2 Primal weight The primal weight adjusts the relative primal and dual step sizes to balance the scales of the primal and dual variables. The convergence of standard PDHG is proven for a constant ωt≡ω>0 _t≡ω>0 and a constant step size η<1/∥A∥2η<1/ A _2 (Chambolle & Pock, 2011; Applegate et al., 2023). Adaptive primal-weight updates are not covered by this fixed-weight convergence guarantee. In the hand-designed adaptive schemes reviewed here, the primal weight is updated at the end of each block of iterations, with block boundaries set by restarts or fixed update intervals. It remains constant within each block. Superscript j indexes blocks, subscript t indexes individual iterations, and ωj>0ω^j>0 is the weight in block j. Let (xj,yj) (x^j,y^j ) denote the starting point of block j. The displacement norms between successive block starting points are dxj:=‖xj+1−xj‖2,dyj:=‖yj+1−yj‖2.d_x^j:= x^j+1-x^j _2, d_y^j:= y^j+1-y^j _2. (63) For dxj,dyjd_x^j,d_y^j, define the log-balance error bj:=log⁡(ωjdxjdyj).b^j:= ( ω^jd^j_xd^j_y ). (64) PDLP, cuPDLP-C/jl, rHPDHG/r2HPDHG, MPAX, and A-PDHG use exponential moving averaging (EMA) in log space, whereas cuPDLPx uses proportional–integral–derivative (PID) control. At the end of block j, the weight for the next block is computed as log⁡ωj+1=log⁡ωj−θωbj,EMA,Kpbj+Ki∑ℓ=0jbℓ+Kd(bj−bj−1),PID, ω^j+1= ω^j- cases _ωb^j,&EMA,\\ K_pb^j+K_i _ =0^jb +K_d (b^j-b^j-1 ),&PID, cases (65) where θω∈[0,1] _ω∈[0,1] is the smoothing coefficient, and Kp,Ki,KdK_p,K_i,K_d are PID gains. Both rules implement feedback control to balance the weighted displacements ωjdxj ω^jd_x^j and dyj/ωjd_y^j/ ω^j, corresponding to bj=0b^j=0. EMA changes log⁡ωj ω^j by −θωbj- _ωb^j, so the correction is proportional to the current error. PID control uses past errors through its integral and derivative terms (Lu et al., 2025). C.3 Restarting Restarting reinitializes PDHG from a selected point z⋆z_ and resets its past averaging or acceleration state. When restarting is combined with averaging or Halpern acceleration, suitable parameter choices and restart rules exploit LP sharpness to achieve accelerated linear convergence (Applegate et al., 2023; Lu & Yang, 2024). Restart timing affects both the theoretical iteration complexity and the practical convergence of restarted PDHG. A fixed-frequency scheme restarts after a prescribed number of iterations. In practice, the restarted solvers in Table 7 use adaptive restart heuristics based on the conditions Mt≤βsufM⋆ M_t≤β^sufM_ (sufficient progress), (sufficient progress), (66a) Mt≤βnecM⋆andMt>Mt−1 M_t≤β^necM_ M_t>M_t-1 (loss of progress), (loss of progress), (66b) s+1≥βartt s+1≥β^artt (long iterations after restart), (long iterations after restart), (66c) Here, Mt:=M⁡(zt)M_t:=M(z_t) is a restart metric evaluated at iterate ztz_t, and M⋆:=M⁡(z⋆)M_ :=M(z_ ) is its value at the restart anchor z⋆z_ . In Eq. (66c), s is the zero-based index within the current epoch, so s+1s+1 counts its completed PDHG steps, while t≥1t≥ 1 counts all completed PDHG steps. These conditions are checked at a prescribed frequency using thresholds satisfying 0<βsuf<βnec<10<β^suf<β^nec<1 and βart∈(0,1)β^art∈ (0,1 ). The metric M, manually tuned constant thresholds (βsufβ^suf, βnecβ^nec, βartβ^art), and the condition-checking intervals differ across solvers, as do initialization rules such as more frequent checks during the first few iterations. C.4 Acceleration An acceleration step constructs the next PDHG input or a restart candidate from the current PDHG update and past iterates. The solvers in Table 7 use Averaging, Halpern, or Anderson accelerations. Averaging. Averaging uses the step-size weighted average z¯t=∑i=t−stηi−1zi∑i=t−stηi−1 z_t= _i=t-s^t _i-1z_i _i=t-s^t _i-1 as the restart candidate. When primal–dual iterates oscillate around a solution, averaging can cancel deviations in opposite directions. Restarted averaging with a fixed step size η and primal weight ω exploits LP sharpness to achieve accelerated linear convergence under the prescribed restart rules based on the normalized duality gap or KKT error (Applegate et al., 2023; Lu & Yang, 2025). Halpern acceleration. Halpern acceleration combines the current PDHG input z^t z_t, output ztz_t, and fixed restart anchor z⋆z_ to form the next PDHG input z^t+1=zt+α⁡(s+1)s+2(zt−z^t)+1s+2(z⋆−zt). z_t+1=z_t+ α(s+1)s+2 (z_t- z_t )+ 1s+2 (z_ -z_t ). (67) Here, α is the reflection coefficient, with α=0α=0 recovering standard Halpern acceleration and α=1α=1 yielding reflected Halpern acceleration (Lu & Yang, 2024). The anchor coefficient 1/(s+2)1/ (s+2 ) decreases as s increases and returns to 1/21/2 at each restart, when s is reset to zero. For feasible and bounded LPs, both variants achieve an O⁡(1/(s+1))O (1/ (s+1 ) ) fixed-point residual bound with fixed η,ω>0η,ω>0 satisfying η∥A∥2<1η A _2<1. When combined with the restart schemes analyzed by Lu & Yang (2024), they achieve accelerated linear convergence. Anderson acceleration. Anderson acceleration corrects the current PDHG output using recent output differences (Zhou et al., 2025). Using the k most recent PDHG output differences, Anderson acceleration takes the form z^t+1=zt+ttAA,t:=[zt−zt−1⋯zt−k+1−zt−k]. z_t+1=z_t+K_t β_t^A, _t:= bmatrixz_t-z_t-1&·s&z_t-k+1-z_t-k bmatrix. (68) The coefficients tAA∈ℝk β_t^A ^k are computed from recent PDHG input–output pairs (z^i,zi)( z_i,z_i). A small least-squares problem selects coefficients that combine past residual changes to approximately cancel the current residual zt−z^tz_t- z_t. The same coefficients are then applied to the output differences to construct the next PDHG input. The resulting history correction has a quasi-Newton interpretation and can reduce error components that decay slowly under basic PDHG. These acceleration methods form linear combinations of the current and past iterates. Averaging uses step-size weights, Halpern acceleration uses coefficients prescribed by convergence analysis, and Anderson acceleration computes coefficients by solving a least-squares problem. Appendix D Further Details on Problem Classes We describe the transportation, auction, minimum-cost flow, robust inventory management, and L∞L_∞ polynomial regression generators used in Section 4.1. The transportation, auction, and minimum-cost flow generators construct primal–dual solutions satisfying the KKT conditions, so their optimal values are known without an auxiliary solve. The inventory generator solves each generated LP, while the regression generator also refines the resulting active set, as detailed below. Nonzero counts refer to the stacked constraint matrix A defined after Eq. (1), with variable bounds stored separately in X. Here nnz⁡(A)nnz(A) counts the nonzero entries, and matrix density is nnz⁡(A)/(mn)nnz (A )/(mn), the complement of the zero fraction. Symbols for problem-specific data in the following subsections are local to each LP class. D.1 Transport Balanced transportation allocates mass from s sources to t destinations at minimum cost. In this subsection only, s,ts,t are problem dimensions. Let a∈ℝ+sa _+^s and d∈ℝ+td _+^t be the supply and demand vectors with ∑i=1sai=∑j=1tdj _i=1^sa_i= _j=1^td_j, and let cij∈ℝc_ij be the unit cost from source i to destination j. For transported amounts xij≥0x_ij≥ 0, the LP is minxij≥0 _x_ij≥ 0 ∑i=1s∑j=1tcijxij _i=1^s _j=1^tc_ijx_ij (69) subject to to ∑j=1txij=ai,i=1,…,s, _j=1^tx_ij=a_i, i=1,…,s, ∑i=1sxij=dj,j=1,…,t. _i=1^sx_ij=d_j, j=1,…,t. The generator obtains the marginals from a feasible plan supported on a spanning tree of the complete bipartite graph, optionally with additional support arcs. After vectorization, the LP has n=stn=st variables, me=m=s+tm_e=m=s+t equality constraints, and no inequality rows. Each column of A=PA=P contains two unit entries, one for its source and one for its destination, giving nnz⁡(A)=2st=2n,nnz⁡(A)mn=2s+t.nnz (A )=2st=2n, nnz (A )mn= 2s+t. (70) All stst transport arcs remain in the LP, so changing the optimal plan’s support or introducing cost ties does not change matrix sparsity. One mass-balance equality is redundant and is retained. D.2 Auction We use a capacitated LP relaxation of combinatorial auction winner determination, where bids request bundles of items. For m items and n bids, let B∈0,1m×nB∈\0,1\^m× n have Bij=1B_ij=1 when bid j includes item i, let v∈ℝnv ^n contain bid values, and let d∈ℝ+md _+^m contain item capacities. The variable xj∈[0,1]x_j∈[0,1] is the accepted fraction of bid j, and the LP is maxx∈ℝnv⊤xsubject toBx≤d,0≤x≤. _x ^n v x to Bx≤ d, 0≤ x 1. (71) In Eq. (1), this corresponds to c=−vc=-v, G=−BG=-B, and h=−dh=-d, with me=0m_e=0 and mi=m_i=m. The generated capacities are not restricted to one. Each bundle contains κ∈1,…,mκ∈\1,…,m\ distinct items sampled without replacement, so nnz⁡(A)=κn,nnz⁡(A)mn=κm.nnz (A )=κ n, nnz (A )mn= κm. (72) Bundle size determines column sparsity, while the item-popularity parameter changes row degrees and overlap among bids without changing the nonzero count. The finite bounds allow both full and partial bid acceptance in the relaxation. D.3 Min-cost Flow Capacitated minimum-cost flow routes a commodity through a directed network while satisfying node balances and arc capacities. Let V and E be the node and arc sets, with m=|V|m=|V| and n=|E|n=|E|, and define the incidence matrix B∈ℝm×nB ^m× n to have −1-1 at the tail and +1+1 at the head of each arc. For arc costs c∈ℝnc ^n, capacities u∈ℝ>0nu _>0^n, and net inflows d∈ℝmd ^m satisfying ⊤d=01 d=0, the flow vector x∈ℝnx ^n solves minx∈ℝnc⊤xsubject toBx=d,0≤x≤u. _x ^n c x to Bx=d, 0≤ x≤ u. (73) The generator randomly orients a spanning tree and adds distinct directed arcs without self-loops, ensuring that the underlying undirected graph is connected. Here A=P=BA=P=B, me=m_e=m, and mi=0m_i=0, and each column contains exactly two nonzeros, giving nnz⁡(A)=2n,nnz⁡(A)mn=2m.nnz (A )=2n, nnz (A )mn= 2m. (74) Adding arcs increases row degrees and the total nonzero count, while column sparsity remains fixed. The incidence matrix has rank m−1m-1, and the redundant conservation row is retained. Capacities affect which box constraints are active without changing matrix sparsity. D.4 Robust Inventory Management We consider N items over T periods with uncertain demand and shared purchasing and storage limits. For item i=1,…,Ni=1,…,N and period t=1,…,Tt=1,…,T, let xt(i)x_t^(i), ut(i)u_t^(i), and bt(i)b_t^(i) denote inventory, purchases, and a nonnegative shortage buffer. The initial inventory x1(i)x_1^(i) is given. The nonnegative parameters ht(i)h_t^(i), ct(i)c_t^(i), pt(i)p_t^(i), st(i)s_t^(i), and v(i)v^(i) are the holding cost, purchase cost, selling price, shortage penalty, and unit volume. Let C and V be the per-period capital budget and storage capacity. Demand lies in a box, and purchases follow an itemwise causal affine policy, =∏i=1N∏t=1T[d¯t(i),d¯t(i)],ut(i)=u^t(i)+∑τ<tkt,τ(i)dτ(i).D= _i=1^N _t=1^T [ d_t^(i), d_t^(i) ], u_t^(i)= u_t^(i)+ _τ<tk_t,τ^(i)d_τ^(i). (75) The nominal purchases u^t(i) u_t^(i), feedback gains kt,τ(i)k_t,τ^(i), and buffers bt(i)b_t^(i) are chosen before demand is observed. Only purchases adapt to past demand, with no feedback between items. The underlying robust problem is minu^,k,b≥0⁡maxd∈ _ u,k,\,b≥ 0\ _d ∑i=1N∑t=1T[ht(i)xt(i)+ct(i)ut(i)−pt(i)dt(i)+(pt(i)+st(i))bt(i)] _i=1^N _t=1^T [h_t^(i)x_t^(i)+c_t^(i)u_t^(i)-p_t^(i)d_t^(i)+ (p_t^(i)+s_t^(i) )b_t^(i) ] (76) subject to to xt+1(i)=xt(i)+ut(i)−dt(i)+bt(i),i=1,…,N,t=1,…,T, x_t+1^(i)=x_t^(i)+u_t^(i)-d_t^(i)+b_t^(i), i=1,…,N, t=1,…,T, xt(i)≥0,ut(i)≥0,∑iv(i)xt(i)≤V,∑ict(i)ut(i)≤C,∀d∈. x_t^(i)≥ 0, u_t^(i)≥ 0, _iv^(i)x_t^(i)≤ V, _ic_t^(i)u_t^(i)≤ C, ∀ d . The last line applies for t=1,…,Tt=1,…,T and, for the itemwise inequalities, every i. Holding costs are charged on xtx_t, the inventory entering period t. The generated problem has no terminal constraint on xT+1x_T+1. The buffer btb_t is a deterministic feasibility slack, rather than the realized positive-part shortage, and has no demand-dependent upper bound. LP formulation. Define the demand-independent inventory at(i)a_t^(i) and demand coefficients βt,τ(i) _t,τ^(i) by at(i) a_t^(i) =x1(i)+∑r<t(u^r(i)+br(i)), =x_1^(i)+ _r<t ( u_r^(i)+b_r^(i) ), βt,τ(i) _t,τ^(i) =∑r=τ+1t−1kr,τ(i)−1,τ<t, = _r=τ+1^t-1k_r,τ^(i)-1, τ<t, (77) xt(i) x_t^(i) =at(i)+∑τ<tβt,τ(i)dτ(i). =a_t^(i)+ _τ<t _t,τ^(i)d_τ^(i). Empty sums are zero. The total cost is Γ0+∑i,τγτ(i)dτ(i) _0+ _i,τ _τ^(i)d_τ^(i), where Γ0 _0 =∑i=1N∑t=1T[ht(i)at(i)+ct(i)u^t(i)+(pt(i)+st(i))bt(i)], = _i=1^N _t=1^T [h_t^(i)a_t^(i)+c_t^(i) u_t^(i)+ (p_t^(i)+s_t^(i) )b_t^(i) ], (78) γτ(i) _τ^(i) =∑t=τ+1T[ht(i)βt,τ(i)+ct(i)kt,τ(i)]−pτ(i). = _t=τ+1^T [h_t^(i) _t,τ^(i)+c_t^(i)k_t,τ^(i) ]-p_τ^(i). For each coefficient q∈kt,τ(i),βt,τ(i),γτ(i)q∈\k_t,τ^(i), _t,τ^(i), _τ^(i)\, introduce q+,q−≥0q^+,q^-≥ 0 with q=q+−q−q=q^+-q^-. For complementary splits satisfying q+q−=0q^+q^-=0, the minimum and maximum of qdτ(i)qd_τ^(i) over its interval are, respectively, Lτ(i)(q±)=d¯τ(i)q+−d¯τ(i)q−,Uτ(i)(q±)=d¯τ(i)q+−d¯τ(i)q−.L_τ^(i) (q^± )= d_τ^(i)q^+- d_τ^(i)q^-, U_τ^(i) (q^± )= d_τ^(i)q^+- d_τ^(i)q^-. (79) Let ℛR be the sum of all split variables, ℛ=∑i=1N[∑τ<t(kt,τ(i),++kt,τ(i),−+βt,τ(i),++βt,τ(i),−)+∑τ=1T(γτ(i),++γτ(i),−)].R= _i=1^N [ _τ<t (k_t,τ^(i),++k_t,τ^(i),-+ _t,τ^(i),++ _t,τ^(i),- )+ _τ=1^T ( _τ^(i),++ _τ^(i),- ) ]. (80) With a free epigraph variable η, the generated LP is minη,u^,b,k±,β±,γ± _η, u,b,k^±,β^±,γ^± η+εℛ η+ (81) subject to to Γ0+∑i=1N∑τ=1TUτ(i)(γτ(i),±)≤η, _0+ _i=1^N _τ=1^TU_τ^(i) ( _τ^(i),± )≤η, ∑i=1Nct(i)[u^t(i)+∑τ<tUτ(i)(kt,τ(i),±)]≤C, _i=1^Nc_t^(i) [ u_t^(i)+ _τ<tU_τ^(i) (k_t,τ^(i),± ) ]≤ C, t=1,…,T, t=1,…,T, ∑i=1Nv(i)[at(i)+∑τ<tUτ(i)(βt,τ(i),±)]≤V, _i=1^Nv^(i) [a_t^(i)+ _τ<tU_τ^(i) ( _t,τ^(i),± ) ]≤ V, t=2,…,T, t=2,…,T, u^t(i)+∑τ<tLτ(i)(kt,τ(i),±)≥0, u_t^(i)+ _τ<tL_τ^(i) (k_t,τ^(i),± )≥ 0, i=1,…,N,t=1,…,T, i=1,…,N, t=1,…,T, at(i)+∑τ<tLτ(i)(βt,τ(i),±)≥0, a_t^(i)+ _τ<tL_τ^(i) ( _t,τ^(i),± )≥ 0, i=1,…,N,t=2,…,T, i=1,…,N, t=2,…,T, b,k±,β±,γ±≥0. b,\ k^±,\ β^±,\ γ^±≥ 0. Here ata_t, Γ0 _0, k, β, and γ are substituted using their definitions above. The definitions of β and γ impose linear equalities on the split variables. The nominal purchases u u are unrestricted. Removing the common part of any split pair preserves these equalities and weakens every inequality, so complementary splits suffice without explicit complementarity constraints. Thus ε=0 =0 gives an exact LP reformulation of Eq. (76). The benchmark uses ε=10−6 =10^-6 and includes this finite regularization in its objective. Generated instances. We use T=10T=10, x1(i)=0x_1^(i)=0, v(i)=1v^(i)=1, ct(i)=0.5c_t^(i)=0.5, ht(i)=0.3h_t^(i)=0.3, st(i)=0.1s_t^(i)=0.1, C=0.25NC=0.25N, and V=0.155NV=0.155N. Selling prices are sampled independently from ⁡[1,2]U[1,2]. For μt=0.25[1+cos⁡(2πt/7)] _t=0.25 [1+ (2π t/7 ) ], demand bounds are d¯t(i)=max⁡0,μt−0.075 d_t^(i)= \0, _t-0.075\ and d¯t(i)=μt+0.075 d_t^(i)= _t+0.075. The generator checks the constant initial-inventory constraints once and omits their rows. With variable bounds stored separately, the LP has n=1+2NT(T+1),me=NT(T+1)2,m=2T+N⁡[T⁡(T+5)2−1].n=1+2NT (T+1 ), m_e= NT (T+1 )2, m=2T+N [ T (T+5 )2-1 ]. (82) At the stated settings, nnz⁡(A)=1075N+1nnz (A )=1075N+1. Each instance is solved with HiGHS to obtain a primal–dual reference solution. D.5 L∞L_∞ Polynomial Regression L∞L_∞ polynomial regression minimizes the largest absolute residual over a finite set of observations. This criterion is also called minimax or Chebyshev regression and differs from least-squares regression, which minimizes the sum of squared residuals. Let (ti,yi)i=1N \ (t_i,y_i ) \_i=1^N be observations with ti∈[−1,1]t_i∈[-1,1], and let ϕjj=1p \ _j \_j=1^p be a basis for polynomials of degree at most p−1p-1. Write y=(y1,…,yN)y= (y_1,…,y_N ) and define the design matrix Φ∈ℝN×p ^N× p by Φij=ϕj(ti) _ij= _j (t_i ). The coefficient vector β=(β1,…,βp)∈ℝpβ= ( _1,…, _p ) ^p solves minβ∈ℝp⁡‖Φβ−y‖∞=minβ∈ℝp⁡maxi=1,…,N|∑j=1pβjϕj(ti)−yi|. _β ^p β-y _∞= _β ^p _i=1,…,N | _j=1^p _j _j (t_i )-y_i |. (83) Introducing an epigraph variable τ≥0τ≥ 0 gives the LP, where ∈ℝN1 ^N denotes the all-ones vector. minβ∈ℝp,τ≥0 _β ^p,\,τ≥ 0 τ τ (84) subject to to −τ≤Φβ−y≤τ. - 1≤ β-y≤ 1. For the convention Gx≥hGx≥ h in Eq. (1), set x=(β,τ)x= (β,τ ) and G=[Φ−Φ],h=[y−y],c=[p1].G= bmatrix &1\\ - &1 bmatrix, h= bmatrixy\\ -y bmatrix, c= bmatrix0_p\\ 1 bmatrix. (85) The coefficients β are free variables, whereas only τ has the lower bound zero. There are no equality constraints or finite upper bounds. Thus the generated LP has n=p+1n=p+1 variables and mi=m=2Nm_i=m=2N inequality rows. Because P is empty, the stacked matrix is A=GA=G. For continuously sampled nodes, this matrix is dense almost surely, with nnz⁡(A)=2N(p+1)=mnnnz (A )=2N (p+1 )=mn. Let TkT_k be the degree-k Chebyshev polynomial of the first kind, defined by Tk(cos⁡θ)=cos⁡(kθ)T_k ( θ )= (kθ ). For an excess degree r≥1r≥ 1 and noise scale σ≥0σ≥ 0, the generator samples a target function in the Chebyshev basis and adds Gaussian observation noise. yi=∑k=0p+r−1akTk(ti)+σϵi,ak,ϵi∼i.i.d.(0,1),y_i= _k=0^p+r-1a_kT_k (t_i )+σ _i,a_k, _i i.i.d. N (0,1 ), (86) We use r=3r=3 and σ=0.1σ=0.1. The nodes are sampled uniformly on [−1,1][-1,1] by default. The generator also supports the Chebyshev nodes ti=cos⁡((2i−1)π/(2N))t_i= ( (2i-1 )π/ (2N ) ) and center-clustered nodes ti=ui|ui|t_i=u_i |u_i | for ui∼Unif⁡[−1,1]u_i [-1,1 ]. The fitted basis is either the monomial basis ϕj(t)=tj−1 _j (t )=t^j-1 or the Chebyshev basis ϕj(t)=Tj−1(t) _j (t )=T_j-1 (t ). Both bases span the same degree-(p−1)(p-1) polynomial space, so paired instances with the same nodes and observations have the same optimal residual τ⋆τ in exact arithmetic. Their numerical conditioning can differ substantially. The monomial design is a Vandermonde matrix and can become ill-conditioned as the degree increases, while the Chebyshev representation is typically better conditioned on [−1,1][-1,1]. The basis choice therefore changes the LP representation without changing the regression function class. Unlike the preceding planted generators, we obtain the reference primal–dual solution by solving Eq. (84) with the HiGHS dual-simplex method and refining the active constraints by least squares. D.6 Training set construction We generate a fixed training set for each family using the sizes and instance counts in Tables 8 and 9. Each instance is generated from its bucket seed plus its within-bucket index, and the same LPs are reused across PPO iterations. The training sets contain 224 transport, 448 auction, 224 min-cost flow, 256 robust inventory, 512 regression, and 384 random LP instances. We train one policy per family, select one checkpoint by held-out evaluation, and compare it on the same four unseen LPs per test size as both baselines. Tables 8 and 9 list the training, validation, and test sizes, whose instances use disjoint seed ranges. Training counts denote LPs per PPO rollout, held-out counts denote LPs per evaluation, and compared counts denote LPs per solver. Size multipliers are relative to the largest training size and measured in n, except for inventory (item count N) and regression (m at fixed n=9n=9). For inventory, m=20+74Nm=20+74N and n=1+220Nn=1+220N at a fixed horizon of 10 periods. Table 8: Training, validation, and test sizes for transport, auction, and min-cost flow. Family Multiplier m n nnz⁡(A)nnz(A) LPs per split Train / Held-out Compared Transport 0.125 × 120 3,200 6,400 64 / – – 0.25 × 160 6,400 12,800 64 / 16 4 0.5 × 240 12,800 25,600 64 / – – 1 × 320 25,600 51,200 32 / 16 4 2.5 × 506 64,000 128,000 – / 8 4 5 × 720 128,000 256,000 – / 8 4 10 × 1,012 256,000 512,000 – / 4 4 12 × 1,112 307,200 614,400 – / – 4 14 × 1,200 358,400 716,800 – / – 4 16 × 1,280 409,600 819,200 – / – 4 25 × 1,600 640,000 1,280,000 – / – 4 36 × 1,920 921,600 1,843,200 – / – 4 64 × 2,560 1,638,400 3,276,800 – / – 4 100 × 3,200 2,560,000 5,120,000 – / – 4 400 × 6,400 10,240,000 20,480,000 – / – 4 Auction 0.125 × 800 3,200 64,000 128 / 16 – 0.25 × 1,600 6,400 128,000 128 / 16 – 0.5 × 3,200 12,800 256,000 128 / 16 – 1 × 6,400 25,600 512,000 64 / 16 4 2.5 × 16,000 64,000 1,280,000 – / – 4 4 × 25,600 102,400 2,048,000 – / 8 – 5 × 32,000 128,000 2,560,000 – / – 4 10 × 64,000 256,000 5,120,000 – / – 4 12 × 76,800 307,200 6,144,000 – / – 4 14 × 89,600 358,400 7,168,000 – / – 4 Min-cost flow 0.125 × 600 3,600 7,200 64 / – – 0.25 × 1,200 7,200 14,400 64 / 16 – 0.5 × 2,400 14,400 28,800 64 / – – 1 × 4,800 28,800 57,600 32 / 16 4 2.5 × 12,000 72,000 144,000 – / 8 4 5 × 24,000 144,000 288,000 – / 8 4 10 × 48,000 288,000 576,000 – / 4 4 12 × 57,600 345,600 691,200 – / – 4 Table 9: Training, validation, and test sizes for inventory, L∞L_∞ regression, and random LP. Family Multiplier m n nnz⁡(A)nnz(A) LPs per split Train / Held-out Compared Robust inventory 0.25 × 390 1,101 5,376 64 / – 4 0.5 × 760 2,201 10,751 128 / 4 4 1 × 1,500 4,401 21,501 64 / 4 4 1.5 × 2,240 6,601 32,251 – / 4 4 2 × 2,980 8,801 43,001 – / 2 4 3 × 4,460 13,201 64,501 – / 2 4 5 × 7,420 22,001 107,501 – / 2 4 8 × 11,860 35,201 172,001 – / – 4 10 × 14,820 44,001 215,001 – / – 4 L∞L_∞ regression 0.5 × 120 9 1,080 128 / 16 4 0.67 × 160 9 1,440 128 / 16 4 0.83 × 200 9 1,800 128 / 16 4 1 × 240 9 2,160 128 / 16 4 2 × 480 9 4,320 – / 16 4 5 × 1,200 9 10,800 – / 16 4 6 × 1,440 9 12,960 – / – 4 7 × 1,680 9 15,120 – / – 4 Random LP 0.2 × 2,000 4,000 512,000 64 / – – 0.4 × 4,000 8,000 1,024,000 128 / – 4 0.8 × 8,000 16,000 2,048,000 128 / 16 – 1 × 10,000 20,000 2,560,000 64 / – 4 2 × 20,000 40,000 5,120,000 – / 4 4 3 × 30,000 60,000 7,680,000 – / – 4 4 × 40,000 80,000 10,240,000 – / – 4 5 × 50,000 100,000 12,800,000 – / 4 – 6 × 60,000 120,000 15,360,000 – / – 4 10 × 100,000 200,000 25,600,000 – / 2 4 15 × 150,000 300,000 38,400,000 – / – 4 Appendix E Full Results on Six Problem Classes Tables 10 and 11 report the complete size-generalization results for all six LP families. Each entry is the median over four unseen LPs at the indicated size, evaluated by all three solvers. Iterations and matrix–vector products are counted to the first iterate satisfying the common relative KKT tolerance ϵtol=10−6 _tol=10^-6. Algorithm wall clock measures execution after compilation and excludes JIT compilation and preconditioning, as in Section 4.1. The dagger marks the largest training size, and a dash denotes failure to reach the tolerance within the iteration limit. Table 10: Full test results for transport, auction, and min-cost flow (four-LP medians). Iterations to ϵtol _tol Matrix–vector products Algorithm wall clock (ms) Family (m,n)(m,n) Ours HT MPAX Ours HT MPAX Ours HT MPAX Transport (160, 6,400) 146 488 550 294 979 1,158 46 119 76 (320, 25,600)† 164 633 548 329 1,268 1,150 39 107 87 (506, 64,000) 155 674 540 312 1,349 1,134 46 143 96 (720, 128,000) 150 628 526 302 1,257 1,108 48 146 122 (1,012, 256,000) 140 554 518 282 1,111 1,090 52 154 156 (1,112, 307,200) 158 628 528 317 1,259 1,110 57 173 170 (1,200, 358,400) 152 621 527 307 1,244 1,108 61 192 188 (1,280, 409,600) 148 712 524 299 1,427 1,101 66 247 203 (1,600, 640,000) 153 533 536 308 1,068 1,128 102 286 308 (1,920, 921,600) 152 594 538 305 1,190 1,130 135 447 445 (2,560, 1,638,400) 149 676 527 300 1,354 1,110 248 965 830 (3,200, 2,560,000) 149 608 532 300 1,217 1,120 385 1,399 1,445 (6,400, 10,240,000) 200 570 524 403 1,142 1,105 1,998 5,222 7,334 Auction (6,400, 25,600)† 190 928 980 381 1,859 2,025 70 256 198 (16,000, 64,000) 190 924 1,108 382 1,851 2,284 81 293 324 (32,000, 128,000) 190 923 994 382 1,848 2,055 84 316 479 (64,000, 256,000) 192 918 1,044 387 1,838 2,158 114 433 895 (76,800, 307,200) 192 924 1,046 387 1,849 2,162 118 453 1,143 (89,600, 358,400) 192 926 1,066 387 1,853 2,198 134 521 1,352 Min-cost flow (4,800, 28,800)† 164 302 440 329 606 928 49 68 76 (12,000, 72,000) 180 319 448 363 640 944 63 85 88 (24,000, 144,000) 234 356 444 469 715 939 80 94 109 (48,000, 288,000) 224 374 446 451 749 938 91 120 139 (57,600, 345,600) 242 370 482 485 741 1,010 115 142 147 Ours: GALLOP; HT: GA-PDHG (HT). Bold: best per metric. †: largest training size. A dash denotes failure to reach ϵtol=10−6 _tol=10^-6 within the iteration limit. Table 11: Full test results for inventory, L∞L_∞ regression, and random LP (four-LP medians). Iterations to ϵtol _tol Matrix–vector products Algorithm wall clock (ms) Family (m,n)(m,n) Ours HT MPAX Ours HT MPAX Ours HT MPAX Robust inventory (390, 1,101) 1,426 2,892 3,774 2,853 5,787 7,682 259 397 309 (760, 2,201) 1,250 4,242 3,072 2,503 8,485 6,256 260 659 482 (1,500, 4,401)† 1,905 4,266 3,875 3,812 8,535 7,872 507 915 432 (2,240, 6,601) 2,140 10,294 5,888 4,282 20,591 11,934 635 2,529 591 (2,980, 8,801) 2,311 7,768 5,348 4,624 15,537 10,854 570 1,494 661 (4,460, 13,201) 2,874 25,798 6,121 5,749 51,599 12,391 858 5,978 895 (7,420, 22,001) 3,833 12,103 6,167 7,668 24,208 12,495 1,164 3,034 975 (11,860, 35,201) 6,308 13,499 8,588 12,618 27,000 17,356 1,920 3,209 1,505 (14,820, 44,001) 7,832 – 11,581 15,666 – 23,370 2,239 – 1,834 L∞L_∞ regression (120, 9) 190 465 702 383 932 1,466 28 42 69 (160, 9) 169 554 907 340 1,110 1,876 25 51 81 (200, 9) 127 533 740 256 1,068 1,544 20 49 72 (240, 9)† 150 586 834 303 1,175 1,734 24 54 82 (480, 9) 270 1,060 786 543 2,121 1,637 44 113 72 (1,200, 9) 236 673 760 475 1,348 1,583 41 78 69 (1,440, 9) 308 741 764 619 1,484 1,590 54 92 76 (1,680, 9) 358 752 782 718 1,505 1,626 62 93 78 Random LP (4,000, 8,000) 264 424 706 531 850 1,468 87 109 221 (10,000, 20,000)† 284 440 664 570 883 1,384 110 140 408 (20,000, 40,000) 270 416 650 541 834 1,359 125 164 773 (30,000, 60,000) 280 421 663 561 844 1,384 151 199 1,175 (40,000, 80,000) 264 429 666 530 860 1,388 169 244 1,612 (60,000, 120,000) 304 426 774 609 854 1,609 243 310 2,443 (100,000, 200,000) 314 455 711 629 912 1,481 363 490 5,283 (150,000, 300,000) 290 428 702 582 859 1,468 450 626 7,202 Ours: GALLOP; HT: GA-PDHG (HT). Bold: best per metric. †: largest training size. A dash denotes failure to reach ϵtol=10−6 _tol=10^-6 within the iteration limit. Appendix F Implementation Details F.1 Network architecture and action distributions The GALLOP policy uses separate tanh MLPs for the restart decision, history coefficients, relaxation and anchor coefficients, step size, and primal weight. A separate network estimates the value function. Each network has two hidden layers of width 64. The restart head produces pte∈(0,1)p_t^e∈(0,1) through a sigmoid and samples et∼Bernoulli⁡(pte).e_t (p_t^e). (87) Conditional on oto_t, the continuous heads sample independent Gaussian latents with the means μtαx,μtαy,tβ,μtζ,μtη,μtω _t^α^x, _t^α^y, μ_t^β, _t^ζ, _t^η, _t^ω output by their respective networks; all but tβ∈ℝk μ_t^β ^k are scalars. Their standard deviations σα,σβ,σζ,ση,σω>0 _α, _β, _ζ, _η, _ω>0 are learned globally and do not depend on oto_t. Their logarithms are optimized with the corresponding actor and clipped to group-specific lower bounds and a common upper bound when the distributions are evaluated. Appendix G.1 gives the initialization values and bounds used for each family. The primal and dual relaxation latents share σα _α, and all history coordinates share σβ _β. We restrict the learned coefficients to αtx,αty>0 _t^x, _t^y>0 to focus the policy search on extrapolation beyond the PDHG output, motivated by potential acceleration rather than damping through underrelaxation. Exponential transformations enforce positive relaxation coefficients and a positive anchor multiplier, utαx u_t^α^x ∼⁡(μtαx,σα2), ( _t^α^x, _α^2 ), αtx _t^x =exp⁡(utαx), = (u_t^α^x ), (88) utαy u_t^α^y ∼⁡(μtαy,σα2), ( _t^α^y, _α^2 ), αty _t^y =exp⁡(utαy), = (u_t^α^y ), (89) utζ u_t^ζ ∼⁡(μtζ,σζ2), ( _t^ζ, _ζ^2 ), ζ~t ζ_t =exp⁡(utζ). = (u_t^ζ ). (90) The effective anchor coefficient ζt _t is computed by Eq. (50), with ζt=0 _t=0 when z⋆=ztz_ =z_t. During training, the history coefficients are sampled from an unbounded Gaussian distribution, t∼⁡(tβ,σβ2Ik). β_t ( μ_t^β, _β^2I_k ). (91) At every iteration, including restarts, the step-size sample is transformed as utη∼⁡(μtη,ση2),ηt=0.1+4.9sig⁡(utη)‖A‖2.u_t^η ( _t^η, _η^2 ), _t= 0.1+4.9 (u_t^η )\|A\|_2. (92) Here sig⁡(u):=1/(1+e−u) (u):=1/(1+e^-u) is the sigmoid. The denominator estimates the spectral norm of the matrix supplied to PDHG. Sparse instances use a converged spectral-norm calculation during both training and evaluation; dense training uses 100 power iterations, while dense evaluation uses a converged calculation. The initialization step uses 0.998/‖A‖20.998/\|A\|_2. The learned step size can exceed 1/‖A‖21/\|A\|_2 and therefore does not enforce the fixed-step nonexpansiveness condition. At every iteration, including non-restart transitions, the primal-weight sample sets an absolute target relative to ω0 _0, utω∼⁡(μtω,σω2),ωt=ω0exp⁡(Lωtanh⁡(utω)).u_t^ω ( _t^ω, _ω^2 ), _t= _0 (L_ω (u_t^ω ) ). (93) Here Lω=2L_ω=2 for auction and Lω=4.6L_ω=4.6 for the other five families. This weight is applied in the next PDHG map regardless of ete_t; it does not multiply the previous weight. The rollout records the raw action at=(et,utαx,utαy,t,utζ,utη,utω).a_t=(e_t,u_t^α^x,u_t^α^y, β_t,u_t^ζ,u_t^η,u_t^ω). Let atacta_t^act consist of ete_t, utηu_t^η, utωu_t^ω, and the acceleration groups enabled when et=0e_t=0. Marginalizing the unused draws gives the density of this effective action, πθ,act(atact∣ot)= _θ,act (a_t^act o_t )= (pte)et(1−pte)1−et(utη,μtη,ση2)(utω,μtω,σω2) (p_t^e)^e_t(1-p_t^e)^1-e_tN (u_t^η; _t^η, _η^2 )N (u_t^ω; _t^ω, _ω^2 ) (94) ×[⁡(utαx,μtαx,σα2)(utαy,μtαy,σα2)×⁡(t,tβ,σβ2Ik)(utζ,μtζ,σζ2)]1−et. × [ aligned &N (u_t^α^x; _t^α^x, _α^2 )N (u_t^α^y; _t^α^y, _α^2 )\\ &×N ( β_t; μ_t^β, _β^2I_k )N (u_t^ζ; _t^ζ, _ζ^2 ) aligned ]^1-e_t. The step-size and primal-weight factors are active on every transition, and the acceleration factors only when et=0e_t=0. Inactive factors contribute zero to the log-density. PPO evaluates the density on the raw latents, so the output transformations require no Jacobian correction in the likelihood ratio. Deterministic inference. At inference, et=[pte≥0.5]e_t=1[p_t^e≥ 0.5], and every Gaussian latent is replaced by its mean before applying the same transformations. The observation map and action transforms are unchanged, and the trained parameters remain fixed throughout the solve. When et=0e_t=0, the transformed coefficients form the accelerated input in Eq. (5a); when et=1e_t=1, the input is ztz_t and the history is reset after the PDHG step. The step size and primal weight are applied on both branches at every iteration. The critic and Gaussian standard deviations are needed for training, but do not determine the inference actions. F.2 GALLOP observations For t≥1t≥ 1, the observation oto_t is evaluated before ata_t using the carried weight ωt−1 _t-1. The initial weight ω0 _0 is used for the initialization step that produces z1z_1. It is ‖c‖2/‖b‖2\|c\|_2/\|b\|_2 in solver coordinates when both norms are positive, and 11 otherwise. The features below follow the convergence state, restart context, and local geometry described in Section 3.2. Convergence state. The primal residual, dual residual, and duality-gap diagnostic in original LP coordinates are defined in Eq. (6). The sign cone is Λ:=∏i=1nΛi := _i=1^n _i, where Λi=ℝ _i=R if both bounds are finite, ℝ≥0R_≥ 0 if only lil_i is finite, ℝ≤0R_≤ 0 if only uiu_i is finite, and 0\0\ if neither is finite. With elementwise [a]±:=max⁡±a,0[a]^±:= \± a,0\, the box minimum in rg(z)r_g (z ) equals l⊤[c¯]+−u⊤[c¯]−l [ c]^+-u [ c]^-, with infinite-bound terms omitted. When y∈Yy∈ Y and rd(z)=0r_d (z )=0, rg(z)r_g (z ) is the primal–dual objective gap. Otherwise it is a finite diagnostic used together with the feasibility residuals. The corresponding log normalized metrics are p⁡(z):=log⁡(max⁡‖rp(z)‖∞,‖[l−x]+‖∞,‖[x−u]+‖∞1+max⁡‖b‖∞,‖Ax‖∞).p(z):= ( \\|r_p(z)\|_∞,\|[l-x]^+\|_∞,\|[x-u]^+\|_∞ \1+ \\|b\|_∞,\|Ax\|_∞ \ ). (95) d⁡(z):=log⁡(max⁡‖rd(z)‖∞,‖[−yin]+‖∞1+max⁡‖c‖∞,‖A⊤y‖∞).d(z):= ( \\|r_d(z)\|_∞,\|[-y^in]^+\|_∞ \1+ \\|c\|_∞,\|A y\|_∞ \ ). (96) g⁡(z):=log⁡(|rg(z)|1+max⁡|c⊤x|,|c⊤x−rg(z)|).g(z):= ( |r_g(z)|1+ \|c x|,|c x-r_g(z)| \ ). (97) All three metrics are evaluated in the original LP coordinates, with log arguments floored at 10−2010^-20. The bound and dual-sign violation terms vanish on X×YX× Y. Thus these formulas also apply to the initial anchor z⋆=0z_ =0, which may lie outside X. For the controller, fixed-point progress uses the executed direction zt−z^tz_t- z_t in solver coordinates, measured in a fixed metric, Mtfp:=‖zt−z^t‖Q⁡(1,0.998/‖A‖2).M_t^fp:=\|z_t- z_t\|_Q (1,0.998/\|A\|_2 ). (98) The measurement metric uses the scaled A when preconditioning is enabled. Its spectral norm is computed by dense SVD for dense training batches and by a converged calculation for sparse batches and evaluation, independently of the dense training step-size estimate. For t≥1t≥ 1, it measures the step executed with (ωt−1,ηt−1) ( _t-1, _t-1 ). The observation includes progress relative to the initial step, the restart baseline M⋆M_ , and the preceding step, log⁡max⁡Mtfp,10−20max⁡M1fp,10−20,log⁡max⁡Mtfp,10−20max⁡M⋆,10−20,log⁡max⁡Mtfp,10−20max⁡Mt−1fp,10−20. \M_t^fp,10^-20\ \M_1^fp,10^-20\, \M_t^fp,10^-20\ \M_ ,10^-20\, \M_t^fp,10^-20\ \M_t-1^fp,10^-20\. (99) It also includes the one-step changes in the KKT metrics, p⁡(zt)−p⁡(zt−1),d⁡(zt)−d⁡(zt−1),g⁡(zt)−g⁡(zt−1).p(z_t)-p(z_t-1), d(z_t)-d(z_t-1), g(z_t)-g(z_t-1). (100) The one-step features are evaluated only for t≥2t≥ 2 and are set to zero at t=1t=1. The initial restart baseline is M⋆=M1fpM_ =M_1^fp. Restart context. The policy observes the log restart age and global iteration count, log⁡(s+1) (s+1) and log⁡t t, and the KKT progress from the anchor, p⁡(zt)−p⁡(z⋆),d⁡(zt)−d⁡(z⋆),g⁡(zt)−g⁡(z⋆).p(z_t)-p(z_ ), d(z_t)-d(z_ ), g(z_t)-g(z_ ). (101) At initialization, the anchor metrics are evaluated at z⋆=0z_ =0. After a restart step, z⋆←zt+1z_ ← z_t+1 and M⋆←Mt+1fpM_ ← M_t+1^fp, so the anchor-relative KKT differences are zero. Primal–dual balance, the relative weight, and its one-step change are measured by bte:=log⁡(ωt−1‖xt−x⋆‖2+10−20‖yt−y⋆‖2+10−20),log⁡(ωt−1ω0),log⁡(ωt−1ωt−2).b_t^e:= ( _t-1\|x_t-x_ \|_2+10^-20\|y_t-y_ \|_2+10^-20 ), ( _t-1 _0 ), ( _t-1 _t-2 ). (102) The weight-change feature is set to zero at t=1t=1. The displacements in these features use the solver coordinates. KKT residuals evaluated at z^t z_t are not supplied to the policy. Local geometry. Partition the history directions into primal and dual blocks, t=(txty),tx∈ℝn×k,ty∈ℝm×k.K_t= pmatrixK_t^x\\ K_t^y pmatrix, _t^x ^n× k, _t^y ^m× k. (103) For z=(x⊤,y⊤)⊤z= (x ,y ) , define the weighted norm ∥z∥ω2:=ω∥x∥22+ω−1∥y∥22. z _ω^2:=ω x _2^2+ω^-1 y _2^2. (104) Using this norm and ϵgeom=10−8 _geom=10^-8, the normalized history Gram matrix and its alignment with the current direction are Γt _t :=ωt−1(tx)⊤tx+ωt−1−1(ty)⊤ty‖zt−z^t‖ωt−12+ϵgeom2∈ℝk×k, := _t-1(K_t^x) K_t^x+ _t-1^-1(K_t^y) K_t^y\|z_t- z_t\|_ _t-1^2+ _geom^2 ^k× k, (105) cthist c_t^hist :=ωt−1(tx)⊤(xt−x^t)+ωt−1−1(ty)⊤(yt−y^t)‖zt−z^t‖ωt−12+ϵgeom2∈ℝk. := _t-1(K_t^x) (x_t- x_t )+ _t-1^-1(K_t^y) (y_t- y_t )\|z_t- z_t\|_ _t-1^2+ _geom^2 ^k. The anchor features measure its relative magnitude and weighted alignment with the current direction, ρtanc _t^anc :=log⁡‖z⋆−zt‖ωt−1+ϵgeom‖zt−z^t‖ωt−1+ϵgeom, := \|z_ -z_t\|_ _t-1+ _geom\|z_t- z_t\|_ _t-1+ _geom, (106) ctanc c_t^anc :=ωt−1(x⋆−xt)⊤(xt−x^t)+ωt−1−1(y⋆−yt)⊤(yt−y^t)(‖z⋆−zt‖ωt−1+ϵgeom)(‖zt−z^t‖ωt−1+ϵgeom). := _t-1(x_ -x_t) (x_t- x_t )+ _t-1^-1(y_ -y_t) (y_t- y_t ) (\|z_ -z_t\|_ _t-1+ _geom ) (\|z_t- z_t\|_ _t-1+ _geom ). The history quantities undergo an elementwise signed logarithm before entering the policy, vech⁡(sign⁡(Γt)log⁡(1+|Γt|)),sign⁡(cthist)log⁡(1+|cthist|), (sign( _t) (1+| _t| ) ), (c_t^hist) (1+|c_t^hist| ), (107) ρtanc,ctanc. _t^anc, c_t^anc. Here vechvech stacks the upper-triangular entries, including the diagonal, in row-major order. The sign, absolute value, logarithm, and their products act elementwise. The observation has 19+k+k⁡(k+1)/219+k+k(k+1)/2 entries, giving 33 inputs for k=4k=4. These inputs are passed directly to the networks without running normalization or input clipping. F.3 Policy optimization At each PPO iteration, the parameters (θold,ϕold)( _old, _old) are frozen while one episode is collected per LP instance. For an episode of T transitions, generalized advantage estimation (Schulman et al., 2016) uses discount factor γ=1γ=1 and GAE coefficient κ=0.99κ=0.99, δt _t :=rt+γ⁡(1−successt)Vϕold(ot+1)−Vϕold(ot), :=r_t+γ (1-success_t )V _old(o_t+1)-V _old(o_t), (108) Ψ^t _t :=∑j=0T−t(γκ)jδt+j,1≤t≤T. := _j=0^T-t(γκ)^j\, _t+j, 1≤ t≤ T. Here successt:=[Mrel(zt+1)≤logϵtol]success_t:=1 [M^rel(z_t+1)≤ _tol ]. The sum ends with the episode. A successful episode has terminal value zero; truncation at TmaxT_ retains Vϕold(oT+1)V _old (o_T+1 ) in the final temporal-difference residual. The critic target and the advantage normalized over the full rollout buffer are R^t:=Ψ^t+Vϕold(ot),Ψ~t:=Ψ^t−mean⁡(Ψ^)std⁡(Ψ^)+10−8. R_t:= _t+V _old(o_t), _t:= _t-mean( )std( )+10^-8. (109) Here mean⁡(Ψ^)mean( ) and std⁡(Ψ^)std( ) are the empirical mean and population standard deviation of all valid advantages in ℬB. Writing ^t E_t for the empirical average over that buffer, the critic loss is ℒV(ϕ):=^t[(Vϕ(ot)−R^t)2].L_V(φ):= E_t\! [ (V^φ(o_t)- R_t )^2 ]. (110) Groupwise PPO. Let πθ,e _θ,e and πθ,β _θ,β denote the restart and history densities, πθ,αζ _θ,αζ the product of the relaxation and anchor Gaussian densities, and πθ,η,πθ,ω _θ,η, _θ,ω the step-size and primal-weight latent densities. Each density below is evaluated at its recorded component of ata_t, conditional on the recorded oto_t; these arguments are suppressed. For the action groups =e,β,αζ,η,ωG=\e,β,αζ,η,ω\, the masked log-densities are ℓt:=(log⁡πθ,e,(1−et)log⁡πθ,β,(1−et)log⁡πθ,αζ,log⁡πθ,η,log⁡πθ,ω)⊤. _t:= ( _θ,e,(1-e_t) _θ,β,(1-e_t) _θ,αζ, _θ,η, _θ,ω ) . (111) For j∈j , let ℓtj _t^j be the corresponding entry of ℓt _t and ℓt,oldj _t,old^j its value under θold _old. Each group uses its own likelihood ratio and clipped surrogate, with clipping parameter εclip=0.2 ^clip=0.2, ρtj _t^j :=exp⁡(ℓtj−ℓt,oldj), := ( _t^j- _t,old^j ), (112) j(θ) _j(θ) :=^t[min⁡(ρtjΨ~t,clip⁡(ρtj,1−εclip,1+εclip)Ψ~t)]. := E_t [ ( _t^j _t, ( _t^j,1- ^clip,1+ ^clip ) _t ) ]. Here clip⁡(v,a,b):=min⁡b,max⁡a,v (v,a,b):= \b, \a,v\\. The five surrogates are summed; the advantages and recorded old log-densities are held fixed during optimization. On restart rows, the acceleration groups have ρtj=1 _t^j=1 and contribute no policy gradient. The restart, step-size, and primal-weight surrogates use all valid transitions in the common minibatches and their full-trajectory, globally normalized advantages Ψ~t _t. The primal-weight actor is therefore also updated when the rollout contains no restarts. Each actor and the critic has its own gradient clipping and Adam state. The critic minimizes cVℒVc_VL_V, with cVc_V specified in Appendix G.1. Auxiliary restart supervision. Before the action is sampled, the conventional sufficient-progress, deterioration, and epoch-length criteria provide a binary label ytpriory_t^prior. The label is one if Mtfp≤0.2M⋆M_t^fp≤ 0.2M_ , if both Mtfp≤0.6M⋆M_t^fp≤ 0.6M_ and Mtfp>Mt−1fpM_t^fp>M_t-1^fp, or if s≥0.2ts≥ 0.2t. At t=1t=1, the preceding residual in this test is taken to equal M1fpM_1^fp. It is zero otherwise and is always set to zero when s=0s=0. These tests use the policy trajectory and do not execute a restart. The auxiliary loss and the resulting actor loss are ℒprior _prior :=−^t[ytpriorlog⁡pte+(1−ytprior)log⁡(1−pte)], :=- E_t [y_t^prior p_t^e+(1-y_t^prior) (1-p_t^e) ], (113) ℒactor _actor :=−∑j∈j+λpriorℒprior. :=- _j J_j+ _priorL_prior. The coefficient λprior _prior is fixed within each run; its family-specific value is given in Appendix G.1. The binary cross-entropy term updates only the restart actor; its pre-action labels are held fixed during optimization. The rollout stores the observations, raw actions, masked log-densities, rewards, success indicators, critic values, and pre-action labels, τ=(ot,at,ℓt,rt,successt,Vϕ(ot),ytprior)t=1T.τ= \ (o_t,a_t, _t,r_t,success_t,V^φ(o_t),y_t^prior ) \_t=1^T. (114) F.4 Preconditioning We apply 10 iterations of Ruiz scaling (Ruiz, 2001), followed by Pock–Chambolle scaling with exponent 11 (Pock & Chambolle, 2011). Let Dr∈ℝm×mD_r ^m× m and Dc∈ℝn×nD_c ^n× n be the resulting positive diagonal row and column scalings, with Dr=(DP00DG)D_r= pmatrixD_P&0\\ 0&D_G pmatrix partitioned by equality and inequality constraints. The scaled problem data and their relation to the original variables are Asc A^sc =DrADc, =D_rAD_c, bsc b^sc =Drb, =D_rb, csc c^sc =Dcc, =D_cc, (115) Psc P^sc =DPPDc, =D_PPD_c, qsc q^sc =DPq, =D_Pq, Gsc G^sc =DGGDc, =D_GGD_c, hsc h^sc =DGh, =D_Gh, lsc l^sc =Dc−1l, =D_c^-1l, usc u^sc =Dc−1u. =D_c^-1u. x=Dcxsc,y=Drysc.x=D_cx^sc, y=D_ry^sc. (116) In this subsection, superscript scsc explicitly identifies scaled quantities, and unsuperscripted data and z=(x⊤,y⊤)⊤z= (x ,y ) refer to the original LP. The solver notation elsewhere omits scsc; calls to p,d,gp,d,g, and MrelM^rel include the coordinate conversion in Eq. (116). PDHG operates in scaled coordinates, while convergence metrics, rewards, and termination are evaluated for the original LP. The residual transformations and gap invariance are rpsc=Drrp,rdsc=Dcrd,r_p^sc=D_rr_p, r_d^sc=D_cr_d, (117) rgsc=rg.r_g^sc=r_g. (118) As in Eq. (6), infinite-bound terms are omitted from the computed gap. Original residual norms can be recovered from scaled residuals, ∥rp∥∞ r_p _∞ =∥Dr−1rpsc∥∞=maxi⁡|(rpsc)i|(Dr)ii, = D_r^-1r_p^sc _∞= _i |(r_p^sc)_i |(D_r)_i, (119) ∥rd∥∞ r_d _∞ =∥Dc−1rdsc∥∞=maxj⁡|(rdsc)j|(Dc)jj. = D_c^-1r_d^sc _∞= _j |(r_d^sc)_j |(D_c)_j. The normalized primal and dual metrics, including bound and dual-sign violations, are therefore p⁡(z)=log⁡(max⁡‖Dr−1rpsc‖∞,‖Dc[lsc−xsc]+‖∞,‖Dc[xsc−usc]+‖∞1+max⁡‖b‖∞,‖Dr−1Ascxsc‖∞).p(z)= ( \\|D_r^-1r_p^sc\|_∞,\|D_c[l^sc-x^sc]^+\|_∞,\|D_c[x^sc-u^sc]^+\|_∞ \1+ \\|b\|_∞,\|D_r^-1A^scx^sc\|_∞ \ ). (120) d⁡(z)=log⁡(max⁡‖Dc−1rdsc‖∞,‖DG[−ysc,in]+‖∞1+max⁡‖c‖∞,‖Dc−1(Asc)⊤ysc‖∞).d(z)= ( \\|D_c^-1r_d^sc\|_∞,\|D_G[-y^sc,in]^+\|_∞ \1+ \\|c\|_∞,\|D_c^-1(A^sc) y^sc\|_∞ \ ). (121) The primal objective is also invariant, (csc)⊤xsc=c⊤x(c^sc) x^sc=c x, so the gap metric is g⁡(z)=log⁡(|rgsc(zsc)|1+max⁡|(csc)⊤xsc|,|(csc)⊤xsc−rgsc(zsc)|).g(z)= ( |r_g^sc(z^sc)|1+ \|(c^sc) x^sc|,|(c^sc) x^sc-r_g^sc(z^sc)| \ ). (122) These are evaluation formulas for the original metrics, with the same 10−2010^-20 log floor as in Appendix F.2. They use only scaled residuals, cached products, and diagonal scalings, requiring no additional multiplication by A or A⊤A . The merit MrelM^rel, tolerance test, and reward in Eq. (8) use these original-coordinate metrics. Preconditioned and unpreconditioned runs are thus evaluated at the same original-LP tolerance. F.5 Matrix–vector product caching For the primal and dual blocks of the vectors used in Eq. (5a), we retain (Axt,Ax^t,Ax⋆,Atx)and(A⊤yt,A⊤y^t,A⊤y⋆,A⊤ty).(Ax_t,\;A x_t,\;Ax_ ,\;AK_t^x) (A y_t,\;A y_t,\;A y_ ,\;A K_t^y). Applying the acceleration coefficients to these products reconstructs Ax^t+1A x_t+1 and A⊤y^t+1A y_t+1 without another application of either matrix. The newest history products are Axt+1−AxtAx_t+1-Ax_t and A⊤yt+1−A⊤ytA y_t+1-A y_t. The remaining columns are shifted. At a restart, the products for the new anchor are copied from those of zt+1z_t+1, and the history matrices and their products are reset together. For t≥1t≥ 1, define the PDHG direction qt:=zt−z^tq_t:=z_t- z_t, with primal and dual blocks qtx:=xt−x^tq_t^x:=x_t- x_t and qty:=yt−y^tq_t^y:=y_t- y_t. The implementation applies the matrices to the blocks of qt+1q_t+1. The only new sparse products are Aqt+1xAq_t+1^x for the dual update and A⊤qt+1yA q_t+1^y for the residual features. The output products are recovered as Axt+1=Ax^t+1+Aqt+1xAx_t+1=A x_t+1+Aq_t+1^x and A⊤yt+1=A⊤y^t+1+A⊤qt+1yA y_t+1=A y_t+1+A q_t+1^y. The primal update uses the reconstructed A⊤y^t+1A y_t+1. For the observation, we also cache the unweighted Gram matrices (tx)⊤tx(K_t^x) K_t^x and (ty)⊤ty(K_t^y) K_t^y. Only their newest row and column require fresh inner products when a history direction is inserted. The remaining submatrices are shifted, and both Gram matrices are zeroed at restarts. GA-PDHG thus uses two new sparse matrix–vector products per iteration, with O⁡((m+n)k+k2)O ( (m+n )k+k^2 ) additional work for acceleration and cached geometry updates. GALLOP additionally evaluates the policy at each iteration. F.6 Full GALLOP algorithm Algorithm 1 collects one training episode of GA-PDHG under policy πθ _θ on an LP instance I, storing its T transitions in τ. The critic VϕV^φ supplies value estimates for training. Algorithm 2 trains the GALLOP policy, with ℬB denoting the buffer of valid transitions from all episodes in a PPO iteration. Feature definitions and policy losses are provided in Appendix F.2 and Appendix F.3, respectively. Initialization takes one unaccelerated PDHG step before the first controlled iteration. The history stores only output differences inserted after initialization or the latest restart, with unavailable columns set to zero. In the history update, t(:,:k−1)K_t(:,1\!:\!k-1) denotes its first k−1k-1 columns and is empty for k=1k=1. The returned VT+1ϕV_T+1^φ is the terminal value or truncation bootstrap used by GAE. The cap TmaxT_ counts controlled iterations and excludes the initial unaccelerated PDHG step. We assume ‖A‖2>0\|A\|_2>0 for the step-size formulas and require ϵtol≥10−20 _tol≥ 10^-20 to respect the numerical floor in the log residuals. Algorithm 1 GALLOP training rollout 1: Positive integer cap TmaxT_ , memory k≥1k≥ 1, tolerance ϵtol _tol, and controls and features specified in Appendix F.1 and Appendix F.2 2: function GALLOP-Rollout(πθ,Vϕ,I _θ,V^φ,I) 3: Initialize z⋆=z^1=z0=0z_ = z_1=z_0=0, weight ω0 _0, s←1s← 1, t←1t← 1, τ←∅τ← 4: z1←ℋ⁡(z^1,ω0,η0)z_1 ( z_1; _0, _0 ), 1←0K_1← 0 5: M⋆←M1fpM_ ← M_1^fp by Eq. (98); initialize observation baselines as in Appendix F.2 6: Form o1o_1 using Appendix F.2 7: while t≤Tmaxt≤ T_ do 8: Evaluate ytpriory_t^prior and Vϕ(ot)V^φ(o_t); sample ata_t using Appendix F.1 9: Record ℓt _t by Eq. (111); set ηt,ωt _t, _t by Eq. (92) and Eq. (93) 10: if et=0e_t=0 then 11: Form z^t+1 z_t+1 by Eq. (5a), using the transformed acceleration coefficients 12: else 13: z^t+1←zt z_t+1← z_t 14: end if 15: zt+1←ℋ⁡(z^t+1,ωt,ηt)z_t+1 ( z_t+1; _t, _t); measure Mt+1fpM_t+1^fp by Eq. (98) 16: if et=1e_t=1 then 17: z⋆←zt+1z_ ← z_t+1, M⋆←Mt+1fpM_ ← M_t+1^fp, t+1←0K_t+1← 0, s←0s← 0 18: else 19: t+1←[zt+1−ztt(:,:k−1)]K_t+1← [z_t+1-z_t _t(:,1\!:\!k-1) ], s←s+1s← s+1 20: end if 21: Evaluate the KKT metrics and rtr_t by Eq. (8) 22: successt←[Mrel(zt+1)≤logϵtol]success_t 1 [M^rel(z_t+1)≤ _tol ] 23: Append the transition to τ as in Eq. (114) 24: t←t+1t← t+1; form oto_t using Appendix F.2 25: if successt−1=1success_t-1=1 then 26: break 27: end if 28: end while 29: T←t−1T← t-1 30: if successT=1success_T=1 then 31: VT+1ϕ←0V_T+1^φ← 0 32: else 33: VT+1ϕ←Vϕ(oT+1)V_T+1^φ← V^φ (o_T+1 ) ⊳ bootstrap at the episode cap 34: end if 35: return τ,VT+1ϕτ,V_T+1^φ 36: end function Algorithm 2 PPO training of GALLOP 1: LP dataset =I(i)i=1ND=\I^(i)\_i=1^N, number of PPO iterations J, history length k, 2: control and standard-deviation bounds, γ,κ,εclip,λprior,cV,χ,ϵtol,Tmax,ϵgeomγ,κ, ^clip, _prior,c_V,χ, _tol,T_ , _geom 3: Initialize five actors, global log standard deviations, critic VϕV^φ, and their optimizers 4: for j=1,…,Jj=1,…,J do 5: (θold,ϕold)←(θ,ϕ)( _old, _old)←(θ,φ), ℬ←∅B← 6: for i=1,…,Ni=1,…,N in batches do 7: τ(i),VT(i)+1ϕold←GALLOP-Rollout(πθold,Vϕold,I(i))τ^(i),V_T^(i)+1 _old← [rgb]0,0,0GALLOP-Rollout( _ _old,V _old,I^(i)) 8: Compute Ψ^t(i) _t^(i) by Eq. (108), using VT(i)+1ϕoldV_T^(i)+1 _old for the final bootstrap 9: Append valid rollout records and their advantages to ℬB 10: end for 11: Compute Ψ~t _t and R^t R_t over the full buffer by Eq. (109) 12: for several epochs of minibatches from ℬB do 13: Update all five actors on ℒactorL_actor in Eq. (113) 14: Update the critic on cVℒVc_VL_V in Eq. (110) 15: end for 16: end for 17: return θ,ϕθ,φ Appendix G Details on Experiments G.1 Implementation details Hardware and numerical precision. All experiments use an NVIDIA GeForce RTX 4090 GPU with 24 GB of memory and an Intel Core i9-13900K CPU. The implementation uses JAX 0.10.2 with CUDA 12 and Flax 0.12.8; the external baseline uses MPAX 0.2.4. The policy networks compute in float32, while the PDHG state, residuals, rewards, and advantage estimates use float64. Optimization settings. Table 12 lists the settings of the runs that produced the selected checkpoints. All networks use Adam with a constant learning rate, four PPO epochs per rollout buffer, and a separate gradient-norm clip of 0.50.5 per network. All families use k=4k=4, γ=1γ=1, κ=0.99κ=0.99, PPO clipping εclip=0.2 ^clip=0.2, and no entropy bonus. Table 12: Settings of the training runs used in the six-family comparison. Family Learning rate PPO minibatch cVc_V λprior _prior LωL_ω Transport 5×10−45× 10^-4 16,384 0.5 0.02 4.6 Auction 5×10−45× 10^-4 4,096 0.5 0 2.0 Min-cost flow 10−410^-4 16,384 0.5 0.02 4.6 Robust inventory 10−410^-4 16,384 0.5 0.05 4.6 L∞L_∞ regression 5×10−45× 10^-4 16,384 0.2 0.01 4.6 Random LP 10−410^-4 16,384 0.5 0.02 4.6 Exploration standard deviations. Each continuous control group has a learned global log standard deviation, shared across observations and updated with its actor. For j∈α,β,ζ,η,ωj∈\α,β,ζ,η,ω\, let ℓj _j denote this learned scalar. The distribution uses σj=exp⁡(clip⁡(ℓj,log⁡σjmin,log⁡1.5)). _j= ( ( _j, _j , 1.5 ) ). (123) Table 13 gives the initialization and lower bounds from the selected runs. The two relaxation latents share σα _α, and all k history coordinates share σβ _β. Inference replaces the Gaussian draws by their means. Table 13: Gaussian exploration settings. Both tuples use the order (α,β,ζ,η,ω) (α,β,ζ,η,ω ); the common upper bound is 1.51.5. Family Initial standard deviations Minimum standard deviations Transport (0.3,0.2,0.4,0.4,0.2) (0.3,0.2,0.4,0.4,0.2 ) (0.2,0.1,0.2,0.2,0.1) (0.2,0.1,0.2,0.2,0.1 ) Auction (0.4,0.2,0.4,0.6,0.2) (0.4,0.2,0.4,0.6,0.2 ) (0.001,0.001,0.001,0.12,0.001) (0.001,0.001,0.001,0.12,0.001 ) Min-cost flow (0.2,0.3,0.3,0.4,0.2) (0.2,0.3,0.3,0.4,0.2 ) (0.05,0.2,0.05,0.05,0.05) (0.05,0.2,0.05,0.05,0.05 ) Robust inventory (0.05,0.05,0.07,0.1,0.07) (0.05,0.05,0.07,0.1,0.07 ) (0.01,0.01,0.01,0.02,0.02) (0.01,0.01,0.01,0.02,0.02 ) L∞L_∞ regression (0.4,0.2,0.4,0.5,0.2) (0.4,0.2,0.4,0.5,0.2 ) (0.2,0.1,0.05,0.3,0.1) (0.2,0.1,0.05,0.3,0.1 ) Random LP (0.1,0.1,0.1,0.1,0.1) (0.1,0.1,0.1,0.1,0.1 ) (0.02,0.01,0.02,0.02,0.01) (0.02,0.01,0.02,0.02,0.01 ) Batches across problem sizes. LPs with the same dimensions are grouped into a size bucket and solved in parallel within that bucket. Each PPO iteration collects one episode per training LP from every bucket before updating the policy. The valid transitions are then combined, their advantages are normalized over the full buffer, and shuffled PPO minibatches are formed using the sizes in Table 12. Thus rollout batches have fixed LP dimensions, while PPO minibatches can mix observations from different problem sizes because the observation and action dimensions depend only on k. We evaluate the same trained policy separately on each held-out and test bucket, including sizes larger than the training sizes listed in Appendix D.6.

Instant research discovery