Linear Programming in MILP Solving

A Computational Perspective

Matthias Miltenberger

Once upon a time…

… I set out to the exciting world of optimization!

From SCIP 2.0 to 7.0…

… from SoPlex 1.5 to 5.0!

SCIP 2.0 to 7.0 and SoPlex 1.5 to 5.0

• SoPlex:
  • bound flipping ratio test
  • LP solution polishing
  • improved row representation
  • better sparsity exploitation
  • increased numerical stability
• SCIP:
  • new LP interface to SoPlex
  • persistent scaling
  • TreeD visualization tool
  • PySCIPOpt Python interface

• bound flipping ratio test is an advanced Simplex component that was missing
• solution polishing is a feature for SCIP but implemented entirely in SoPlex
• row representation is available by default
• SoPlex runs more stable and can make better use of sparse problems
• new LP interface enables better communication between SCIP and SoPlex
• persistent scaling improves performance and stability
• TreeD is a novel visualization tool
• PySCIPOpt provides better usability and enables rapid prototyping

Introduction to the problem

Linear programming

\[ \begin{align} \max && c^Tx & \\ \text{subject to} && Ax & \leq b \\ && x &\geq 0 \\ && \phantom{x_j} & \phantom{\in \mathbb{Z}, \forall j \in \mathcal{I}} \end{align} \]

• \(A\in\mathbb{R}^{m,n},\; b\in\mathbb{R}^m,\; c\in\mathbb{R}^n\)
• LP in canonical form
• two common ways to solve:

Simplex algorithm:

📑 Dantzig, 1947

Barrier method:

📑 Kamarkar, 1984

• simplex walks from one vertex of the polyhedron towards the next
• barrier follows an interior path until it’s close enough to the optimal face
• vast simplification: needs simplex-like crossover to generate a vertex solution

Mixed-integer linear programming

\[ \begin{align} \max && c^Tx & \\ \text{subject to} && Ax & \leq b \\ && x & \geq 0 \\ && \color{blue}{x_j} & \color{blue}{\in \mathbb{Z}, \forall j \in \mathcal{I}} \end{align} \]

• additional integrality constraints
• common way to solve: branch-and-cut

Branch-and-bound:

📑 Land and Doig, 1960

• only the two main components of solvers

Mixed-integer linear programming

\[ \begin{align} \max && c^Tx & \\ \text{subject to} && Ax & \leq b \\ && x & \geq 0 \\ && \color{RoyalBlue}{x_j} & \color{RoyalBlue}{\in \mathbb{Z}, \forall j \in \mathcal{I}} \end{align} \]

• additional integrality constraints
• common way to solve: branch-and-cut

Branch-and-bound:

📑 Land and Doig, 1960

Cutting plane separation:

📑 Gomory, 1958

• only the two main components of solvers

Visualization of solving process

• branch-and-cut tree with node LP objective value as height value
• bounds moving closer together
• no new nodes outside the primal and dual bounds
• cut progress in the root can be seen
• makes it clear that this is an LP-based branch-and-cut algorithm

SCIP 7.0 solving the MIPLIB3 instance lseu

SCIP Optimization Suite

• SCIP provides more functionality beyond MILP solving
• MILP still substantial part of the code base

Benchmark sets for computational experiments

LP (1336 instances):

• Netlib benchmark set
  http://www.netlib.org/lp/
• Csaba Mészáros’s LP collection
  http://www.sztaki.hu/~meszaros/public_ftp/lptestset/
• COR@L w/o integrality
  http://coral.ie.lehigh.edu/data-sets/mixed-integer-instances/
• MIPLIB¹, MIPLIB 3², MIPLIB 2003³, MIPLIB 2010⁴, and MIPLIB 2017⁵ w/o integrality

MILP (240 instances):

• MIPLIB 2017 benchmark set
  https://miplib.zib.de/

📑 1: Bixby, Boyd, and Indovina, 1992
📑 2: Bixby, Ceria, et al., 1998
📑 3: Achterberg, Koch, and Martin, 2006

📑 4: Koch, Achterberg, et al., 2011
📑 5: Gleixner, Hendel, et al., 2021

• multiple random seeds to reduce noise and get more reliable results

SoPlex

Spotlight on some selected features

• SoPlex is the default LP solver in SCIP and part of the SCIP Opt Suite
• implements the revised primal and dual simplex algorithm

Row representation

• define the Simplex basis not as subset of columns but of rows
• basic rows are active constraints defining the current vertex

A:

Basis:

	column rep	row rep
entering type	primal simplex	dual simplex
leaving type	dual simplex	primal simplex

• new features often need to be implemented twice

• smaller basis size has computational advantages like smaller LU
• both representations are equivalent and interchangeable

📑 Wunderling, 1996

Row representation

• define the Simplex basis not as subset of columns but of rows
• basic rows are active constraints defining the current vertex

A:

Basis:

	column rep	row rep
entering type	primal simplex	dual simplex
leaving type	dual simplex	primal simplex

• new features often need to be implemented twice

• smaller basis size has computational advantages like smaller LU
• both representations are equivalent and interchangeable

📑 Wunderling, 1996

Computational results

• LP test set
• minimum solving times of either row or column representation
• in relation to the row-column ratio (x-axis)
• trend towards more wins for col-rep for suited problem instances
• row rep more often wins very easy instances below 1 second

Computational results

• always column representation:
  8% slowdown
• always row representation:
  43% slowdown
• automatic switch to row when row/column ratio > 1.2

• minimum solving times of either row or column representation
• trend towards more wins for col-rep for suited problem instances
• row rep more often wins very easy instances below 1 second
• I worked on making the row rep more reliable and stable to allow for automatic switching (too buggy before)

Bound flipping ratio test

• also known as long step rule
• allows increasing the progress made in a single iteration in the dual simplex
• combine several iterations into one

• impact on LP benchmark:
  • 17% fewer iterations
  • 7% speedup

• idea has been around since late 70s
• ratio test determines step length; go as far as possible
• BFRT allows to go even further by flipping bounds when sign changes
• only works with boxed variables and ranged rows
• ratio test is crucial to stability of simplex
• hard to get right and easy to get wrong
• implemented for row and column version of dual simplex
• 17% iteration reduction and 7% speedup on allLP

📑 Kirillova, Gabasov, and Kostyukova, 1979

📑 Fourer, 1994

Impact of LP in solving MILPs

Persistent scaling

• one scaling factor for every row and every column: \(A'=RAC,\; R,C\) diagonal
• keep LP scaling factors for the entire MILP solving process
• pass all LP modifications through scaling layer
• better numerics and performance for cuts and node LPs
• solves 9 more instances on MIPLIB 2017

• update scaling factors for added variables and constraints
• was one reason to write a new LP interface between SCIP and SoPlex
• persistent scaling leads to largest performance improvement
• BFRT performance does not carry over to SCIP

LP Solution Polishing

• not all solutions are equal
• additional iterations to find less fractional solutions (more non-basic integer variables)
• performed internally by SoPlex
• SCIP provides integrality information

Not to be confused with [MIP] Solution Polishing by Ed Rothberg, An evolutionary algorithm for polishing mixed integer programming solutions, INFORMS Journal on Computing (2007).

• runs primal Simplex
• customized pricing and ratio test to favor non-basic integer variables

📑 Gamrath, Berthold, and Salvagnin, An exploratory computational analysis of dual degeneracy in mixed-integer programming, 2020

📑 Achterberg, pump-reduce, 2013

Computational results on MIPLIB 2017

• downward pointing bars are good
• mixed results: some outliers with strongly increased runtimes
• half of first LPs benefit from reduced fractionality
• more than 80% reduction on some instances

First LP effect on MIPLIB 2017

• most increase in time and iters when also fractionality was reduced

Tabular results

• 3 more instances solved, 6% overall speedup, 8% node reduction

Different LP solvers in SCIP (MIPLIB 2017)

1. cumulative comparison combining solved instances and timeouts
2. right section shows the respective remaining gap after one hour
3. Clp and MOSEK are performing worse than the “big 3” and even SoPlex
4. Clp suffers from many instable solves that are not shown here
5. next chart might be surprising…

Pure LP comparison

1. SoPlex stands out with the worst performance
2. can also be seen on Hans Mittelmann’s LP benchmarks but way less pronounced (factor 20 vs 2 or 3)
3. larger problems are harder for SoPlex
4. MILPs are often smaller than LPs
5. not pictured: time spent in LP is very comparable among solvers

Effect on the root gap (MIPLIB 2017)

LP solver	wins
CPLEX	79
SoPlex	59
MOSEK	55
Gurobi	53
Xpress	51
Clp	48

177 instances finished root node after one hour

1. instances are sorted by best gap over all solvers (nodelimit 1)
2. already big differences after the root node processing
3. SoPlex takes the second place in finding a good root gap
4. pure LP speed is not enough

• air05: CPLEX finds better primal bound in root
• net12: CPLEX has fewer fracts
• gmu-35-50: different primal heuristics

Numerics in
Branch & Bound & Cut

Joint work with Ted Ralphs and Dan Steffy

A difficult question…

… is it Laurel or Yanny?

• tough question with no clear answer

A numerically difficult question…

Consider this feasibility problem (no objective function):

\[ \begin{align} x + 10^{-8}y &= 10^{-7}\\ x, y &= 0 \end{align} \]

• mathematically infeasible (because \(0 + 0 \neq 10^{-7}\))

• using a feasibility tolerance of \(\epsilon=10^{-6}\):

\(x\) is basic

• \(y=0\) non-basic
• \(x=10^{-7}\) basic
• feasible solution (within \(\epsilon\))

\(y\) is basic

• \(x=0\) is non-basic
• \(y=10\) basic
• infeasible solution

-> The result depends on the chosen Simplex basis!

• certainly important to deal with numerical issues
• ideally, identifying problematic instances before running into problems

How to measure numerical difficulty?

Definition of condition number \(\kappa\) of a square matrix \(A\)

\(\kappa(A) = \|A\|\cdot\|A^{-1}\|\) provides upper bound on input error amplification in the solution of a linear system \(Ax=b\). It is the normed reciprocal of the distance to singularity and applicable to LP/MILP because of their reliance on solving linear systems with the Simplex basis matrix.

\(\kappa\):

\(\;1 \dots 10^7\)
fine

\(\;10^7 \dots 10^{10}\)
suspicious

\(\;10^{10} \dots 10^{14}\)
unstable

\(\;>10^{14}\)
ill-posed

• upper bound is typically not reached
• machine precision \(10^{-16}\) and tolerance \(10^{-6}\) to handle \(\kappa<10^{10}\)
• ideal: instance-specific, solver-independent metric for MILPs

• usual strategy: sample Kappa in the tree
• dependent on the solver and its behavior
• condition number computation is at least the cost of solving the problem
• usually approximated
• the condition of the condition is the condition number

LP approach

• compute a \(\kappa_\text{LP}\) value that captures the entire problem data \(d=(A,b,c)\): \[ \begin{align} \kappa_\text{LP}(d) &:= \frac{\left\|d\right\|}{\min\left\{\rho_P(d), \rho_D(d)\right\}}, \\[1em] \|d\| &:= \max\left\{\|A\|_{\infty,1}, \|b\|_1, \|c\|_1\right\},\\[.5em] \rho_P(d) &:= \min_{i\in\left\{1,\dots,m\right\},\;j\in\left\{-1,1\right\}} \min_{y,s,v} \max\left\{\left\|A^T y+s-\right\|_1, \left|b^T y - v\right|\right\} \\[-1em] & \hspace{8.8em}\text{subject to}\quad y_i=j, \\[.5em] \rho_D(d) &:= \min_{i\in\left\{1,\dots,n\right\},\;j\in\left\{-1,1\right\}} \min_{x,p,g} \max\left\{\left\|Ax - p \right\|_1, \left|c^T x + g\right|\right\} \\[-1em] & \hspace{8.6em}\text{subject to}\quad x_i=j \end{align} \]

• normed reciprocal to the minimum of the distances to primal or dual feasibility
• the bounds of the variables are also considered but left out better readability
• they showed a relation to the number of barrier iterations required to solve the instance

📑 Ordóñez and Freund: “Computational Experience and the Explanatory Value of Condition Measures for Linear Optimization”. SIAM Journal on Optimization (2003)

Applied to MILP (ignoring integrality)

• both CPLEX and Xpress implement MIP-\(\kappa\)
• hard to make any reliable statements
• too many problems are “ill-conditioned” \((\kappa_\text{LP}\geq 10^{20})\)
• attention level learning appears to be a promising alternative
  📑 Berthold, Numerics IV: Learning to pay attention, 2020

1. CPLEX and Xpress implement a MIP kappa to sample the condition numbers of optimal basis matrices
2. attention level is a weighted sum of percentages of kappa values
3. matrix plot does not reveal any relationship for kappa_LP
4. attention level and max log kappa shows some information
5. use exact or high-precision solvers to handle numerically challenging models

TreeD as a side product

• started as a tool to study numerical features of MILP trees
• LP-focused visualization keeps true to the theme of the thesis
• uses custom EventHandler in PySCIPOpt to collect LP data
• Python package scikit-learn transforms high-dimensional solutions to 3D or 2D (Multi dim. scaling)
• can provide additional insight into how SCIP solves a problem (LP iters, condition numbers)
• reveals impact of cutting planes
• artistic value: quite pretty to look at
• ported to Amira for 3D-Studio DaVinci at ZIB

https://github.com/mattmilten/TreeD

Mittelmann benchmark plots

• interactive solver-by-solver comparisons¹
• largest publicly available benchmark dataset for linear, mixed-integer, nonlinear, and combinatorial problems by Hans Mittelmann²
• published results are updated on a daily basis
• easier overview of latest ranking
• additional virtual best or portfolio solver to reveal performance potential

-> MIPLIB 2017 benchmark using SCIP/spx as base

1. https://github.com/mattmilten/mittelmann-plots

2. http://plato.asu.edu./bench.html

• adds some more visibility to those benchmarks and the optimization community
• more accessible than raw tables

Conclusion

• Simplex-oriented perspective on MILP solving
• performance and stability improvements for LP in MILP
  • persistent scaling
  • solution polishing
• impact of different LP solvers in SCIP
  • pure LP performance is not decisive factor
• LP-based numerical study on MILP
  • goal: find holistic measure to determine when to branch
  • resulting in TreeD visualization
• PySCIPOpt extending the accessibility and usability of SCIP

Thank you!

1. “LP does not matter for MILP”
2. expect continued competition in the optimization landscape, fueled by Hans’ benchmarks
3. numerical study is one piece in this vast area of research that deserves more attention and understanding
4. TreeD demonstrates how PySCIPOpt can be used

Algorithms

Simplex

Solution polishing

Figures

1. Introduction

Figure 1.5: Size distribution of LP test set, including solving time colorization for SoPlex with default settings. Sparsity is computed as nonzeros/\((n*m)\)

Figure 1.6: Number of iterations and time to optimality (using a time limit of one hour) with respect to the number of nonzero elements in the problem matrices

Figures

3. Implementational Aspects of the Simplex Algorithm

Figure 3.1: Minimum solving times for row and column representation depending on row to column ratio on the allLP test set.

Figure 3.4: Pricing statistics

Figure 3.6: Distribution of variability scores > \(10^{−6}\) for solving time and iterations of SoPlex 4.0.2 across three different seeds; allLP test set

Figures

4. Impact of Linear Programming in MILP

Figure 4.2: Fraction of time spent solving LPs during MILP solving; MIPLIB 2017 benchmark

Figure 4.3: Comparison of root LP performance with SCIP/MOSEK; MIPLIB 2017 benchmark

Figure 4.4: Histogram of fraction of time spent solving LPs during MILP solving; MIPLIB 2017 benchmark

Figure 4.5: Distribution of different LP types to be solved during MILP optimization; MIPLIB 2017 benchmark

Figure 4.9: Illustrative example of MDS transformation and associated Shepard error plot

Figure 4.10: TreeD visualization of MIPLIB3 instance lseu (89 binary variables, 28 constraints)

Figure 4.11: Shepard plot of original and transformed pairwise distances of all encountered node LP solutions after applying the MDS transformation

Figure 4.12: Overall performance comparison of different LP solvers in SCIP 6.0.2 with default settings on MIPLIB 2017 benchmark.

Figure 4.13: Distribution of solving time and iteration ratios for different LP solvers in SCIP to solve the root LPs on MIPLIB 2017 benchmark.

Figure 4.14: Development of node LP solving times and simplex iterations in the tree, relative to the root LP. MIPLIB 2017 instances solved within one hour by all solvers (54 instances)

Figure 4.15: Comparison of the gap after root node processing. Node limit of one, time limit of one hour on MIPLIB 2017 benchmark (117 instances).

Figure 4.16: Comparison of the node throughput, time limit of one hour MIPLIB 2017 benchmark. No cuts, no primal heuristics. Results are restricted to instances with a tree size of at least 100.

Figures

5. LP Solution Polishing

Figure 5.2: Effect of LP solution polishing on the condition numbers of the optimal basis matrix of first LP relaxations of the MIPLIB 2017 benchmark instances. Condition number changes are differences in log10 with and w/o polishing.

Figure 5.3: Effect of solution polishing on the fractionality of the first LP solution for the MIPLIB 2017 benchmark instances. 3 seeds with SCIP 7

Figure 5.4: Effect of solution polishing on the total solving time for the MIPLIB 2017 benchmark instances. Factors below 1 indicate a speedup. 3 seeds with SCIP 7.

Figures

6. Numerics in Branch & Bound & Cut

Figure 6.3: Condition number development (vertical axis, in log scale) for every simplex iteration in the root node (horizontal axis) including re-optimizations after adding cutting planes in multiple rounds (vertical lines). A plot of objective values at each iteration is overlaid as a dashed gray line with the scales given to the right of each plot.

Figure 6.4: Root node comparison of condition numbers of the original LP and after including cutting planes. Most instances show a higher condition number during the cutting stage and when this process is completed

Figure 6.5: Aggregated condition numbers (in log10 scale) for different tree depths. MMM test set; solved in 1 hour and min tree depth of 10. Mean value of condition numbers and the respective 95% confidence interval. Breadth-first search with different cutting plane separation strategies.

Figure 6.6: Expected developments of condition numbers (in log10 scale): aggressive separation (sepaaggr) leads to an increased condition number across the entire tree while disabled cutting plane separation generates a larger tree but has smaller condition numbers.

Figure 6.7: Unexpected developments of condition numbers (in log10 scale): Disabling cutting plane generation (sepaoff) leads to worse condition numbers. The left instance exhibits decreasing condition numbers going down the tree while for the right instance, the numbers are growing.

Figure 6.8: Statistics of the κ\(_\text{LP}\) values and the tree condition numbers for various MILPs. Note that the κ\(_\text{LP}\) values are in each case the lower bounds from the estimation. The lower right subplot shows the κ\(_\text{LP}\) distribution which is independent of the solver.

PySCIPOpt

The Python interface to SCIP

Python interface to SCIP

• easy entryway to using and extending SCIP
• publicly developed on GitHub¹ under the MIT license since 2016
• conda install pyscipopt also installs the SCIP Opt Suite

• requires no manual compilation
• no C knowledge needed to work with SCIP

1. https://github.com/scipopt/PySCIPOpt

Modeling example

Mathematical formulation \[\begin{align*} &\text{minimize} & x + 3y & \\ &\text{subject to} & 2x - y^2 & \geq 10 \\ & & x, y & \geq 0 \\ & & x & \in \mathbb{R} \\ & & y & \in \mathbb{Z} \\ \end{align*}\]

Python code

from pyscipopt import Model

m = Model()
x = m.addVar("x")
y = m.addVar("y", vtype="I")
m.setObjective(x + 3*y)
m.addCons(2*x - y*y >= 10)

• allows natural and human-readable model formulations
• inspired by Gurobi’s Python interface

• rapid prototyping without recompilation
• easy access to vast library of Python modules
• design principles:
  • wrapper around SCIP’s C functions
  • additional convenience functions