A New Crossover Algorithm for LP Inspired by the Spiral Dynamic of PDHG

PDHG was studied in Chambolle and Pock (2011) for image processing applications. Applegate et al. (2021) developed a restarted PDHG algorithm for solving linear programs (LPs), which, together with a few practical enhancements, leads to the PDLP solver. Compared with other popular FOM solvers such as SCS and ABIP, there is no linear system to solve in PDLP, making it a more suitable choice for large-scale LPs. cuPDLP, the GPU implementation of PDLP, demonstrated strong numerical performance on par with commercial LP solvers, and has attracted attention from both the academic and solver industry (Lu and Yang, 2023a; Lu et al., 2023). Recently, additional numerical enhancements have been proposed to speed up the algorithm further. Xiong and Freund (2024) pointed out that the central-path-based scaling can markedly improve the convergence rate of PDHG. Lu and Yang (2024a) embedded Halpern iteration (Halpern, 1967) in PDLP, achieving accelerated theoretical guarantees and improved computational results.

The practical success of PDLP motivates a stream of theoretical analysis. Lu and Yang (2022) showed PDHG has a linear convergence rate when solving LP. Applegate et al. (2023) revealed that an accelerated linear convergence rate can be achieved with restarts, and such an accelerated linear rate matches the lower bound. Applegate et al. (2024) illustrated that for infeasible/unbounded LP, PDHG iterates diverge towards the direction of infeasibility/unboundedness certificate; thus, one can detect the infeasibility/unboundedness of LP using PDHG without additional effort. Lu and Yang (2024b) discovered that PDHG has a two-stage behavior when solving LP. In the first stage, PDHG uses finite iterations to verify the optimal active variables; in the second stage, PDHG solves a homogeneous linear inequality system with a significantly improved linear rate. This explains the slow behaviors of PDHG when the LP is near-degeneracy. Xiong and Freund (2024) introduced level-set geometry that characterizes the difficulty of PDHG for LP, which shows that PDHG converges faster for LP with a regular non-flat level-set around the optimal solution set.

For more general fixed point problems, trajectories of various FOMs were studied in Poon and Liang (2019, 2020), which showed that eventually (after a finite number of iterations), these algorithms follow a straight line or a spiral structure if the problem is non-degenerate. Unfortunately, this non-degeneracy assumption is generally never satisfied for real-world LP instances.

Crossover, or purification, refers to the process that moves a feasible or optimal solution to a vertex solution, also known as the extreme point, basic feasible solution, or corner solution. Purification methods (Charnes and Kortanek, 1963; Charnes et al., 1965) were proposed even earlier than those non-vertex LP methods, including the ellipsoid methods (Khachiyan, 1979) and IPMs (Karmarkar, 1984).

Given a basis $B$ and a feasible solution $\mathbf{x}$, the purification algorithm (Charnes and Kortanek, 1963; Charnes et al., 1965) generates a null space vector step $\bm{\alpha}$ same as the direction in simplex methods and then tries to move $\mathbf{x}$ along it without worsening the objective or update $B$ for feasible directions, which will finally return a feasible vertex with the objective value no worse in finite steps. Given a basis $B$ and a strictly interior feasible point, Kortanek and Jishan (1988) improved Charnes’s method by drawing the vertex back to the interior if it is not optimal. An optimal extreme point will be reached in a finite number of steps. If the feasible point is not entirely in the interior, Kortanek and Jishan (1988) also provided another purification method that considers dual feasibility. Combining Kortanek’s second method with a proper cleanup will return an optimal basis. Given a dual optimal solution $\bm{\pi}$, Amor et al. (2006) constructed an auxiliary dual LP by adding box constraints around $\bm{\pi}$ to the original dual problem. Dualizing back to obtain an auxiliary primal LP, Amor claimed that the auxiliary primal LP can be efficiently solved by simplex methods, and then an optimal vertex can be easily recovered.

Starting from an optimal primal-dual solution pair, Megiddo (1991) designed a strongly polynomial crossover approach to an optimal vertex. Megiddo also pointed out that given only the optimal primal or dual solution, solving a general LP is theoretically as hard as solving it from scratch. Megiddo’s crossover method lays down the algorithmic and theoretical foundation for modern crossover algorithms in modern LP solvers. Our randomized crossover modifies Megiddo’s crossover framework, which will be introduced in Section 3.

With the boom of various IPMs and MIP applications, crossover gradually plays a significant role in the LP solvers. The key task becomes refining IPM solutions to optimal vertices. Mehrotra (1991) added a controlled random cost perturbation to eliminate dual degeneracy at the sacrifice of losing only a little optimality. Especially for nondegenerate LPs, an indicator (Tapia and Zhang, 1991) can be calculated from the interior solution to identify the basis. Bixby and Saltzman (1994) implemented a modified simplex method, which accepts an interior feasible point as a super-basic solution. Then similar to a normal simplex method, Bixby’s algorithm uses the primal (or dual) ratio test to push primal (or dual) variables to bound if possible or pivots it into (or out of) the bound. Andersen and Ye (1996) used the strictly complementary property of IPM solutions to define a perturbed LP problem. Andersen’s method obtains a basis and then recovers feasibility by pivoting and taking a weighted average with the given IPM solution. Ge et al. (2021) also proposed a perturbation crossover method that combines variable fixing and random cost perturbation. They proved that as long as the perturbation is sufficiently small, their crossover approach will not bring objective loss.

Our crossover modifies Megiddo’s framework and adopts random perturbation, together with PDLP and ordinary least squares (OLS) solver as subroutines. More details will be discussed in Section 3.

Our primal push blends two approaches: solving one auxiliary LP and a few sequential OLS problems that identify the ray of PDHG iterates.

After recognizing the basic set $B$, we build an auxiliary LP with random cost perturbation $\tilde{\mathbf{c}}_{B}$,

and solve (16) with PDLP. This auxiliary LP has multiple optimal solutions w.p. 0 (see Lemma 1 in Ge et al. (2021)), so with proper perturbation to avoid the unbounded case, we have the unique optimal solution $\mathbf{x}_{B}$ to (16). Then $(\mathbf{x}_{B},\mathbf{0}_{E})$ can be the primal part of an optimal vertex of (1).

Unfortunately, due to the numerical residuals in practice, we may not correctly fix all the variables on bounds at once. To reduce the negative impact of the issue, the first thing to consider is to control the amplitude of disturbance, see Mehrotra (1991); Ge et al. (2021) for more specific discussions. Besides, a verification scheme needs to be designed.

Therefore, after solving (16), we solve a series of OLS problems (17) with various randomly generated costs $\tilde{\mathbf{c}}_{B}$ to calculate the PDHG primal ray directions (7).

From Proposition 10, since the current primal solution $\mathbf{x}_{B}$ is still feasible for (16) with cost perturbation, it is located at the primal spiral center of the rotation in the current phase of PDHG. We can then start from $\mathbf{x}_{B}$ and move along $(\bm{\delta}_{\mathbf{x}})_{B}$ until reaching new boundaries, mirroring how the PDHG solves the auxiliary LPs. At the end of the primal push, no more primal variable can be pushed to its bound when columns of $\mathbf{A}_{B}$ are linearly independent. OLS can be solved by the direct QR factorization, LSMR (Fong and Saunders, 2011), or the conjugate gradient (CG) method.

The reason why we only solve one auxiliary LP is that the auxiliary LP (16) is in most cases enough to obtain a primal solution close to a vertex and that the OLS has lower complexity to reduce the support and verify the termination of primal push. Directly following the primal spiral ray avoids unnecessary oscillation. However, when the original LP is too large, or the resulting OLS is too hard to solve, using multiple auxiliary LPs can be beneficial for the primal push.

In the dual push, we still combine the auxiliary LP and sequential OLS.

For the auxiliary LP, we sum up the non-negative dual slacks as a penalty to activate more dual constraints.

Then similarly, we compute a series of OLS (19) with perturbed right-hand side $\tilde{\mathbf{b}}$ to calculate the dual ray directions of PDHG (7).

More constraints are fixed to be active step by step until rows of $\mathbf{A}_{D}$ are linearly independent when $\mathbf{y}$ will be uniquely determined by $\mathbf{A}_{D}^{\top}\mathbf{y}=\mathbf{c}_{D}$.

After the primal and dual push, $(\mathbf{x},\mathbf{y})$ should reach some vertex. All we need is to identify $m-|B|$ columns from $\mathbf{A}_{C}=\mathbf{A}_{D\setminus B}$ to form an optimal basis with $\mathbf{A}_{B}$. Here we apply LU factorization twice. We first decompose $\mathbf{A}_{B}$ as

where $(R_{1},R_{2})$ divides the rows according to the decomposition. The basic candidates become

Then, we only need to select $m-|B|$ linearly independent columns from $\mathbf{A}_{R_{2},C}-\mathbf{L}_{R_{2},B}\mathbf{L}_{R_{1},B}^{-1}\mathbf{A}_{R_{1},C}$, where the second LU decomposition is involved.

So far, we have built our theory and algorithm design on the standard LP (1). In practice, the general LP formulation (20) is more convenient and thus widely used.

To extend the concept of the basis on the general form, we reformulate it as (21) via introducing the artificial variable $\mathbf{w}$.

where $\delta=10$ aims to avoid the right-hand-side term becoming $\mathbf{0}$, as the presence of $\mathbf{0}$ may render the LP highly sensitive to perturbations, causing the optimal solution to change drastically. This reformulation guarantees the constraint matrix $\begin{bmatrix}\mathbf{A}&-\mathbf{I}\end{bmatrix}$ to be full row rank, where we can easily apply the definition of the basic solution.

At the beginning of primal and dual push, we need to determine $B$ and $(D,N)$ to construct the auxiliary LP and OLS. For an IPM solution, which tends to stay away from bounds, the central path and strict complementarity are crucial for identifying non-basic variables. Identification is typically done by comparing the orders of magnitude of $\mathbf{x}$ and $\mathbf{c}-\mathbf{A}^{\top}\mathbf{y}$. In contrast, PDHG can benefit from the sparsity of its solution, induced by the projection, to more aggressively fix variables to their bounds.

We take

and after $B$ is updated in primal push

with $\gamma=1$ and $\epsilon=1e-8$. This criterion works well in our experiments for purifying the PDLP solutions.

The cost and right-hand-side perturbations guide the direction of crossover. They are perturbed as

where $\bm{\delta}_{\mathbf{c}},\bm{\delta}_{\mathbf{b}}$ are independently and uniformly distributed over the interval $[0,1]$. The first normalized terms function as a center to stabilize the disturbance, preventing the unfixed non-basic variables from occasionally escaping boundaries.

Currently, all the OLS problems are solved independently from scratch. Note that the coefficient matrices of two adjacent OLS problems differ in only a few columns or rows; thus, further development on reusing the factorization (probably via low-rank updates) may largely enhance the efficiency of our crossover to the level of the classic simplex-based crossover.

This crossover scheme can be applied to solutions obtained by IPMs or PDLP. However, their structures differ as discussed below.

IPMs find the analytical center of the optimal solution face, and the obtained solutions tend to be denser, while PDLP solutions are usually much sparser because of the projection operator.

IPMs and PDHG both have a strong ability to distinguish non-basic variables. IPMs can identify the strictly complementary solution pair with a sufficiently small duality gap (Andersen and Ye, 1996). PDHG can correctly classify variables as basic and non-basic in the scope of PDHG basis within finite iterations (Lu and Yang, 2024b).

Therefore, even though the PDLP solution may be less accurate than the IPM solution, it can still be refined to the vertex successfully in practice.