MIP: presolve fixes a column whose bounds are equal only within tolerance without bounding the objective error, returning a suboptimal solution as Optimal

[EamonHetherton]

Summary

Note: Authored in conjunction with Claude AI

HPresolve::checkColBounds() treats a column as fixed when its bounds are equal only within the
feasibility tolerance. The guard bounds the resulting perturbation of the row activities, but
not of the objective. A column can be nearly fixed in the constraints while still being worth a
great deal in the objective, and colPresolve() then pins it at its lower bound, losing up to
|cost| * (ub - lb).

On the attached 40x33 model that loses 2496.5 — presolve reduces the model to empty, reports
Presolve: Optimal, and returns an objective 6.3% away from the true optimum, with default
options and 0 branch-and-bound nodes.

Provenance, in the interests of not overselling this. I did not hit this in my own models. It
turned up in a randomly generated model, from a generator I wrote to produce minimal reproducers
for #3170 and #3171. On the two real models that motivated those two issues, this condition is
reached zero times, so this is not a defect I have observed in practice.

That said, I do not think it is only an artefact of odd generated data. The column being fixed has
a cost of -7.29e10, but the model's own cost range is [1e-06, 1e+06] — the coefficient is about
73,000x larger than anything in the input, built up by presolve's own substitutions before the
column reaches checkColBounds(). So reaching this state does not require an extreme cost in the
MPS file, only an awkward one plus a few substitutions.

Reproducing

3173-1.mps.txt is attached — rename it to 3173-1.mps. 40 rows, 33 columns, 138 nonzeros,
8 binaries.

$ highs 3173-1.mps

Presolve reductions: rows 0(-40); columns 0(-33); nonzeros 0(-138) - Reduced to empty
Presolve: Optimal

  Status            Optimal
  Primal bound      -37333.0762329
  Dual bound        -37333.0762329
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Nodes             0

The true optimum is -39829.5816585168. presolve=off returns it.

The optimum was verified independently of the MIP search, by enumerating all 2^8 binary assignments
and solving the resulting LP for each — so the claim does not rest on trusting branch-and-bound:

import itertools, highspy

h = highspy.Highs(); h.setOptionValue("output_flag", False)
h.readModel("3173-1.mps")
lp = h.getLp()
ints = [j for j in range(lp.num_col_)
        if lp.integrality_[j] == highspy.HighsVarType.kInteger]

best = float("inf")
for bits in itertools.product((0.0, 1.0), repeat=len(ints)):
    g = highspy.Highs(); g.setOptionValue("output_flag", False)
    g.readModel("3173-1.mps")
    for j, v in zip(ints, bits):
        g.changeColBounds(j, v, v)
        g.changeColIntegrality(j, highspy.HighsVarType.kContinuous)
    g.run()
    if g.getModelStatus() == highspy.HighsModelStatus.kOptimal:
        best = min(best, g.getObjectiveValue())

h.run()
print("brute-force optimum:", best)                  # -39829.5816585168
print("HiGHS reports      :", h.getObjectiveValue())  # -37333.0762503876

Also reproduces on 1.13.1 (-37333.076250).

Cause

highs/presolve/HPresolve.cpp, checkColBounds():

double boundDiff = model->col_upper_[col] - model->col_lower_[col];
if (boundDiff <= primal_feastol &&
    (boundDiff <= options->small_matrix_value ||
     getMaxAbsColVal(col) * boundDiff <= primal_feastol)) {
  ...
  if (isFixed != nullptr) *isFixed = true;
}

and colPresolve(), which acts on it:

  HPRESOLVE_CHECKED_CALL(checkColBounds(col, &isFixed));
  if (isFixed) {
    postsolve_stack.removedFixedCol(col, model->col_lower_[col], ...);
    removeFixedCol(col);
    return checkLimits(postsolve_stack);
  }

  switch (colsize[col]) {
    case 0:
      return emptyCol(postsolve_stack, col);   // never reached for this column

The state on the attached model, at the point the column is removed:

lb           1.3502174835205079
ub           1.3502175177853328
boundDiff    3.426482e-08        <= mip_feasibility_tolerance (1e-7)  -> "fixed"
maxAbsColVal 0                   <- the column is empty by this stage
cost         -72859132125.478302

|cost| * boundDiff = 2496.5050, against an observed objective error of 2496.5054.

Two things compound here:

1. the column is empty by this stage, so getMaxAbsColVal(col) * boundDiff is 0 <= 1e-7,
   the row-activity guard is vacuously satisfied no matter how large boundDiff is;
2. checkColBounds() runs before the switch (colsize[col]) that would dispatch to
   emptyCol(), which exists precisely to fix an empty column at its cost-optimal bound. The
   tolerance check pre-empts the routine that would have got it right, and pins at the lower bound.

Dropping mip_feasibility_tolerance below boundDiff, so the column is no longer treated as
fixed, removes the error entirely:

mip_feasibility_tolerance	objective
1e-6	-37333.0762329 (wrong)
1e-7	-37333.0762329 (wrong)
1e-8	-39829.5816345 (correct)
1e-9	-39829.5816345 (correct)

Version

master @ 04024d701f and latest @ bc30ccf1a1 (1.15.1); also reproduces on 1.13.1.

Related

#3170 and #3171 are two invalid-cut defects in MIP cut generation, found on the same class of
model. This one is separate: it is decided entirely in presolve, before branch-and-bound runs, and
the fixes for those two leave it unchanged.

I have a branch carrying this reproducer as a failing regression test, and a candidate change, but
the appropriate threshold for the missing condition is a judgement call I would rather leave to
you, so I have kept it out of this report. It is written up on the PR instead.