# Scylla: a matrix-free fix-propagate-and-project heuristic for mixed-integer optimization

Gioni Mexi, Mathieu Besançon, Suresh Bolusani, Antonia Chmiela, Alexander Hoen, Ambros Gleixner

► To cite this version:

Gioni Mexi, Mathieu Besançon, Suresh Bolusani, Antonia Chmiela, Alexander Hoen, et al.. Scylla: a matrix-free fix-propagate-and-project heuristic for mixed-integer optimization. Operations Research Proceedings 2023, Sep 2023, Hamburg, Germany. pp.65-72, 10.1007/978-3-031-58405-3_9. hal-05216987

HAL Id: hal-05216987

https://hal.science/hal-05216987v1

Submitted on 20 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Scylla: a matrix-free fix-propagate-and-project heuristic for mixed-integer optimization

Gioni Mexi*¹[0000-0003-0964-9802], Mathieu Besançon¹[0000-0002-6284-3033], Suresh Bolusani¹[0000-0002-5735-3443], Antonia Chmiela¹[0000-0002-4809-2958], Alexander Hoen¹[0000-0003-1065-1651], and Ambros Gleixner¹,²[0000-0003-0391-5903]

¹ Interactive Optimization and Learning, Zuse Institute Berlin, Germany

² HTW Berlin, Germany

{mexi, besancon, bolusani, chmiela, hoen, gleixner}@zib.de

**Abstract.** We introduce Scylla, a primal heuristic for mixed-integer optimization problems. It exploits approximate solves of the Linear Programming relaxations through the matrix-free Primal-Dual Hybrid Gradient algorithm with specialized termination criteria, and derives integer-feasible solutions via fix-and-propagate procedures and feasibility-pump-like updates to the objective function. Computational experiments show that the method is particularly suited to instances with hard linear relaxations.

**Keywords:** Mixed-integer optimization · heuristics · matrix-free.

## 1 Introduction

Mixed-integer optimization (MIP) is an essential and challenging class of problems found in numerous applications including transport, power systems, engineering, or manufacturing. Primal heuristics are used as stand-alone methods to compute good-quality solutions or as subroutines in branch-and-bound solvers, potentially exploiting the information extracted in the rest of the solution process [5]. Solutions derived from primal heuristics can be returned immediately when the corresponding gap is small, they help to prune nodes of the branch-and-bound tree, and are also leveraged by other subroutines of the solver, e.g., cut selection [12].

When integrated within a branch-and-bound scheme, many primal heuristics assume that an optimal solution of the linear programming (LP) relaxation is readily available and that the LP can be re-optimized quickly after few modifications. These repeated calls to LP solvers can become a bottleneck in several scenarios. First, while it is true that standard MIP solvers produce and rely on optimal LP relaxation solutions, they are not yet available during presolving. Second, for some applications, re-optimization via the simplex method is expensive [14]. For primal heuristics in a stand-alone setting, the computation of an

<footnote>* Corresponding author.</footnote>
initial solution to the LP relaxation by standard methods that require matrix factorizations, such as the simplex or interior point methods, may be expensive for large-scale applications.

To address this problem, we present *Scylla*, a *matrix-free* primal heuristic, in the sense that it uses only iterated matrix-vector multiplications to compute approximate LP solutions without factorizing the constraint matrix nor directly solving linear systems. We supplement this with LP-free constraint propagation to round such approximate reference solutions to integer feasibility and objective function updates inspired by the well-known feasibility pump [9,4,2].

The paper is organized as follows. In Section 2, we present the algorithm design and highlight its key features. In Section 3, we report on our computational experiments. We summarize our observations and discuss future research directions in Section 4.

## 2 The algorithm

The overall structure of Scylla is outlined in Algorithm 1.1 and is based upon several subroutines:

- a limited-precision call to the Primal-Dual Hybrid Gradient (PDHG) algorithm [17,8] to compute an approximate solution of the LP relaxation,

- a fix-and-propagate procedure to attempt rounding this fractional solution to a fully, i.e., also integer, feasible solution,

- an objective function update similar to that of the feasibility pump [9,4,2] in order to iteratively drive the PDHG solution towards integer feasibility.

The outer loop of the algorithm is relatable to a feasibility pump [9,4] where the rounding step uses constraint propagation [10] and the projection step is replaced with a limited number of PDHG iterations until relaxed termination criteria are met. The inner iterations counted by $k$ are dominated by matrix-vector multiplications with the constraint matrix and its transpose; the outer iterations counted by $K$ are dominated by a fix-and-propagate call and update the objective function at the end. Note that before starting the actual algorithm, we presolve the instance to reduce its size using the SCIP solver [7].

In the following sections, we describe the principal components of the algorithm and their interaction.

### 2.1 Approximate LP solutions by matrix-free linear optimization

Heuristic and exact methods for mixed-integer problems often rely on repeated calls to a linear optimization solver. These calls represent a substantial amount of the runtime and require, for the most common LP solving techniques, a matrix factorization. As highlighted in [3], first-order methods applied to the saddle-point formulations of LPs open the door to a trade-off between solution accuracy and runtime which was rarely reasoned about in the linear and mixed-integer
<table><caption>**Algorithm 1.1:** High-level structure of the Scylla heuristic</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>Input</td><td colspan="2">: MIP instance min{c<sup>T</sup>x : Ax &ge; b, l &le; x &le; u, x<sub>j</sub> &isin; Z, &forall; j &isin; I}</td></tr><tr><td>Output</td><td colspan="2">: MIP-feasible solution x or &emptyset; if none found</td></tr><tr><td>1</td><td colspan="2">initialize k = 1, K = 1, &alpha; = 0.9</td></tr><tr><td>2</td><td colspan="2">initialize primal-dual step sizes &tau;<sub>x</sub> &gt; 0, &tau;<sub>y</sub> &gt; 0 and x&#773;<sup>(k)</sup> = 0, y&#773;<sup>(k)</sup> = 0</td></tr><tr><td>3</td><td colspan="2">initialize &hat;c = c</td></tr><tr><td>4</td><td colspan="2">while limits not reached do</td></tr><tr><td>5</td><td colspan="2">repeat</td></tr><tr><td>6</td><td colspan="2">update step sizes &tau;<sub>x</sub>, &tau;<sub>y</sub></td></tr><tr><td>7</td><td colspan="2">d<sub>x</sub> = A<sup>T</sup>y&#773;<sup>(k)</sup> &minus; &hat;c and x&#773;<sup>(k+1)</sup> = proj<sub>[l,u]</sub>(x&#773;<sup>(k)</sup> + &tau;<sub>x</sub>d<sub>x</sub>)</td></tr><tr><td>8</td><td colspan="2">d<sub>y</sub> = b &minus; A(2x&#773;<sup>(k+1)</sup> &minus; x&#773;<sup>(k)</sup>) and y&#773;<sup>(k+1)</sup> = proj<sub>R<sub>+</sub><sup>n</sup></sub>(y&#773;<sup>(k)</sup> + &tau;<sub>y</sub>d<sub>y</sub>)</td></tr><tr><td>9</td><td colspan="2">k = k + 1</td></tr><tr><td>10</td><td colspan="2">until x&#773;<sup>(k)</sup>, y&#773;<sup>(k)</sup> &asymp; LP optimal or bounds stalling or k &gt; iterlimit</td></tr><tr><td>11</td><td colspan="2">if x&#773;<sup>(k)</sup> MIP-feasible then return x&#773;<sup>(k)</sup></td></tr><tr><td>12</td><td colspan="2">x&#775;<sup>(K)</sup> = fix-and-propagate(x&#773;<sup>(k)</sup>)</td></tr><tr><td>13</td><td colspan="2">if x&#775;<sup>(K)</sup> MIP-feasible then return x&#775;<sup>(K)</sup></td></tr><tr><td>14</td><td colspan="2">if cycling on x&#775;<sup>(K)</sup> detected then x&#775;<sup>(K)</sup> = perturb(x&#775;<sup>(K)</sup>)</td></tr><tr><td>15</td><td colspan="2">&hat;c<sup>T</sup>x = &alpha;<sup>K</sup> &radic;(|I|)/||c|| c<sup>T</sup>x + (1 &minus; &alpha;<sup>K</sup>) &Delta;(x, x&#775;<sup>(K)</sup>)</td></tr><tr><td>16</td><td colspan="2">K = K + 1</td></tr><tr><td>17</td><td colspan="2">return &emptyset;</td></tr></tbody></table>

setting. In our algorithm, these approximate LP solutions are used as a starting point for a fix-and-propagate procedure.

To leverage this trade-off, we use the PDHG [8] algorithm and its recent PDLP implementation (PDHG specialized for linear optimization) [3]. Some notable algorithmic improvements were introduced in [3] compared to the standard PDHG including adaptive step size strategies, restarts, and matrix preconditioning. We use the PDLP implementation available in OR-Tools<fnref n="1" />. The modifications allow us to warm-start the solution process from a previously computed primal-dual pair and to add custom stopping criteria detailed in Section 2.2.

For conciseness of the algorithm description, we only included the core PDHG iterations and did not detail the step size computation procedure or restarting scheme from PDLP. For Lines 5 to 10 of Algorithm 1.1, Scylla calls PDLP to solve similar LPs where only the objective function vector changes. We use the fact that the state of PDLP is fully described by its primal-dual pair and primal and dual residuals to warm-start all calls after the first one, which may help to reduce the total number of iterations after the first solve. Using PDLP also means that Scylla remains a matrix-free heuristic only requiring sparse matrix-vector multiplications.

<footnote n="1">¹ The package is available at https://github.com/google/or-tools.</footnote>
## 2.2 Custom stopping criteria

The standard stopping criterion for PDLP is a maximum error $\epsilon$ on the primal and dual feasibilities. This error can be relaxed compared to what would be expected in “standard” LP solving since the resulting fractional solution is then guiding a rounding procedure. In particular, this error can also be adapted dynamically to allow a bigger tolerance at the beginning of the heuristic. In addition to increasing the initial error tolerance (to $\epsilon^{(1)} = 0.01$ in our experiments), we also use a progressive refinement scheme in which the error is reduced geometrically by a factor $\beta = 0.98$, resulting in the following.

$$ \epsilon^{(K+1)} = \max\{\beta \epsilon^{(K)}, 10^{-8}\} $$

## 2.3 The fix-and-propagate procedure

Fix-and-propagate algorithms have been used successfully as part of branch-and-bound solvers [6,11]. They are also an essential step in the most successful variants of feasibility pump algorithms [10], and we use them with the same motivation: to “round” the PDHG solution $\bar{x}^{(k)}$ to integer feasibility.

“Rounding” is achieved by fixing an integer variable with fractional value to an integer within its domain, followed by a propagating step. This is repeated until all such variables are fixed or infeasibility is detected by some domain becoming empty. In the latter case, fix-and-propagate continues in order to produce an integer vector by ignoring any constraint that would lead to empty domains. At the end of the propagation, all remaining unfixed variables are fixed to their values in the fractional reference solution or its projection on to their domain. The procedure always produces an integer-feasible, but not necessarily LP-feasible, solution.

One major issue of feasibility pump-like algorithms and hence also Scylla is cycling: after some iterations, it may happen that $\hat{x}^{(K)} = \hat{x}^{(K')} \text{ for } 1 \le K' < K$. Similar to [9], this issue is handled by randomly perturbing some of the variables in the current rounded vector $\hat{x}^{(K)}$ (Line 14).

## 2.4 Objective function update

Inspired by objective feasibility pump approaches [9,4], we update the objective function as follows. After computing the integer point $\hat{x}^{(K)}$ and before the subsequent PDLP call, the objective function is updated (Line 15 of Algorithm 1.1) as a convex combination of the initial objective function $c^\top x$ and an $\ell_1$-norm distance to the integer vector $\Delta(x, \hat{x}^{(K)})$ with $\alpha \in [0, 1)$ chosen as constant 0.9 in our experiments.

This convex combination represents a trade-off between the initial objective (towards high-quality solutions) and the distance to an integer vector (towards integer feasibility). After each iteration $K$, the weight shifts further towards integer feasibility in the attempt to produce a PDLP solution $\bar{x}^{(k)}$ that fix-and-propagate can successfully convert to a MIP-feasible solution.
## 3 Computational results

As a baseline for our implementation, we modified an open source implementation of the Feasibility Pump [16] (FP)<fnref n="2" /> (that uses SoPlex 6.0.2 [7] as the underlying LP solver) and extended it to support approximate LP solves with PDLP.

Our test set consists of the benchmark library MIPLIB 2017 [13] and the instances used in [10]. We excluded all instances that were solved to optimality during presolving with SCIP, leaving us with 283 instances.

All results were obtained on a cluster of Intel Xeon Gold 6338 @ 2.0 GHz machines with 1024 GB of RAM and a time limit of 3600 seconds. To mitigate the impact of performance variability [15], both FP and Scylla are tested on the complete test set with ten different random seeds.

Table 1 shows the summarized results for 2830 instances (283 instances × 10 random seeds). On the whole testset, FP outperforms Scylla w.r.t. both solution quality and running time. This is an expected result since the LPs in the FP iterations are for most instances easily solved to optimality, and warm-starting of simplex LP solvers is very effective. Scylla therefore loses the benefit of exact LP solutions while still requiring similar costs per iteration. We see this when looking at the two different subsets of instances.

For the instances where the root LP solves in [1, 60] seconds, FP is able to find more solutions than Scylla and provides better quality solutions. Scylla is however faster on average, despite suffering under a greater model construction overhead which increases with the instance size. Note that this overhead is tied to the current implementation and not to the heuristic itself. More specifically, the current OR-Tools interface through which PDLP is accessed requires re-building the model at every iteration of the feasibility pump, a restriction that penalized our approach, in particular for larger instances.

For instances in [60, lim], Scylla significantly improves over FP in terms of both the number of solved instances and running time. Notably, Scylla is capable of finding feasible solutions for 42 more instances than FP, and it is nearly four times faster. These results suggest that for MIPs with difficult LP relaxations, Scylla is the preferred heuristic choice to compute first primal solutions.

## 4 Conclusion

In this paper, we presented Scylla, a primal heuristic for mixed-integer linear optimization that combines matrix-free computation of approximate LP solutions with a fix-and-propagate procedure and an objective function update scheme similar to the feasibility pump. Computational experiments on a large set of instances from the MIPLIB 2017 collection show that the method is particularly suited to instances with hard linear relaxations. Our current implementation does not exploit the potential for parallelization that lies in the sparse matrix-vector multiplications of PDHG, which dominate the PDLP runtime. A robust,

<footnote n="2">² https://github.com/GioniMexi/FeasPumpCollection</footnote>
<table><caption>**Table 1.** Aggregated results comparing FP with Scylla. Row [1, 60] ([60, lim]) only considers instances on which solving the LP takes between 1 and 60 seconds (more than 60 seconds). Column ‘#’ (‘only’) refers to the number of instances for which a solution was found (exclusively by only one heuristic), ‘gap’ refers to the relative gap to the best known solution, ‘wins’ refers to the number of instances with a performance improvement of 10% over the other heuristic, and ‘overhead’ refers to the overhead time to build the model.</caption><thead><tr><th rowspan="2">testset</th><th rowspan="2">algorithm</th><th colspan="2">found sol</th><th colspan="2">#wins</th><th colspan="2">shift. geom. mean</th><th rowspan="2">overhead(s)</th></tr><tr><th>#</th><th>only</th><th>gap</th><th>time(s)</th><th>gap</th><th>time(s)</th></tr></thead><tbody><tr><td rowspan="2">all</td><td>FP</td><td>2268</td><td>152</td><td>1154</td><td>1557</td><td>27.46</td><td>10.06</td><td>0.02</td></tr><tr><td>Scylla</td><td>2270</td><td>154</td><td>485</td><td>797</td><td>40.92</td><td>17.30</td><td>1.93</td></tr><tr><td rowspan="2">[1, 60]</td><td>FP</td><td>660</td><td>62</td><td>341</td><td>295</td><td>26.75</td><td>43.16</td><td>0.01</td></tr><tr><td>Scylla</td><td>635</td><td>37</td><td>130</td><td>381</td><td>41.78</td><td>36.75</td><td>1.99</td></tr><tr><td rowspan="2">[60, lim]</td><td>FP</td><td>163</td><td>40</td><td>98</td><td>64</td><td>28.11</td><td>974.66</td><td>0.01</td></tr><tr><td>Scylla</td><td>205</td><td>82</td><td>105</td><td>180</td><td>35.37</td><td>234.80</td><td>6.09</td></tr></tbody></table>

efficient, and multithreaded sparse BLAS library will thus accelerate this core operation of PDLP and could be key to a broader success especially on large-scale instances.

One aspect to be explored in future research is the fact that PDHG is by design suited to solve more general convex differentiable optimization problems, and PDLP handles quadratic problems where the objective Hessian is diagonal [1]. On the one hand, this makes Scylla a candidate heuristic for the broader class of MIQPs. On the other hand, this can be exploited also in the linear setting by choosing the $l_2$-norm as the distance function $\Delta(\cdot, \hat{x}^{(K)})$ and avoiding auxiliary variables needed to model the $l_1$-norm on general integer variables [9,4].

## Acknowledgements

Research reported in this paper was partially supported through the Research Campus Modal funded by the German Federal Ministry of Education and Research (fund numbers 05M14ZAM,05M20ZBM) and the Deutsche Forschungsgemeinschaft (DFG) through the DFG Cluster of Excellence MATH+.

## References

1. Mathematical background for PDLP, Google OR-tools documentation. https://developers.google.com/optimization/lp/pdlp.math, accessed: 2023-06-27

2. Achterberg, T., Berthold, T.: Improving the feasibility pump. Discrete Optimization 4(1), 77–86 (2007)

3. Applegate, D., Díaz, M., Hinder, O., Lu, H., Lubin, M., O’Donoghue, B., Schudy, W.: Practical large-scale linear programming using primal-dual hybrid gradient. Advances in Neural Information Processing Systems 34 (2021)
4. Bertacco, L., Fischetti, M., Lodi, A.: A feasibility pump heuristic for general mixed-integer problems. Discrete Optimization **4**(1), 63–76 (2007)

5. Berthold, T.: Heuristic algorithms in global MINLP solvers. Ph.D. thesis, Technische Universität Berlin (2014)

6. Berthold, T., Hendel, G.: Shift-and-propagate. Journal of Heuristics **21**(1), 73 – 106 (2014). https://doi.org/10.1007/s10732-014-9271-0

7. Bestuzheva, K., Besançon, M., Chen, W.K., Chmiela, A., Donkiewicz, T., van Doornmalen, J., Eifler, L., Gaul, O., Gamrath, G., Gleixner, A., et al.: Enabling research through the scip optimization suite 8.0. ACM Transactions on Mathematical Software (2023)

8. Chambolle, A., Pock, T.: A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of mathematical imaging and vision **40**(1), 120–145 (2011)

9. Fischetti, M., Glover, F., Lodi, A.: The feasibility pump. Mathematical Programming **104**(1), 91–104 (2005)

10. Fischetti, M., Salvagnin, D.: Feasibility pump 2.0. Mathematical Programming Computation **1**(2), 201–222 (2009)

11. Gamrath, G., Berthold, T., Heinz, S., Winkler, M.: Structure-driven fix-and-propagate heuristics for mixed integer programming. Mathematical Programming Computation **11**(4), 675–702 (2019)

12. Gleixner, A., Bastubbe, M., Eifler, L., Gally, T., Gamrath, G., Gottwald, R.L., Hendel, G., Hojny, C., Koch, T., Lübbecke, M.E., Maher, S.J., Miltenberger, M., Müller, B., Pfetsch, M.E., Puchert, C., Rehfeldt, D., Schlösser, F., Schubert, C., Serrano, F., Shinano, Y., Viernickel, J.M., Walter, M., Wegscheider, F., Witt, J.T., Witzig, J.: The SCIP Optimization Suite 6.0. Tech. rep., Optimization Online (2018)

13. Gleixner, A., Hendel, G., Gamrath, G., Achterberg, T., Bastubbe, M., Berthold, T., Christophel, P., Jarck, K., Koch, T., Linderoth, J., et al.: MIPLIB 2017: data-driven compilation of the 6th mixed-integer programming library. Mathematical Programming Computation pp. 1–48 (2021). https://doi.org/10.1007/s12532-020-00194-3

14. Koch, T., Achterberg, T., Andersen, E., Bastert, O., Berthold, T., Bixby, R.E., Danna, E., Gamrath, G., Gleixner, A.M., Heinz, S., Lodi, A., Mittelmann, H., Ralphs, T., Salvagnin, D., Steffy, D.E., Wolter, K.: MIPLIB 2010. Mathematical Programming Computation **3**(2), 103–163 (2011). https://doi.org/10.1007/s12532-011-0025-9

15. Lodi, A., Tramontani, A.: Performance variability in mixed-integer programming. In: Theory Driven by Influential Applications, pp. 1–12. INFORMS (2013). https://doi.org/10.1287/educ.2013.0112

16. Mexi, G., Berthold, T., Salvagnin, D.: Using multiple reference vectors and objective scaling in the feasibility pump. EURO Journal on Computational Optimization **11**, 100066 (2023)

17. Zhu, M., Chan, T.: An efficient primal-dual hybrid gradient algorithm for total variation image restoration. UCLA Cam Report **34**, 8–34 (2008)