Practical Large-Scale Linear Programming using Primal-Dual Hybrid Gradient
Authors: David Applegate
Mateo Diaz Diaz
Oliver Hinder
Haihao Lu
Miles Lubin
Brendan O'Donoghue
Warren Schudy
Venue: NeurIPS 2021 Poster (2021)
Keywords: linear programming
convex optimization
operations research
first-order methods
large scale optimization

Abstract: We present PDLP, a practical first-order method for linear programming (LP) that can solve to the high levels of accuracy that are expected in traditional LP applications. In addition, it can scale to very large problems because its core operation is matrix-vector multiplications. PDLP is derived by applying the primal-dual hybrid gradient (PDHG) method, popularized by Chambolle and Pock (2011), to a saddle-point formulation of LP. PDLP enhances PDHG for LP by combining several new techniques with older tricks from the literature; the enhancements include diagonal preconditioning, presolving, adaptive step sizes, and adaptive restarting. PDLP improves the state of the art for first-order methods applied to LP. We compare PDLP with SCS, an ADMM-based solver, on a set of 383 LP instances derived from MIPLIB 2017. With a target of $10^{-8}$ relative accuracy and 1 hour time limit, PDLP achieves a 6.3x reduction in the geometric mean of solve times and a 4.6x reduction in the number of instances unsolved (from 227 to 49). Furthermore, we highlight standard benchmark instances and a large-scale application (PageRank) where our open-source prototype of PDLP, written in Julia, outperforms a commercial LP solver.

Practical Large-Scale Linear Programming using Primal-Dual Hybrid Gradient
David Applegate
Mateo Díaz California Institute of Technology"
Oliver Hinder
Google Research
Google Research University of Pittsburgh ohinder&pitt. edu
dapplogate@google..o
mateodd@caltech. edu
Haihao Lu
Miles Lubin Google Research mlubin@google.. com
Brendan O'Donoghue
University of Chicago
DeepMind bodomoghue@ddeepminddccoo com
haihao. ww@chicagobooth.ccomm edu
Warren Schudy Google Research wachudy8google. com
Abstract
We present PDLP, practical first-order method for linear programming (LP) that can solve to the high levels of accuracy that are expected in traditional I.P applications. In addition, it can scale to very large problems because its core operation is matrix-vector multiplications. PDLP is derived by applying the primal- dual hybrid gradient (PDHG) method, popularized by Chambolle and Pock (2011), to a saddle-point formulation of LP. PDLP enhances PDHG for LP by combining several new techniques with older tricks from the literature; the enhancements include diagonal preconditioning, presolving, adaptive step sizes, and adaptive restarting. PDLP improves the state of the art for first-order methods applied LP. We compare PDLP with SCS, an ADMM based solver, on a set of 383 LP instances derived from MIPLIB 2017. With target of 10-8 relative accuracy and hour time limit, PDLP achieves a 6.3x reduction in the geometric mean solve times and a 4.6x reduction in the number of instances unsolved (from 227 to 49). Furthermore, we highlight standard benchmark instances and large-scale application (PageRank) where our open-source prototype PDLP, written Julia, outperforms a commercial LP solver.
Introduction
First-order methods (FOMs), which use gradient and not Hessian information, are now applied as standard practice in many areas of optimization [12]. A known weakness of FOMs is the tailing-off off effect, where FOMs quickly find moderately accurate solutions, but progress towards an optimal solution slows down over time. While moderately accurate solutions are often sufficient for large machine learning applications, other applications traditionally demand higher precision. One such area is Linear Programming (LP), the focus of this work.
LP a fundamental class of optimization problems in applied mathematics, operations research, and computer science with a huge range of applications, including mixed-integer programming. scheduling, network flow, chip design, budget allocation, and many others [17, 22, 65,69) 69). Software
Part the work was done Cornell and Google Research.
Part of the work was done at Google Research.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

for solving LP problems, called LP solvers, originated in the earliest days of computing, predating the invention of operating systems [55]. The state-of-ther methods for LP, namely Dantzig's simplex method [22,23] 23] and interior-point (or barrier) methods [5]], are quite mature and reliable at delivering highly accurate solutions. These widely successful methods have left little room for FOMs to make inroads. Furthermore, practitioners who use LP solvers are not accustomed to reasoning about the trade-off between accuracy and computing times typically intrinsic to FOMs.
In this paper, we provide evidence that, if properly enhanced, FOMs can obtain high quality solutions to LP problems quickly. Indeed, there's 's reason to expect this, as authors have developed FOMs for LP with linear rates of convergence [24, 32, 70,7000777 71]. On the other hand, the linear rates depend on potentially loose and hard-to-con card-to- constants; hence, tailing off may still be observed in practice. To our knowledge, ours is the first work to combine both theoretical enhancements with practical heuristies, demonstrating their combined effectiveness with extensive computational experiments on standard benchmark instances. In fact, our experiments will expose substantial gap between algorithms presented in the literature and what's needed to obtain good performance.
Starting from a baseline pril hybrid gradient (PDHG) method [19] applied to saddle point formulation of LP, we develop series of algorithmic improvements. These enhancements include adaptive restarting [7], dynamic primal-dual step size selection [36,37] 37], presolving techniques [1], and diagonal preconditioning (data equilibration)
[33]. Most of these enhancements, while inspired by existing literature, are novel. We name our collection of enhancements PDLP (PDHG for LP).
The impact of these improvements is substantial. For example, on 383 LP instances derived from the MIPLIB 2017 collection [34], our implementation of a baseline version of PDHG solved only 50 problems to 10-$ relative accuracy given a limit of approximately 100,000 iterations per problem. By contrast, PDLP solves 283 of the 383 problems under the same conditions. We demonstrate that PDLP outperforms FOM baselines and, in a small number of cases, obtains performance competitive with commercial LP solver.
Although not the focus of this paper, we believe that our results open the door to new set of possibilities and computational trade-offss when solving LP problems. PDLP has the potential to solve extremely large scale instances where the simplex method and interior-point methods are unable to run because of their reliance on matrix factorization. Since PDLP uses matrix-vector operations at its core, it can effectively run on multi-threaded CPUs, GPUs [68], or distributed clusters [26]. Furthermore, a GPU implementation of PDLP could efficiently solve batches of similar problems, setup that has already been successfully applied with other optimization algorithms in applications like strong branching [46] and training neural networks that contain optimization layers [5].
The remainder of this section focuses on related work. Section 2 introduces LP and PDHG. Section 3 describes the set of enhancements that define PDLP. Section presents numerical experiments, and Section 5 concludes and outlines future directions.
Outline.
1.1 Literature review
PDHG PDHG was first developed by Zhu and Chan [72], with subsequent analysis and extension by number of authors [3.18.19..2 19,21,27,38 38.. PDHG is closely related the Hurwizz method [8]. PDHG is a form of operator-splitinng [11, 64] and can be interpreted as a variant of the alternating directions method of multipliers (ADMM) and Douglas- Rachford splitting (DRS) [16, 56), which themselves are both instantiations of the proximal point method [25. 62). As opposed to ADMM or DRS, PDHG is 'matrix-fred in that the data matrix is only used for matrix-vector multiplications. This allows PDHG to scale to problems even larger than those tackled by these other techniques, and to make better use of parallel and distributed computation.
FOM-based solvers Recent interest in large-scale cone programming has sparked the development several first-order solvers based on competing methods. ProxSDP [66] is solver for semidefinite programming based on PDHG. Solvers based on Nesterov's accelerated gradients [49] include TFOCS [14], and FOM which is a suite of solvers employing both gradient and proximal algorithms [13]. Solvers based on operator splitting techniques like ADMM include SCS [52-54], OSQP [67], POGS [28], and COSMO [30]. Of these both SCS and POGS offer matrix-fres implementation where the linear system, that arises from the proximal operator used in ADMM, is solved using the conjugate gradient method. However, we shall show experimentally that our method can be significantly faster and more robust than this approach. Finally, [4] considers applying a truncated semismooth Newton method to the system of equations defining fixed point of the SCS operator.

---

FOMs for LP Lan, Lu and Monteiro [40] and Renegar [61] develop FOMs for LP as special case of semidefinite programming, with sublinear convergence rates. The FOM-based solvers above all apply to more general problem classes like cone programming or quadratic programming. In contrast, some of the enhancements that constitute PDLP are specialized, either theory or practice, for LP (namely restarts [7] and presolving). A number of authors [24, 42, 70, 71] have proposed linearly convergent FOMs for LP: to our knowledge, none have been subject of a comprehensive computational study. ECLIPSE [10] solves huge-scale industrial LP problems by accelerated gradient descent, without presenting comparisons on standard test problems. Lin et al. [42] propose an ADMM-based interior point method. In contrast with PDLP which solves to high accuracy (i.e., 10 relative error), [42] perform experiments with 10 and 10 10-5 relative error. SNIPAL [41] semismooth Newton method based the proximal augmented Lagrangian. SNIPAL has fast asymptotic convergence, yet, to get good performance, the authors use ADMM for warm-starts. Given PDLP's favorable comparisons with SCS, it's plausible that PDLP could provide more effective warm-start. Finally, Pock and Chambolle [59] apply PDHG with diagonal preconditioning to a limited set of test LP problems and Applegate al. [6] show how to extract infeasibility certificates when applying PDHG to LP.
2 Preliminaries
In this section, introduce the notation throughout the paper, summarize the LP formulations we solve, and introduce the baseline PDHG algorithm.
Notation. Let R denote the set of real numbers, R+ the set of nonnegative real numbers, and R the set of nonpositive real numbers. Let N denote the set of natural numbers (starting from one). Let Ip denote the Ep norm for a vector, and let || Il2 denote the spectral norm for a matrix. For a vector u E IR", we use wt and v for their positive and negative parts, i.e., 2 = max(0, e. and v; = m0, v. The symbol Vim denotes the vector with the first m components v. The symbols K2, and K. correspond to the ith column and jth row of the matrix K, respectively. The symbol denotes the vector of all ones. Given a convex set X, we use projx to denote the map that projects onto X.
Linear Programming.
We solve primal-dual LP problems of the form:
minimize #ER"
c
maximize yER "",,,
q y -
-
subject to:
Gr h
subject to:
c- K KTy= y=
(1)
Ax h
V1-m: 20 ΛEA,
where G E R X A E m2 c TR" h E Rm b E RRm L E (
KT = (G A A q = b and
{0} 1 = -0, Ui = 00,
R- 1, = 00, E IR
A = A .... An
:=
R R. DC
R otherwise
the set of variables À such that the dual objective is finite. This pair of primal-dual problems is equivalent to the saddle-point problem:
- min max L(x, v) = c..T y -o--yKke Ka Tx +ou+q" T
(2)
zEX wEY
with := fx B = u}. and == { E R : Vism, 2 0}.
PDHG. When specialized to (2), the PDHG algorithm takes the form:
x k+1 = proj(a*
T(C - K g)))))
x
(3)
++1 A+1 = y
proj(y* + + K K (2x k+1 1- -
k+1
3

---

where T, 7 > 0 are primal and dual step sizes, respectively. PDHG is known to converge to optimal solution when rar K = [20, 21]. We reparameterize the step sizes by
- m// and
o = wi
with E (0, 0o) and
w E (0)
(4)
We call w E (0, o)) the primal weight, and E (0, oo) the step size. Under this reparameterization PDHG converges for all 1) [|2 This allows us control the scaling between the primal and dual iterates with a single parameter w. We use the term primal weight to describe w because weights the primal variables in the following
|
V w wl|
This norm plays a role in the theory for PDHG [20] and later algorithmic discussions.
For the baseline PDHG algorithm that we use for comparisons, we consider two simple choices for 7 and w. For the step size, set 0.9/ K||2 where ||/]|2 is estimated via power iteration, and for the primal weight we set w 1; this similar the default parameters in the standard PDHG implementation in ODI. [2].
Practical algorithmic improvements
In this section, we detail these enhancements, and defer further experimental testing of them to Section 4 and ablation studies Appendix C. While enhancements inspired by theory, focus is on practical performance. The algorithm as whole has no convergence guarantee, although some individual enhancements do; see Section 3.6 for further discussion.
Algorithm 1 presents pseudo-code for PDLP after preprocessing steps. We modify the step sizes (Sec- tion 3,1), add restarts (Section 3.2), and dynamically update the primal weights (Section 3.3). Before running Algorithm I we apply presolve (Section 3.4) and diagonal preconditioning (Section 3.5). There are some minor differences between the pseudo-code and the actual code. In particular, we only evaluate the restart or termination criteria (Line 10) every 40 iterations. This reduces the associated overheads with minimal impact on the total number of iterations. also check the termination criteria before beginning the algorithm if we detect numerical error.
PDLP (after preconditioning and presolve) Input: An initial solution 20 2 Initialize outer loop counter n + 0, total iterations + O, step size 170.0 < 1/|| primal weight w +-
Algorithm 1:
q);
repeat
+ 0:
repeat
n n,t+1 ,i jnttt1
IdaptiveSSeeOOFPPH( "mt ,k)
gnd+1 E.1
26 n1++1 t- GeGeeeestartCCCndddteeeeeeeeee
gn,et 2n.0);
t 1-t++ 1;;
10
until restart or termination criteria holds;
restart the outer loop. =+10 + n.t t- 1;
11
12
w Primall We wegeehhhttttteeeee n00 =n-1,0 10=_1,,,"-1- "-1)
13
until termination criteria holds;
14
Output: 00
3.1
Step size choice
The convergence analysis [20, Equation (15) PDHG (equation (3)) relies on small constant step
size
|| k+ 2k12
1 5 2 - y*) T K (x++1 x*
(5)

---

Algorithm 2: One step of PDHG using our step size heuristic
Function AdaptiveStepOfPDHG (2 w", sim.t k):
(ær y) 2 <- f/m
3
for = I,..., 00 do - projx
2 (c K y))
wn
y projy (y + +
(2x' =))));
[[ -2 2.0 2)1E
7 2w [(w'-x)
min ((1 (k: + 1)-0033) y7, (1 +k+ 1))
if 1 S rj then
return (a', (2, n. n.11 end
11
17 + r/:
12
end
where z² (at ))) Classically one would ensure (5) by picking = 2 This is overly pes-
simistic and requires estimation of ||2 Instead AdaptiveStep@fPDHG adjusts v dynamically to ensure that (5) is satisfied. If (5) isn't satisfied, we abort the step: i.c., we reduce 7), and try again. If (5) is satisfied we accept the step. This is described in Algorithm Note that in Algorithm
1-o(1) holds as ? 2 || holds always, and from this one can show the resulting step size rj 2 | k + D0.
Our step size routine compares favorably in practice with the line search by Malitsky and Pock [43] (See Appendix C.I).
3.2 Adaptive restarts
In PDLP, we adaptively restart the PDHG algorithm in each outer iteration. The key to our restarts at the n-th outer iteration is the normalized duality gap at which for any radius (0, oo) is defined by
1 maximize
(2) = r (3.)) [L(x, 0) L(F, v)),
introduced by [7]. Unlike the standard duality gap
maximize[C(z, (2:
v),.
the normalized duality gap is always a finite quantity. Furthermore, for any value of r and w", the normalized duality gap (=) ) is if and only if the solution is an optimal solution to (2) [7]: thus, provides a valid metric for measuring progress towards the optimal solution. The normalized duality gap is computable in linear time [7]. For brevity, define pn pn(22, zve) the normalized duality gap at with radius || Zeellum, i.e.,
Me (,, 2vet) Pfiz-zal wn (2),
where žee is user-chosen reference point.
Choosing the restart candidate.
To choose the restart candidate z: n,,+1
we call
(2,.211141, 2000 < (2n.4+1 20,)
GetRestartCandidatee (.11+1 gn.t+ 2,00 200
gn,441 otherwise.
This choice is justified in Remark of [7].
Restart criteria. We define three parameters: Amuliccent E (0, 1), Inecessary E (0, Ssullicie) and Baiitciall E (0, 1). In PDLP we use Baulicient 0.9, Buecessary 0.1, and Bartificial 0.5. The algorithm
restarts if one of three conditions holds:
(Sufficient decay in normalized duality gap) ph (= n, n.)) (autfiiien/ - ==-1))
(ii) (Necessary decay + no local progress in normalized duality gap)
(nnttt 22 < Puoocssary pn (2 n2 2-11,0) and Mn (2 ,000 > Mn( ,2 200))
5

---

(iii) (Long inner loop) 2 uunniiaa
The motivation for (i) is presented in [7]: guarantees the linear convergence of restarted PDHG on LP problems. The second condition (ii) is inspired by adaptive restart schemes for accelerated gradient descent where restarts triggered if the function value increases [57]. The first inequality in (ii) provides safeguard for the second one, preventing the algorithm restarting every inner iteration or never restarting. The motivation for (iii) relates to the primal weights (Section 3.3). In particular, primal weight updates only occur after restart, and condition (iii) ensures that the primal weight will be updated infinitely often. This prevents a bad choice of primal weight in earlier iterations causing progress to stall for a long time.
3.3 Primal weight updates
The primal weight is initialized using
| Illla
| 11q122 > Ezero
Initiializeerrimmlleiiiig q) =
otherwise
where Ezcro is small nonzero tolerance. This primal weight update scheme guarantees scale in- variance. particular, in Appendix A we consider PDHG with Ezo = Initializen rrriaaizeerimaleeeee
17 0.9/1 |||2 and (c, q). In this simplified setting, we prove that if we multiply the objective, constraints, or the right hand side and variable bounds by a scalar then the iterate behaviour remain identical (up to scaling factor).
Algorithm 3: Primal weight update
Function PrimalWeightUpdarn (zn.0 nn -- A || x" x=-1.016 12,
2
A" = 1l y" yn11001 1.0
if A > Ezero and A > Ezem A"
then
return exp log A
+ (1 - 0) log(w""11)
else
return w-
end
Algorithm aims to choose the primal weight w such that distance optimality in the primal and dual is the same, i.e., "1( x², 0)(lum Re [I(0, y0, y²) [wm. By definition of I
[1 - x", x0",0) 1) = w m - 1212,, [1(0, y - -y') wn
1 y
Setting these terms equal yields nn 1 | -
1y""-="'" is unknown "y""-="'le 122 Of course, the quantity || -x*xxx"
beforehand, but we attempt to estimate it using A /// However, the quantity Ag A can change wildly from restart another, causing wn oscillate. To dampen variations in "h we first move to a log-scale where the primal weight is symmetric, i.e., log(1/w") (/m) = log(w" ()) and perform a exponential smoothing with parameter 0 [0, 1]. In PDLP, we use 0.5.
There are several important differences between our primal weight heuristic and literature [36, 37]. For example, [36,37] make relatively small changes the primal weights each iteration, attempting to balance the primal and dual residual. These changes have to be diminishingly small because, in our experience, PDHG may be unstable if they are too big. In contrast, in method the primal weight only updated during restarts, which in practice allows for much larger changes without instability issues. Moreover, scheme tries balance the weighted distance traveled in the primal and dual rather than the residuals [36.37].
3.4 Presolve
Presolving refers to transformation steps that simplify the input problem before starting the opti- mization solver. These steps span from relatively easy transformations such as detecting inconsistent bounds, removing empty rows and columns of K, and removing variables whose lower and upper bounds are equal, to more complex operations such as detecting duplicate rows in K and tightening bounds. Presolve is standard component of traditional I.P solvers [44]. We are not aware of presolve being combined with PDHG for LP. However, [41, 42] combine presolve with other FOMs.

---

As an experiment to measure the impact of presolve, we used PaPILO [29]. an open-source presolving library. For technical reasons, was easier to use PaPILO standalone executable than library. We simulate its effect by simply solving the preprocessed instances. Convergence criteria are evaluated with respect to the presolved instance, not the original problem.
3.5 Diagonal Preconditioning
Preconditioning is a popular heuristic in optimization for improving the convergence of FOMs. To avoid factorizations, we only consider diagonal preconditioners. Our goal is rescale the constraint matrix K (G, A) to K (G, A) Dj K D2 with positive diagonal matrices D₁ and D2, that the resulting matrix K is "well balanced". Such preconditioning creates new LP instance that replaces A, G, c, h, u, and in (1) with G, A, & Dz" xr ė D2c, (6, h) Di Di (b, h), û D2"
and D Common choices for D1 and D2 include:
No scaling: Solve the original LP instance (1) without additional scaling, namely D1 = D2 =
Sock-Chabbolle [59]: Pock and Chambolle proposed a family of diagonal preconditioners? for PDHG parameterized by a, where the diagonal matrices are defined by (D1)ss = VIl in PDLP (we also
for j 1,... M, + my and (D2) VIIK,, ala for i 1, ..... n. We tested a = 0 and a 2). This the baseline diagonal preconditioner in the PDHG literature.
Ruiz [63]: Ruiz scaling is a popular algorithm in numerical linear algebra to equilibrate matrices. In an iteration of Ruiz scaling, the diagonal matrices are defined as (Di) )as = V for 1,.... m7 m and (D2)a Vi for i L, ..... Ruiz [63] shows that if this rescaling is applied iteratively, the infinity norm of each row and each column converge to 1.
For the default PDLP settings, we apply a combination of Ruiz rescaling [63] and the preconditioning technique proposed by Pock and Chambolle [59]. In particular, we apply 10 iterations of Ruiz scaling and then apply the Bock-Chambolle scaling. To illustrate the effectiveness of our proposed scaling technique, compare it against these three common techniques in Appendix C.5.
3.6 Theoretical guarantees for the above enhancements
While PDLP's enhancements are motivated by theory, some of them may not preserve theoretical guarantees discussed below:
We do not have proof of convergence for the adaptive step size rule (Section 3.1).
One can show our restart criteria (Section 3.2) preserve convergence guarantees by modifying the proof of [7] to more general setting.
Primal weight updates (Section 3.3) do not readily preserve convergence guarantees, but conjecture that proof of convergence is possible they are updated infrequently.
Presolve (Section 3.4) and diagonal preconditioning (Section 3.5) preserve theoretical guarantees because they can be viewed as upplying PDHG to an LP instance with different data.
4 Numerical experiments
Our numerical experiments study the effectiveness of PDLP primarily with respect to traditional LP applications and benchmark sets. Section 4.1 describes the setup for the experiments. Section 4.2 demonstrates PDLP's s improvements over baseline PDHG. Section 4.3 compares PDLP with other FOMs. Section 4.4 highlights benchmark instances where PDLP outperforms commercial LP solver. Finally, Section 4.5 illustrates the ability PDLP to scale to a large application where barrier and simplex-based based solvers out of memory. The supplemental materials contain extensive ablation studies and additional instructions for reproducing the experiments.
DDiagonal preconditioning is equivalent to changing to a weighted l, norm in the proximal step of PDHG (weight defined by D2 and D; for the primal and dual respectively). Pock and Chambolle use this weighted perspective.
7

---

4.1 Experimental setup
Optimality termination criteria. PDLP terminates with an approximately optimal solution when the primal-dual iterates & E X, y E Y. À E A, satisfy:
la y - u _ - c 2| $ 1 - u T +
(6a)
Az - b || (11 Gz)+
< + l|all2)
(6b)
||c- y A112 S c1 + l|d|22
(6c)
where e E (0, oo) is the termination tolerance. Note that if (6) is satisfied with e 0, then by LP duality we have found an optimal solution [35]. Indeed, (6a) is the duality gap, (6b) is primal feasibility, and (6c) is dual feasibility. We use these criteria to be consistent with those of SCS [54]. The PDHG algorithm does not explicitly include a reduced costs variable X. Therefore, to evaluate the optimality termination criteria we compute À = Proja (c- KTy), All instances considered have an optimal primal-dual solution. We use c = 10 as a benchmark for high-quality solutions and 10 for moderately accurate solutions.
Benchmark datasets. We use three datasets to compare algorithmic performance. One is the LP benchmark dataset of 56 problems, formed by merging the instances from "Benchmark of Simplex LP Solvers", "Benchmark of Barrier LP solvers", and "Large Network-L. Benchmark" from [45]. We also created larger benchmark of 383 instances curated from LP relaxations of mixed-integer programming problems from the MIPLIB2017 collection [34] (see Appendix B) that we label MIP Relaxations. MIP belaxations was used extensively during algorithmic development, e.g., for hyperparameter choices; we held out LP benchmark as a test set. Finally, we also performed some experiments on the Netlib LP benchmark [31], an historically important benchmark that is no longer state the art for large-scale LP.
Software. PDLP is implemented in an open-source Julia [15] module available at https: //github. com/ researcrreeeechhheeeee Fiirstrrderrp..ii The module also contains baseline im- plementation of the extragradient method with many of the same enhancements as PDLP (labeled 'Enh. Extragradient') We compare with two external packages: SCS [54] version 2.1.3, an open- source generic cone solver based on ADMM, and Gurobi version 9.0.1, a state-of-the-art commercial IP solver. SCS supports two modes for solving the linear system that arises each iteration, direct method based on a cached LDL factorization (which is the default 'SCS') and an indirect method based on the conjugate gradient method (which label 'SCS (matrix-free))). All solvers are run single-threaded. SCS and Gurobi are provided the same presolved instances as PDLP.
Computing environment. We used two computing enviromments for our experiments: I) e2-nighmem-2 virtual machines (VMs) on Google Cloud Platform (GCP). Each VM provides two virtual CPUs and 16GB RAM. 2) A dedicated workstation with an Intel Xeon E5-2669 v3 proces- sor and 128 GB RAM. This workstation has a license for Gurobi that permits most one concurrent solve. Total compute time on GCP for all preliminary and final experiments was approximately 72, 000 virtual CPU hours.
Initialization. All first-order methods use all-zero vectors as the initial starting points.
Metrics. We use the term KKT passes to refer to the number of matrix multiplications by both K and KT. Given that the most expensive operation in algorithm is matrix- multiplication, this metric is less noisy than runtime for comparing performance between matrix-free solvers. SGM10 stands for shifted geometric mean with shift 10, which is computed by adding 10 to all data points, taking the geometric mean, and then subtracting 10. Unsolved instances are assigned values corresponding to the limits specified in the next paragraph.
Time and KKT pass limits. For Section 4.2 impose limit on the KKT passes of 100, 000. For Section 4.3 we impose a time limit of hour.
4.2 Impact of PDLP's improvements
The y-axes of Figure I display the SGM10 of the KKT passes normalized by the value for baseline PDHG. We can see, with the exception of presolve for LP benchmark at tolerance 10-4, each of our modifications described in Section 3 improves the performance of PDHG.
8

---

Irperveeeds
taksanss 26:04
10"
breeanc 15-4
Sermalad MT pasna 50/18
10
lourares 14.24
takzarce
toeeance 12:00
lourance
7a
10
10
Tx22-1
+iNdans eecsing eotrrll
KHH -redaat +scalig
+Stera
Fle
wiiaans eccling
-perral
(a) MIP Relaxations
(b) LP benchmark
(c) Netlib
Figure I: Summary of relative impact of PDLP's improvements
1.
50 Relasationa Saacires Tel 15:04
2.0
N.P Relanetions Sase hr Tel 15:04
1.0 M.P Relaxatior Bassites Tel 15:00
3.0
Hit Relanation Basoiies Tel 15-08
0.5
0.8
0.2
127
12
12
19
13"
12"
72"
19
12.
127
1.0
LP Derchmark Bese tree Tel 10:01
30-
LP Senchmerk Deedinee Tel 15-60
Bordhark Bosciigs fol 18 08
1.0
LP Derchmert Deselinee Tel
Tragtto peotere
0.5
107
32
12
10
22'
107
1.4
Vetto Seadines Tol 18-08
Net Sautthes To 1E:04
Net Sesetines Tol 18:00
Tertic Denel nea Tel 30 00
0.5
04
600
15'
101
23t
10
S5M15
PDLP
Enh. Extragradient
PDHG
SCS (matrix-free)
SCS
Figure 2: Number of problems solved for MIP Relaxations (top), LP benchmark (middle), and Netlib (bottom) datasets.
4.3 Comparison with other first-order baselines
We compared PDLP with several other first-order baselines: SCS [54], in both direct (default) mode and matrix-free mode, and our enhanced implementation of the extragradient method 339.48]. For SCS in matrix-free free mode, direct mode there is no reasonable measure of KKT passes for the factorization and direct solve, so we only measure running time. The comparisons are summarized in Figure 2.
include the KKT passes from the conjugate gradient solves; for SCS in
4.4 PDLP versus simplex and barrier
In this section, we test the performance of PDLP against the three methods available in Gurobi: barrier, primal simplex, and dual simplex. By default when provided multiple threads, Gurobi runs these three methods concurrently and terminates when the first method completes. We used default termination for Gurobi and set 10-8 for PDLP. We ran experiments with instances from the MIP Relaxe delaxations and LP benchmark. Although, for most instances, Gurobi outperforms PDLP, we found problems for which PDLP exhibits moderate to significant gains. Table gives examples of instances where our prototype implementation is within a factor of two of the best of the three Gurobi methods. While further improvements are needed for PDLP to truly compete with the portfolio of methods that Gurobi offers, we interpret these results as evidence that PDL. itself could be of value in this portfolio.

---

Table 1: Instances from MIP Relaxations (top) and benchmark (bottom) where PDLP is within a factor of 2 of the best of all Gurobi methods. Time to solve in seconds.
Instance
PDLP
Gurobi Barrier
Gurobi Primal Simp.
Gurobi Dual Simp.
ex9
1.6
102.6
181.3
47.6
genus-sym--g222
2.1
10.7
6.7
33.2
higiicoooll-aigio
72.6
243.8
>3600
>3600
neos-578379
1.4
0.7
1.7
1.8
rwth-timetable
1870.3
>3600
>3600
>3600
ex
4.9
63.1
16.8
7.9
nug08-3rd
2.2
3.2
2219.2
24.1
savsched1
35.9
25.9
56.0
261.3
Table Solve time for PageRank instances. Gurobi barrier has crossover disabled, 1 thread. PDLP and SCS solve to 10 relative accuracy. SCS is matrix-free free. Baseline PDHG is unable to solve any instances. Presolve not applied. OOM Out Memory. The number of nonzero coefficients per instance is x nodes) 18.
nodes
PDLP
SCS
Gurobi Barrier
Gurobi Primal Simp.
Gurobi Dual Simp.
10a
7.4 sec.
1.3 sec.
36 sec.
37 sec.
14 sec.
10a
sec.
38 sec.
7.8 hr.
9.3
>2 hr.
106
11 min.
25 min.
OOM
>24 hr.
107
5.4 hr.
3.8 hr.
4.5 Large-scale application: PageRank
Nesterov [50, equation (7.3)) gives LP formulation of the standard "PageRank" problem. Although the LP formulation is not the best approach to computing PageRank, it is a source of very large instances. For a random scalable collection of PageRank instances, we used Barabási-Albert [9] preferential attachment graphs with approximately three edges per node; see Appendix D for details. The results summarized in Table 2.
Conclusions and future work
We find our experimental results encouraging for the application of FOMs like PDHG LP. At minimum, they provide evidence against the claim that FOMs are useful only when moderately accurate solutions are desired. The practical success of our heuristics that lack theoretical guarantees provides fresh motivation for theoreticians to study these methods. It is important, well, to understand what drives the difficulty of some instances and how they could transformed solve more quickly. We hope the community will use the benchmarks and baselines released with this work as a starting point for further investigating new FOMs for LP. With additional algorithmic and implementation refinements, we believe that PDLP or similar approaches could become part of the standard toolkit for linear programming.
Acknowledgments and Disclosure of Funding
We thank Yura Malitsky for advice on parameter choices for the linesearch rule of [43].
The authors have no third-party funding or competing interests to declare.
10

---

References
[1] T. Achterberg, R. E. Bixby, Z. Gu, E. Rothberg, and D. Weninger. Presolve reductions in mixed integer programming. INFORMS Journal on Computing, 32(2):473-506, 2020.
[2] Adler, Kohr, and O. Öktem. Operator discretization library (ODL), Jan. 2017.
[3] A. Alacaoglu, O. Fercoq, and Cevher. On the convergence of stochastic primal-dual hybrid gradient. arXiv preprint
2019.
[4] A. Ali, E. Wong, and J. Z. Kolter. A semismooth Newton method for fast, generic convex programming. In International Conference on Machine Learning, pages 70-79. PMLR, 2017.
[5] B. Amos and J. Z. Kolter. OptNet: Differentiable optimization layer in neural networks. In
D. Precup and W. Tch, editors, Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 136-145. PMLR, 06-11 Aug 2017.
[6] D. Applegate, M. Díaz, H. Lu, and M. Lubin. Infeasibility detection with primal-dual hybrid gradient for large-scale scale linear programming. arXiv preprint arXiv:2102.04552) 2021.
[7] D. Applegate, O. Hinder, H. Lu, and M. Lubin. Faster First-Order Primal-Dual Methods for Linear Programming using Restarts and Sharpness. arXiv preprint wrXiv:210S.12715, 2021.
[8] K. J. Arrow, Hurwicz, and Uzawa. Studies linear and non-linear University Press, 1958.
ron-iineurpnggrammiig programming. Stanford
[9] A.-L. Barabási and R. Albert. Emergence of scaling in random networks. Science, 286(5439):509 512, 1999.
[10] K. Basu, A. Ghoting, R. Mazumder, and Y. Pan. ECLIPSE: An extreme-scale linear program solver for web-applications. In H. D. III and A. Singh, editors, Proceedings of the 37th Inter- national Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 704-714, Virtual, 13-18 Jul 2020. PMLR.
[11] H. H. Bauschke and P. I.. Combettes. Convex analysis and monotone operator theory in Hilbert spaces, volume 408. Springer, 2 edition, 2017.
[12] A. Beck. First-Order Methods in Optimization. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2017.
[13] A. Beck and N. Guttmann BBeck. FOM-a MATLAB toolbox of first-order methods for solving convex optimization problems. Optimization Methods and Software, 34(1):172-193, 2019.
[14] S. R. Becker, E. J. Candès, and M. C. Grant. Templates for convex cone problems with applications to sparse signal recovery. Matherattcal programming computation, 363) 2011.
[15] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah. Julia: A fresh approach to numerical computing. SIAM Review, 59(1655-99 2017.
[16] Boyd, N. Parikh, E. Chu, B. Peleato, and Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine learning. 3(1):1-122, 2011.
[17] S. Boyd and I. Vandenberghe. Convex optimization. Cambridge university press, 2004.
[18] A. Chambolle, M. J. Ehrhardt, P. Richtárik, and C.-B. Schonlieb. Stochastic primal-dual hybrid gradient algorithm with arbitrary sampling and imaging applications. SIAM Journal on Optimization, 2874):2783-2808 2018.
[19] A. Chambolle and T. Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of mathematical imaging and vision, 40(1): 201
[20] A. Chambolle and T. Pock. On the ergodic convergence rates of first-order primal-dual algorithm. Mathematical Programming, 159(1):253-28) 2016.

---

[21] L. Condat. A primal- dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms. Journal of Optimization Theory and Applications. 158(2):460 2013.
[22] G. Dantzig. Linear programming and extensions. Princeton university press, 2016.
[23] G. B. Dantzig. Origins of the simplex method. In A history of scientific computing, pages
141-151. Association for Computing Machinery, New York, NY, USA, 1990.
[24] J. Eckstein and D. P. Bertsekas. An alternating direction method for linear programming. Tech- nical Report LIDS-P-1967 Laboratory for Information and Decision Systems, Massachusetts Institute of Technology, 1990.
[25] Eckstein and D. P. Bertsekas. On the Douglas-Rachford point algorithm for maximal monotone operators. Mathematical Programming, 55(13))299- 318, 1992.
Rachford splitting method and the proximal
[26] J. Eckstein and G. Matyasfalvi. Efficient distributed- memory parallel matrix-vector multi- plication with wide or tall unstructured sparse matrices. arXiv preprint arXiv: 1812.00904, 2018.
[27] E. Esser, X. Zhang, and T. F. Chan. A general framework for a class of first order primal-dual algorithms for convex optimization in imaging science. SIAM Journal on Imaging Sciences, 3(4):1015-1046, 2010.
[28] C. Fougner and S. Boyd. Parameter selection and preconditioning for a graph form solver. In Emerging Applications of Control and Systems Theory, pages 41-61. Springer, 2018.
[29] G. Gamrath, D. Anderson, K. Bestuzheva, W.-K. Chen, L. Eifler, M. Gasse, P. Gemander, A. Gleixner, L. Gottwald, K. Halbig, et al. The SCIP optimization suite 7.0. ZIB-Report 20-10,
Zuse Institut Berlin, 2020.
[30] M. Garstka, M. Cannon, and P. Goulart, COSMO: A conic operator splitting method for large convex problems. In European Control Conference, 2019.
[31] D. M. Gay. Electronic mail distribution of linear programming test problems. Mathernatical Programming Society COAL Newsletter, 13:10-12, 1985.
[32] A. Gilpin, Pena, and T. Sandholm. First-order algorithm with (In(1/c) convergence for e-equilibrium in two-person zero-sum games. Mathematical programming, 133(1):279-298, 2012.
[33] P. Giselsson and S. Boyd. Linear convergence and metric selection for Douhlas-Rachford splitting and ADMM. IEEE Transactions on Automatic Control, 62(2):532-544, 2016
[34] A. Gleixner, G. Hendel, G. Gamrath, T. Achterberg, M. Bastubbe, T. Berthold, P. M. Christophel, K. Jarck, T. Koch, J. Linderoth, M. Lübbecke, H. D. Mittelmann, D. Ozyurt, T. K. Ralphs, D. Salvagnin, and Y. Shinano. MIPLIB 2017: Data-Driven Compilation of the 6th Mixed-Integer Programming Library. Mathematical Programming Computation, 2021.
[35] A. J. Goldman and A. W. Tucker. Theory of linear programming. Linear inequalities and related systems, 38:53-97, 1956.
[36] T. Goldstein, M. Li, and X. Yuan. Adaptive primal-dual splitting methods for statistical learning and image processing. In Advances in Neural Information Processing Systems, pages 2089-2097, 2015.
[37] Goldstein, M. Li, X. Yuan, E. Esser, and R. Baraniuk, Adaptive primal-dual hybrid gradient methods for saddle-poii point problems. arXiv preprint arXiv:1305.0546, 2013.
[38] B. He and X. Yuan. Convergence analysis of primal-dual algorithms for a saddle-poin problem: from contraction perspective. SIAM Journal on Imaging Sciences, 5(1):119-144) 2012.
[39] G. M. Korpelevich. The extragradient method for finding saddle points and other problems. Matecon, 12:747-756, 1976.
[40] G. Lan, Z. Lu, and R. D. C. Monteiro. Primal-dual first-order methods with (11/e) iteration- complexity for cone programming. Mathematical Programming, 126(1):1-29, Jan 2011.
12

---

[41] X. Li, D. Sun, and K.-C. Toh. An asymptotically superlinearly convergent semismooth new- ton augmented lagrangian method for linear programming. SIAM Journal on Optimization. 2020.
30(3):2410-244
[42] T. Lin. S. Ma, Y. Ye, and S. Zhang. An ADMM-based interior-point method for large-scale linear programming. Optimization Methods and Software, 36(2-3):3994244, 2021.
[43] Y. Malitsky and T. Pock. A first-order primal-dual algorithm with linesearch. SIAM Journal on Optimization, 28(1):41 1-43,, 2018.
[44] I. Maros. Computational Techniques of the Simplex Method. International Series in Operations Research & Management Science. Springer US, 2002.
[45] H. Mittelmann. Decision tree for optimization software. http://plato. asu edua..dd//aa diduguiide.. html, 2021.
[46] V. Nair, S. Bartunov, F. Gimeno, I. von Glehn, P. Lichocki, I. Lobov, B. O' O'Donoghue, N. Son- nerat, C. Tjandraatmadja. P. Wang. R. Addanki, T. Hapuarachchi, T. Keck, Keeling. P. Kohli, I. Kiena, Y. Li, O. Vinyals, and Y. Zwols. Solving Mixed Integer Programs Using Neural Networks. arXiv preprint urXiv:2012. 13349, Dec. 2020.
[47] I. Necoara, Y. Nesterov, and F. Glineur. Linear convergence of first order methods for non- strongly convex optimization. Mathematical Programming, 175(1):69-107, May 2019.
[48] A. Nemirovski. Prox-method with rate of convergence O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems. SIAM Journal on Optimization, 15(1):229-251 2004.
[49] Y. Nesterov. A method of solving convex programming problem with convergence rate ((1/k2). Soviel Mathematics Doklady, 27(2):372-376, 1983.
[50] Y. Nesterov. Subgradient methods for huge-scale optimization problems. Mathematical Programming, 146:275-297, 2014.
[51] Y. Nesterov and A. Nemirovskii. Interior-point polynomial algorithms in convex programming, volume 13. SIAM, 1994.
[52] B. O'Donoghue. Operator splitting for a homogeneous embedding of the monotone linear complementarity problem. arXiv preprint arXiv:2004.02177, 2020.
[53] B. O'Donoghue, F. Chu, N. Parikh, and S. Boyd. Conic optimization via operator splitting and homogencous self-dual embedding. Journal of Optimization Theory and Applications, 16463):1022-1088, June 2016.
[54] B. O'Donoghue, E. Chu, N. Parikh, and S. Boyd. SCS: Splitting conic solver, version 2.1.0.
https://githut com/ http//vxgrppsseee scs, Nov. 2017.
[55] W. Orchard- Hays. History of mathematical programming systems. IEEE Annals of the History of Computing, 6(3):296-311 2122 1984.
[56) D. O'Connor and L.. Vandenberghe. On the equivalence of the primal-dual hybrid gradient method and Douglas Douhlas-Rachford splitting. Mathematical Programming. (79(1):85- 2020.
[57] B. O' D'nnoghue and E. Candes. Adaptive restart for accelerated gradient schemes. Foundations of computational mathematies, 15(3):715-732 2015.
[58] N. Parikh and S. Boyd. Proximal algorithms. Foundations and Trends in Optimization, 133):127-23, 2014.
[59] T. Pock and A. Chambolle. Diagonal preconditioning for first order primal-dual algorithms in convex optimization. In / International Conference Computer Vision, pages 1762-1769. IEEE, 2011.
[60] T. Pock, D. Cremers, H. Bischof, and A. Chambolle. An algorithm for minimizing the mumford- shah functional. In 2009 IEEE 12th International Conference on Computer Vision, pages 1133 1140. IEEE, 2009.
[61] Renegar. Accelerated first-order methods for hyperbolic programming. Mathematical Programming, 173(1):1-35, Jan 2019.
13

---

[62] R. T. Rockafellar. Monotone operators and the proximal point algorithm. SIAM journal on control and optimization, 14(5):877-898, 1976.
[63] D. Ruiz. A scaling algorithm to equilibrate both rows and columns norms in matrices. Technical report, CM-P00040415, 2001.
[64] E. K. Ryu and S. Boyd. Primer on monotone operator methods. Appl. Comput. Math 15(1):3-43, 2016.
[65] A. Schrijver. Theory of linear and integer programming. John Wiley & Sons, 1998.
[66] M. Souto, D. Garcia, and Á. Veiga. Exploiting low-rank structure in semidefinite programming by approximate operator splitting. Optimization, pages 1-28, 2020.
[67] B. Stellato, G. Banjac, Goulart, A. Bemporad, and S. Boyd. OSQP: an operator splitting solver for quadratic programs. Mathematical Programming Computation, 12(4):637-6) 671-6277-... 2020.
[68] Y. M. Tsai, T. Cojean, and H. Anzt. Sparse linear algebra AMD and NVIDIA GPUs the race is on. In P. Sadayappan, B. L. Chamberlain, G. Juckeland, and H. Ltaief, editors, High Performance Computing, pages 309-327, Cham, 2020. Springer International Publishing.
[69] R. J. Vanderbei et al. Linear programming, volume 3. Springer, 2015.
[70] S. Wang and N. Shroff. A new alternating direction method for linear programming. In Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
[71] Yang and Q. Lin. RSG: Beating subgradient method without smoothness and strong convexity. The Journal of Machine Learning Research, 19(1):236-266 2018.
[72] M. Zhu and Chan. An efficient primal-dual hybrid gradient algorithm for total variation image restoration. UCLA CAM Report, 34:8-34, 2008.
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] Our comparisons with basclines show that PDLP is not always the best, hence demonstrating its limitations. We also note which of our heuristics are lacking theoretical guarantees.
(c) Did you discuss any potential negative societal impacts of your work? [No] As a purely algorithmic paper, we do not believe such a discussion relevant. Linear programming a mature area whose societal impact is well understood.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yss]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] The only theoretical result appear in the appendix (Proposition 1).
(b) Did you include complete proofs of all theoretical results? [Yes]
3. If you experiments...
(a) Did you include the code, data, and instructions needed reproduce the main ex- perimental results (either in the supplemental material or URL)? [Yes] Code, data, and instructions needed to reproduce the main experimental results are publicly released in open source repository https: : h/tgttttttt h/ggtttt.b.. com/ google-resese FirstOrderLp. j1. Additional instructions are included the supplementary materi- als.
(b) Did you specify the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] We discussed that used MIP Relaxations to build our algorithm
14

---

and LP benchmark and Netlib as held out evaluation sets. We rum a large ablation study justify algorithmic decisions and most hyperparameter settings.
(c) Did you report error bars (e.g.. with respect to the random seed after running experi- ments multiple times)? [No] We used large datasets that would have been computa- tional intensive to run multiple times. When possible, we use the "KKT pass" metric that's reproducible and not subject measurement noise. Only the PageRank instances use random seed, and these are too large to run multiple times because we have single concurrent license for Gurobi.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, cloud provider)? [Yes] See Section 4.1.
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] We use LP Benchmark, MIP Relaxations, and Netlib datasets. We cite the sources for the datasets.
(b) Did you mention the license of the assets? [No] Although these datasets are widely used the community, we were unable find explicit licenses.
(c) Did you include any new assets either in the supplemental material or as URL? [Yes] We include our code and data processing scripts in the supplemental material.
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [No]
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid participants and the total amount spent on participant compensation? [N/A]
15


Reviews:
[{"review_id": "yuUKDwFm6Y1", "reviewer": "Reviewer_oPcr", "summary": "The paper applies the primal-dual hybrid-gradient (PDHG) method (also sometimes referred to as Chambolle-Pock algorithm) to\nsolve very large linear programs. The motivation is to solve problems which are out of reach for commerical LP solvers\nthat are based on simplex or interior-point approaches. Moreover, the algorithm is admissible to highly parallel computation due to the simple update rules of PDHG. \n\nThe main contribution is a careful combination of existing heuristic which are known to speed up PDHG in practice: restarts, \"line search\" and preconditioning.\n\nIn experiments, this combination of heuristics is shown to substantially improve a basic PDHG baseline. Contrary to existing wisdom,\nit is shown that FOMs, when enhanced with good heuristics, can indeed find highly accurate solutions to large-scale problems. \n", "questions": "", "limitations": "", "rating": 8, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "First of all, I found this paper to be very well written and enjoyable to read. The main strength of the paper is\nthe quality and clarity of the writing and the experimental setup. It is convincingly demonstrated\nthat the combination of proposed heuristics provide a strong improvement over the baseline method and existing solvers like SCS.\n\nThe introduction gives a beautiful overview and motivation of linear programming with first-order methods.\nI fully agree with the points raised in the extensive and complete literature review. Also from my experience with these methods,\nPDHG typically tends to perform much better than \"matrix-free\" implementations of ADMM or DRS.\n\nA comprehensive comparison to other FOMs for LPs (eg. [21,28,41,63,64] mentioned in the introduction) would have been nice to see, but likely\nout of scope for this work given the already very extensive experiments. Still, it would be nice to know whether PDLP is expected to perform\nbetter than, e.g., the ECLIPSE method.\n\nHowever, I see this more of a problem with the field of first-order methods of convex optimization, where countless algorithmic variants and tweaks\nare proposed but not properly evaluated on challenging standard benchmarks. This paper performs a strong experimental evaluation, which can serve\nas a baseline for further algorithmic developments. \n\nIn the introduction, the paper mentions that FOMs for LPs might not have a \"tailing\" effect due to the linear convergence rate. I found this to be\na bit misleading. To best of my knowledge, such convergence rates are all based on Hoffmann's constant which is a very loose worst case bound which\ndoes not coincide with the practical performance at all. Rather, the fast convergence could be due to a quick identification\nof the active set of constraints at the solution. Once the active set is identified, typically much faster linear convergence rates\nthan Hoffmann's bound apply. \n\nAfter reading the paper, I was wondering why PDHG was considered, since there are many other matrix-free first order methods. Was it simply the one\nthat performed best in practice or is there something special about it? I feel that the fast gradient method (with restarting) on the homogeneous self-dual embedding as considered in the paper \nby Neocara, Nesterov, Glineur, \"Linear convergence of first order methods for non-strongly convex optimization\", 2019, Sec 6.4 with similar enhancements\ncould also be strong candidate for linear programming. \n\nFor the experiments, the number of KKT passes is limited to 100.000 / the solve time to 1 hour. Was there a specific reason for this limitation?\nCan more problems in the MIP library be solved by increasing this limit, by say, factor 10x, or are there some problems which are completely \"out of reach\" for first-order methods? \n\nOne might argue that the weak point of the paper is perhaps its originality, since most considered heuristics are already known. \nBut there are some novel ideas, for example the considered combination of Ruiz+Pock/Chambolle diagonal preconditioning, which \nhas a surprisingly strong effect on the performance in my opinion. \n\nThis is a matter of taste, but as mentioned, I would have liked to see a longer concluding discussion on what have we learned about\nlinear programming or PDHG/first-order methods from this study. Does the paper generate any new insights or shed new light on the nature of PDHG?\nDoes it the paper open up directions for theoreticians to work on?\n\nFor example, as shown in experiments, preconditioning of PDHG often drastically improves the\nperformance. But preconditioning for PDHG is far less understood than preconditioning for linear systems. For linear systems, a direct connection\nof the condition number to the convergence speed can be made. But for nonlinear problems like linear-programming,\nwhat is a good notion of \"condition number\"?  Can we develop a theory which guides\nthe design of better preconditioners, similar to the theory of preconditioning for linear systems? \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The limitations of the approach are clearly mentioned throughout the paper; it is acknowledged that there are no convergence guarantees\nand that the paper presents mainly a heuristic approach which works well in practice.\n\nI don't believe that the paper needs to address possible negative societal impacts. \n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5 hours"}, "extra_text": {}}, {"review_id": "p84vJeNhSCb", "reviewer": "Reviewer_uQiX", "summary": "The authors investigate several heuristics to accelerate the PDHG algorithm, when applied to solve large-scale linear programs.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "1) My main concern is that there is no formal proof of convergence, let alone linear convergence. linear convergence of the PDHG for linear programs (and a wide range of other problems with metric subregularity) has been shown in\nAlacaoglu et al. \"On the convergence of stochastic primal-dual hybrid gradient\". So, you should explain for each heuristic whether it keeps or not the convergence guarantees. Preconditioning preserves convergence [53], but I don't know for the other heuristics.\n2) Your proposed line search is said to be faster that the one in [37], but there is no formal statement that is is valid.\n3) This is a paper on convex optimization, the relationship to machine learning is not clear enough to me to justify publication in NeurIPS. A journal on optimization would be more appropriate.\n4) I am not sure, but I think that the stochastic PDHG in \"Stochastic Primal-Dual Hybrid Gradient Algorithm with Arbitrary Sampling and Imaging Applications\" could be applied to LPs, using at every iteration activation of a random row of A or G. I would be much more convinced by such a study than by the present one, in the context of NeurIPS. \n\nA remark: PDHG converges when tau.sigma||K||^2 <= 1, not <1, as proved in [18]. And fastest convergence is obtained in practice when tau.sigma||K||^2 = 1, so you should set eta=1/||K||.\n\n=== update after author response period\nI read all reviews and corresponding authors' responses, and I think that the authors did a very good job in responding, in a professional way, to all concerns raised.\nYou convinced me about the value of the paper for the practitioners willing to solve large LPs in practice. So, I revise my judgement about Neurips not being the most appropriate venue for the paper, and I recognize that the contributions are interesting for the Neurips community. Also, your reply to my point 4. is well argued. Therefore, I am increasing my score from 5 to 6. With some formal statements about convergence, I would have put an even higher score.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Not applicable", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "1.5"}, "extra_text": {}}, {"review_id": "Ze8qOEL7rrB", "reviewer": "Reviewer_eniT", "summary": "This paper presents PDLP, a practical first-order method for linear programming (LP) that can solve to the high levels of accuracy that are expected in traditional LP applications. The proposed algorithm is matrix-free and can scale to very large problems. The core idea is based on the saddle-point reformulation of LP and the celebrated primal-dual hybrid gradient (PDHG) method (cf. Chambolle and Pock, 2011). Combined with several new techniques, e.g., diagonal preconditioning, presolving, adaptive stepsizes and adaptive restarting, PDLP improves the state-of-the-art approaches, including SCS, by a significant margin. Experimental results on PageRank applications are encouraging. \n", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Pros: As far as I could see, the mathematics appears to be intuitive, theoretically sound, although I did not check the details thoroughly. Despite the lack of the convergence results for several heuristics, the extensive numerical experiments were presented and supported the practical benefits of this new approach. That's the thing that I liked most about this paper since it aims at some useful approaches in real scenario. The strategies for tuning stepsizes as well as the restart criteria and diagonal preconditioning make a lot of sense to me and I appreciate that the authors put effort to explain why these techniques are useful from a practical perspective. This is an important contribution of the paper.\n\nCons: One of the reasons that I did not give this paper a higher score is that some relevant literatures are missing. For example, due to a very specific structure of LP, some specialized second-order algorithms exhibit superior practical performance in real applications. (also with theoretical guarantee)\n\nX. Li, D. Sun, K. C. Toh. \"An asymptotically superlinearly convergent semismooth Newton augmented Lagrangian method for Linear Programming\", SIOPT, 2020. \n\nA. Ali, E. Wong, J. Z. Kolter, \"A semismooth Newton method for fast, generic convex programming\", ICML 2017. \n\nThe presolving and diagonal preconditioning techniques are also combined with some first-order interior-point methods for LP. The resulting approach outperforms SCS on NETLIB collections and some other ML datasets thanks to these heuristics. \n\nT. Lin, S. Ma, Y. Ye, S. Zhang, \"An ADMM-based interior-point method for large-scale linear programming\", Optimization Methods and Software, 2021. \n\nMoreover, on account of the unique role of NETLIB collections in LP history, it would be better to evaluate PDLP on some of these datasets and see if this approach is sufficient robust in practice. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "1. I encourage the authors to do further literature review and include other recent methods for LP, especially FOMs. \n\n2. It would be better to evaluate PDLP on NETLIB collections and see if it still behaves well on hard LP instances. \n\n3. There are a couple of typos in the paper and the descriptions and language could be clearer in places.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2 hours"}, "extra_text": {}}, {"review_id": "JitxxeDriN5", "reviewer": "Reviewer_nQvW", "summary": "This paper proposes a heuristic PDLP for large-scale LP by combining some existing techniques into PDHG. Emprical results are good, but no convergence guarantees are provided. ", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper makes beneficial attempts for applying first-order methods into LP. But I found that all the tricks and heuristics in this paper are already known in existing papers --- the algorithm in this paper just simply combines these parts. As a result, this paper is more like instructions of a software package, not like an academic paper. For engineering implementation, this is a good guideline; but for academic research, this is not enough. \n\n\nSome other comments:\n\nIn my understanding, a linear convergence result is not enough to prevent the tailing-off effect, because the condition number can be arbitrarily large so contracting by 1/2 is still very expensive. I believe that the reason why we say IPMs can provide high-accurate solution is they have superlinear convergence at the end of the stage which is independent of the condition number and make us get an almost exact solution in finite steps. \n\nAlgorithm 2 is the main part for the algorithm, but the name is \"step size heuristic\", which is not proper. \n\nIn the line 133-134, the two uses of the notation $\\rho_{\\cdot}^{\\cdot}$ is not consistent. The first subscript is a constant, while the second one is a variable. \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "No Societal Impact", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}]

Author rebuttal:

Meta-review: Thank you for your submission to NeurIPS. All four reviewers agree that the paper is valuable because it solves an important practical problem well. Initially, there were some concerns regarding the novelty of the methods, and the heuristic nature of the methods. Fortunately, most of these concerns were adequately addressed after discussion with the authors. Please ensure that the camera-ready version is adequately revised to capture the points raised during the discussion period, as these were important in helping the reviewers come to their final decisions. 
Decision: Accept (Poster)