Condition Number

Definition

From \[B(x+\delta x) = b + \delta b \rightarrow \delta x= B^{-1}\delta b\] Cauchy-Schwartz inequalities: \[\|\delta x\| \leq \| B^{-1}\|\|\delta b\|, \|b\| \leq \| B\|\|x\|\] Combine: \[\frac{\|\delta x\|}{\|x\|} \leq \underbrace{\| B\|\| B^{-1}\|}_{\kappa}\frac{\|\delta b\|}{\|b\|}\]

Condition Number

Definition

Condition Number

Value Range

Condition Number

Interpretations

Condition Number

Geometrical interpretation

Example 1

Large matrix coefficient ranges

\[ \begin{align} & B=\begin{pmatrix} M & 1\\ 0 & \frac{1}{M} \end{pmatrix} \text{ has inverse } B^{-1}=\begin{pmatrix} \frac{1}{M} & -1\\ 0 & M \end{pmatrix}\\ & \Rightarrow \|B\| \geq M \text{ and } \|B^{-1}\| \geq M\\ & \Rightarrow \kappa(B) = \|B\|\cdot\|B^{-1}\| \color{red}{\geq M^2} \end{align} \]

Example 1.2

Large matrix coefficient ranges

\[ \begin{align} & B=\begin{pmatrix} M & 1\\ 0 & M \end{pmatrix} \text{ has inverse } B^{-1}=? \end{align} \]

Example 2

Inconsistent or unnecessarily truncated data

• Values that should be equal have slight differences
• Example:

    c306: 0.416666667 x80 – 100 x90 = 0
   [...]
  c11360: 0.416666666 x80 – 100 x90 = 0

• Values appear to be calculated inconsistently
• Values are only calculated to 9 digits → bigger perturbation than necessary

Example 3

Hidden mixtures of large and small coefficients

• \(\kappa(B)=\|B\|\cdot\|B^{-1}\|\) is large if \(\|B\|\) or \(\|B^{-1}\|\) is large
• \(\|B^{-1}\|\) can also be large for other reasons
• Cascade or staircase: \[ B=\begin{pmatrix} 1 & -2 & & & \\ & 1 & -2 & & \\ & & \ddots & \ddots & \\ & & & 1 & -2 \\ & & & & 1 \\ \end{pmatrix} ,\; B^{-1}=\begin{pmatrix} 1 & 2 & 4 &\dots & 2^{n-1} \\ & 1 & 2 & \dots & 2^{n-2} \\ & & \ddots & \ddots & \vdots \\ & & & 1 & 2 \\ & & & & 1 \\ \end{pmatrix} \]

Quick summary

Why is this relevant?

LPs

A solution to an LP requires:

1. A basis
2. Feasibility verification or proof of infeasibility
3. Optimality verification or proof of unboundedness

Requires solving two linear systems of equations (with basis matrix \(B\)):

• Primal solution: \(\ \;x_B=B^{-1}(b-A_N x_N)\)
• Dual solution: \(\quad\pi^T = c_B^T B^{-1}\)

Why is this relevant?

LPs

For verification, \(B\) is used.

Why is this relevant?

LPs

Why is this relevant?

MIPs

graph TD;

    root(( ))-->l;
    root-->r;
    l(( ))-->ll(( )):::green;
    classDef green fill:green;
    l-->lr(( ));
    r(( ))-->rl(( )):::orange;
    r-->rr(( )):::orange;
    classDef orange fill:orange;
    lr-->lrl(( ));
    lr-->lrr(( )):::orange;
    lrl-->lrll(( )):::orange;
    lrl-->lrlr(( )):::orange;

• The LP relaxation is solved after every branching decision and cut.
• Pruning, optimality checks (these are hard!).

Floating Point Numbers

1 + pow(2,-52)

1.0000000000000002

1 + pow(2,-53)

1.0

pow(2,-53) + pow(2,-53) + 1 == 1 + pow(2,-53) + pow(2,-53)

False

(1 +  pow(2,-53)) + pow(2,-53) == 1 + (pow(2,-53)  + pow(2,-53))

False

Floating Point Numbers

s = [1, 2, 3, 4, 5]
w = [11.11, 22.22, 33.33, 44.44, 55.55]

# Multiply component wise
r1 = 0.0
for i in range(5):
    r1 += s[i] * w[i]

# Multiply component wise (in reverse order)
r2 = 0.0
for i in range(4, -1, -1):
    r2 += s[i] * w[i]

print(f"{r1=}, {r2=}")
print(f"{r1==r2=}")

r1=611.05, r2=611.0500000000001
r1==r2=False

Representation of Numbers

Representation of Numbers

float representation of 1234.5 (calculator)

Analyzing

1. Solution and Formulation
2. Optimization
3. Different runs

Analyzing

Solution & Formulation

• Model statistics: model.printStats()
• Solution quality: model.printQuality()
  • Many more Quality Attributes
• Check the condition number.
  • We will see how later…

Analyzing

Optimization

Prior to optimization

Optimize a model with 306204 rows, 579547 columns and 1563587 nonzeros
Model fingerprint: 0x21cdd0c7
Coefficient statistics:
  Matrix range     [1e-05, 2e+08]
  Objective range  [8e-03, 1e+05]
  Bounds range     [1e-03, 3e+09]
  RHS range        [5e-01, 2e+04]
Warning: Model contains large matrix coefficient range
Warning: Model contains large bounds
         Consider reformulating model or setting NumericFocus
         parameter to avoid numerical issues.

Analyzing

Optimization

During LP solving (simplex or barrier)

Iteration    Objective       Primal Inf.    Dual Inf.      Time
[...]
   56658   -3.3179831e+03   4.239720e+05   0.000000e+00     20s
Warning: Markowitz tolerance tightened to 0.0625
   67557   -3.2801023e+03   1.974011e+04   0.000000e+00     25s
   69517    3.7944349e+02   6.207667e+04   0.000000e+00     30s
Warning: 1 variables dropped from basis
Warning: switch to quad precision
   78320    1.7027054e+08   1.149199e+04   0.000000e+00     35s

Barrier performed 49 iterations in 0.03 seconds (0.02 work units)
Numerical trouble encountered

Analyzing

Optimization

During the MIP node process

   Nodes    |    Current Node    |     Objective Bounds      |     Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
[...]
  3252     0 278173.332    0 2120 275179.027 277183.668  0.73%   211  401s
  3253     0 278173.020    0 2060 275179.027 277115.904  0.70%   211  405s
Relaxation Markowitz Tolerance tightened to 0.5
  3254     0 278173.020    0 2060 275179.027 277115.904  0.70%   211  406s
[...]
  3575   295  postponed   43      275179.027 277115.904  0.70%   321  461s
  3599   314 infeasible   44      275179.027 277115.904  0.70%   332  467s
  3673   391 276640.457   55 1881 275179.027 277115.904  0.70%   328  471s

Analyzing

Optimization

Post Optimization

Solved in 131511 iterations and 75.99 seconds (66.11 work units)
Optimal objective  1.730551793e+08
Warning: unscaled dual violation = 0.0141201 and residual = 0.0111682

Analyzing different runs

• What happens if we run the same model using different seeds?
  • Does the solution remain the same?
  • Does the time remain similar?
    • This may be due to other factors

Analyzing different runs

Two different solutions for harp2 (using MIPGap=0).

2c2
< # Objective value = -73899798
---
> # Objective value = -73899799.3523847
41,42c41,42
< x39 0
< x40 1
---
> x39 1
> x40 0
451,452c451,452
< x449 1
< x450 0
---
> x449 0
> x450 1
697,698c697,698
< x695 1
< x696 0
---
> x695 0
> x696 1

Model Reformulation

The problem with big-M

import gurobipy as gp

M = 1e10
opts = {"OutputFlag": 0}

with gp.Env(params=opts) as env, gp.Model(env=env) as m:
    c = m.addVar(name="c", lb=0)
    b = m.addVar(name="b", vtype=gp.GRB.BINARY)
    # Fix variables
    c.UB = 1e4
    c.LB = 1e4
    b.UB = 1e-6 # This is technically zero based on the IntFeasTol
    b.LB = 1e-6
    m.addConstr(c <= M*b) # So we really want c = 0 since b = 0
    m.optimize()
    print(f"Is the model optimal? {m.Status==gp.GRB.OPTIMAL}")

Is the model optimal? True

Model Reformulation

Alternative 1

import gurobipy as gp

M = 1e10
opts = {"OutputFlag": 0}

with gp.Env(params=opts) as env, gp.Model(env=env) as m:
    c = m.addVar(name="c", lb=0)
    b = m.addVar(name="b", vtype=gp.GRB.BINARY)
    # Fix variables
    c.UB = 1e4
    c.LB = 1e4
    b.UB = 1e-6 # This is technically zero based on the IntFeasTol
    b.LB = 1e-6
    m.addConstr((b == 0) >> (c == 0))
    m.optimize()
    print(f"Is the model optimal? {m.Status==gp.GRB.OPTIMAL}")

Is the model optimal? False

Model Reformulation

Alternative 2

import gurobipy as gp

M = 1e10
opts = {"OutputFlag": 0, "IntFeasTol": 1e-7}

with gp.Env(params=opts) as env, gp.Model(env=env) as m:
    c = m.addVar(name="c", lb=0)
    b = m.addVar(name="b", vtype=gp.GRB.BINARY)
    # Fix variables
    c.UB = 1e4
    c.LB = 1e4
    b.UB = 1e-6 # This is technically zero based on the IntFeasTol
    b.LB = 1e-6
    m.addConstr(c <= M*b) # So we really want c = 0 since b = 0
    m.optimize()
    print(f"Is the model optimal? {m.Status==gp.GRB.OPTIMAL}")

Is the model optimal? False

There are many parameters

… and choosing the right ones is no easy task

Parameter recommendations

NumericFocus = [0, 1, 2, 3]

Meta parameter that affects many aspects of the solution process:

Presolve

• Tighter tolerances when fixing variables
• Avoid small coefficients in bound strengthening
• More restrictive checks in aggregator

Parameter recommendations

Presolve = [-1, 0, 1, 2]

The main culprits are:

• Aggregate - may introduce accumulation of numerical errors
  • Try setting this to 0 to deactivate aggregation

Parameter recommendations

Root relaxation or pure LP: Simplex

• MarkowitzTol: pivoting tolerance (for LU factorization)
  • Larger values increase numerical stability at the expense of sparsity

Remember the warning:

Iteration    Objective       Primal Inf.    Dual Inf.      Time
[....]
   56658   -3.3179831e+03   4.239720e+05   0.000000e+00     20s
Warning: Markowitz tolerance tightened to 0.0625
[....]

• Gurobi dynamically adjusts this to balance performance and stability

Parameter recommendations

Root relaxation or pure LP: Simplex

• Quad = [-1, 0, 1]
  • Higher precision at the expense of speed and memory
  • -1 switches dynamically to quadruple precision if numerical issues arise

   69517    3.7944349e+02   6.207667e+04   0.000000e+00     30s
Warning: 1 variables dropped from basis
Warning: switch to quad precision
   78320    1.7027054e+08   1.149199e+04   0.000000e+00     35s
   88724    1.7034414e+08   1.389671e+04   0.000000e+00     40s

Parameter recommendations

Root relaxation or pure LP: Barrier

• BarHomogeneous = [-1, 0, 1]
  • Default or homogeneous algorithm
  • Homogeneous algorithm is usually a bit slower
  • Homogeneous is better in recognizing infeasibility or unboundedness, or dealing with models that are at the boundary of infeasibility

Parameter recommendations

MIP

• Cuts = [-1, ..., 3]
  • 0: Disable all cuts (together with CutPasses = 0 to be completely shut off, may affect heuristics)
  • 1: Moderate cut generation
  • 2, 3: Aggressive and even more aggressive

Parameter recommendations

MIP

• IntegralityFocus = [0, 1]
  • Enables extra computations to avoid trickle flow (like in the big-M example)
  • Generate a “cleaner” solution with fewer tiny solution values
  • Typically safe to enable but may be expensive

Parameter recommendations

MINLP

• FuncNonlinear = [0, 1]
  • Use outer approximation instead of piecewise-linear approximation
  • More accuracy, potentially more expensive
  • When enabled, parameters FuncPieces, FuncPieceError, FuncPieceLength, FuncPieceRatio are obsolete
  • only available in Gurobi 11 and when using nonlinear function contraints

Practical recommendations:

Summary

Analyze:

Practical recommendations:

Summary

• Try stabilizing the solving behavior via different parameters
• Run multiple seeds and slight variations of the model

Short intermission: FeasRelax

Basic idea

Approximiate the distance to singularity

Apply FeasRelax to minimize the sum of infeasibilities of:

Challenges

• The FeasRelax problem is an LP whose structural columns comprise an ill-conditioned basis matrix
  • Set NumericFocus=3 and ScaleFlag=2 if row and column ratios are large
  • Provide more custom parameters via PRM file

• Too many rows or columns in the output; hard to determine source of ill-conditioning
  • Try different explanation method (row, column, angle)
  • Sometimes one method provides a much smaller explanation

Quick Start Guide

• Installation: python -m pip install gurobi-modelanalyzer

• Basic usage:

  import gurobipy as gp
  import gurobi_modelanalyzer as gma
  m = gp.read("badmodel.mps")
  m.optimize()
  print(m.kappa)                         # print kappa of final basis
  gma.kappa_explain(m)                   # row-based explanation
  gma.kappa_explain(m, expltype="COLS")  # column-based explanation
  print(gma.angle_explain(m))            # angle explanation

Output on maliciously modified afiro

• kappa_explain() will generate a new LP or MPS file, containing the ill-conditioning certificate:
• afiro_kappaexplain.lp:

Minimize
  0 X36 + 0 X04 + 0 X15 + 0 X16 + 0 X26 + 0 X38 + 0 X37
Subject To
 GRB_Combined_Row: 0.0303868836044176 X23 + 4.80518e-10 X01
   - 4.65661e-10 X03 = 0
 (mult=2696322.968477607)R09x: - 0.9999999000000001 X01 + X03 = 0
 (mult=-2696322.6896988587)R09: - X01 + X03 = 0
 (mult=0.2787787486643817)X46: - X03 + 0.109 X22 <= 0
 (mult=0.030386883604417606)R19: X23 - X22 + X24 + X25 = 0
 (mult=0.030386883604417606)X45: - X25 <= 0
 (mult=0.030386883604417606)X48: 0.301 X01 - X24 <= 0
Bounds
End

Angle explanation on modified afiro instance

• angle_explain() returns a list of tuples:
  1. almost parallel rows
  2. almost parallel columns
  3. associated model

  (
    [(<gurobi.Constr R09>, <gurobi.Constr R09x>)],
    [],
    <gurobi.Model Continuous instance basismodel:
      28 constrs, 28 vars, No parameter changes>
  )

“Real” example: neos-1603965 (MIPLIB 2017)

• Use model.relax() to relax integrality conditions
• Final condition number 1.60000000016e+21 (original model)
• angle_explain() does not reveal any almost parallel rows or columns

Subject To
 GRB_Combined_Row: = 1.000000015655
 (mult=0.999999999815)R24991: C8008 - C8009 >= 0
 (mult=9.9999999981499e-11)R9004: - 1e+10 C8008 <= 0
 (mult=9.9999999981499e-11)R13004: - 0.15 C7002 + 1e+10 C8009 <= 1e+10
 (mult=1.4999999997225e-11)R1030: C3030 - C7001 <= 0
 (mult=-1.4999999997225e-11)R1966: C3966 - C7001 <= 0
 (mult=-1.4999999997225e-11)R2966: C4966 - C7002 <= 0
 (mult=-1.4999999997225e-11)R28014: C0030 + C3030 = 7165
 (mult=1.4999999997225e-11)R28950: C0966 + C3966 + C4966 = 8221
 (mult=1.12499999982e-11)R0030: 1.3333333333 C0030 = 0
 (mult=-1.12499999982e-11)R0966: 1.3333333333 C0966 = 0

• Don’t ignore rows with small \(y\): big-M values of \(10^{10}\) are the issue here

More functionality

• converttofractions():
  • takes a list of decimal values and tries to represent them as fractions
  • useful to clean-up “dirty” coefficients
• matrix_bitmap():
  • wrapper for spy() (matplotlib) to inspect sparsity pattern
  • useful to identify cascades

Summary

• Track down reasons for numerical issues

• Open-source (Apache-2.0) on GitHub:

  https://github.com/Gurobi/gurobi-modelanalyzer

• Install via pip:

  python -m pip install gurobi-modelanalyzer

• Currently only supports LPs (use model.relax() to inspect root LP of MIPs)

Further Material

• Klotz, E. (2014). Identification, assessment, and correction of ill-conditioning and numerical instability in linear and integer programs. In Bridging Data and Decisions (pp. 54-108). INFORMS.
• Gurobi Tech Talk: Converting Weak to Strong MIP Formulations, YouTube, Materials.
• Gurobi Holiday Tech Talk: Prevent the Grinch from Stealing Your Optimization Projects, YouTube, Materials.