From 68d99773f00b805b84602be7439cfb85799cb95d Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Thu, 6 Aug 2026 01:46:47 -0700 Subject: [PATCH 1/8] BHW coefficient reduction pass --- cpp/src/mip_heuristics/CMakeLists.txt | 1 + .../presolve/bhw_coeff_reduce.cpp | 582 ++++++++++++++++++ .../presolve/bhw_coeff_reduce.hpp | 101 +++ .../presolve/third_party_presolve.cpp | 2 + cpp/tests/internal/CMakeLists.txt | 1 + cpp/tests/mip/bhw_coeff_reduce_test.cpp | 171 +++++ 6 files changed, 858 insertions(+) create mode 100644 cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp create mode 100644 cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp create mode 100644 cpp/tests/mip/bhw_coeff_reduce_test.cpp diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index 6ad1009d84..f54619cfe3 100644 --- a/cpp/src/mip_heuristics/CMakeLists.txt +++ b/cpp/src/mip_heuristics/CMakeLists.txt @@ -13,6 +13,7 @@ set(MIP_LP_NECESSARY_FILES ${CMAKE_CURRENT_SOURCE_DIR}/local_search/rounding/simple_rounding.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/third_party_presolve.cpp ${CMAKE_CURRENT_SOURCE_DIR}/presolve/gf2_presolve.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bhw_coeff_reduce.cpp ${CMAKE_CURRENT_SOURCE_DIR}/solution/solution.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/conflict_graph/clique_table.cu ) diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp new file mode 100644 index 0000000000..862a5b61cb --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp @@ -0,0 +1,582 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "bhw_coeff_reduce.hpp" + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Gordon H. Bradley, Peter L. Hammer, Laurence Wolsey (1974) Coefficient reduction for inequalities +// in 0-1 variables. Mathematical Programming 7:263-282. +// +// Over binaries many inequalities share the same 0/1 feasible set. Theorem 2.5 characterizes that +// set of "equivalent" inequalities by a system over the ceilings and roofs of the row. The test +// here is the same characterization taken over the plain maximal feasible and minimal infeasible +// points, which are supersets of BHW's ceilings and roofs because we skip their ordering condition. +// That costs nothing: with all coefficients positive the activity is monotone, so the maximum of +// w.x over the feasible points is attained at a maximal one and the minimum over the infeasible +// points at a minimal one. Checked exhaustively over 9.4M (row, weight) pairs against full 2^k +// equivalence, in both directions, before landing. +// +// BHW's minimization is not implemented. Section 5 obtains the minimum equivalent inequality by LP +// over the polytope. We instead search weight vectors by increasing max|w| up to a cap, which is +// minimal by construction below the cap, and fall back to two heuristic candidates above it. That +// search is affordable at every width we enumerate only because conditions N1 and N3 restrict +// candidates to non-negative non-increasing sequences: at most sum_{m=1..6} C(11+m, m) = 18563 +// vectors for a 12-entry row, measured at 187us for the worst shape. Lemma 3.6 supplies the lower +// bound that seeds and prunes the search. +// +// TODO: extend the fallback path to the row generation of Section 6. On the PaPILO-reduced corpus +// the LP reaches 1489 reducible rows where these heuristics match the optimum on roughly 60%, and +// unlocks int8 on 321 of them. The separation oracle in step 3 is a 0-1 knapsack, max w.x subject +// to a.x <= b, which solve_knapsack_problem in cuts/cuts.hpp already solves by DP. Note that +// routine is the wrong tool for the LP-strength gate below, which needs the continuous relaxation. +// Boyd (1993) Generating Fenchel Cutting Planes for Knapsack Polyhedra, SIAM J. Optim. +// 3(4):734-750, treats the same separation problem as cut generation. + +namespace cuopt::mathematical_optimization::mip { + +namespace { + +// One row in BHW's normalized frame: condition N1 (every coefficient positive, negatives +// complemented by x_i = 1 - y_i) and condition N3 (coefficients sorted descending). N1 is what +// makes the row activity monotone in x, and it also makes the weights of any equivalent inequality +// non-negative, since a variable that tightens the row must tighten every inequality with the same +// feasible set. Both normalizations are undone before any reduction is returned. +struct norm_row_t { + int k = 0; + std::array coef{}; // descending, all > 0 + std::array slot{}; // position of this entry in the caller's row arrays + std::array flipped{}; // whether this entry was complemented by N1 + int64_t rhs = 0; +}; + +struct partition_t { + std::vector maximal_feasible; + std::vector minimal_infeasible; + // Lemma 3.6: where coef[i] > coef[i+1] and the two variables are not symmetric, every equivalent + // inequality has w_i >= w_{i+1} + 1. Chaining those steps down to w_{k-1} >= 0 bounds max|w|. + std::array strict{}; + std::array suffix_strict{}; + int lemma36_bound = 0; +}; + +int64_t weight_activity(const int64_t* w, uint32_t mask) +{ + int64_t sum = 0; + while (mask != 0u) { + sum += w[std::countr_zero(mask)]; + mask &= mask - 1u; + } + return sum; +} + +// Splits {0,1}^k and collects the maximal feasible and minimal infeasible points. Returns false for +// a degenerate row (all points feasible or all infeasible), which is left to other presolvers. +bool build_partition(const norm_row_t& row, partition_t& out) +{ + const int k = row.k; + const uint32_t n_pat = 1u << k; + cuopt_assert(k >= 2 && k <= BHW_MAX_LEN, "row length outside the enumerable range"); + + std::vector activity(n_pat, 0); + std::vector feasible(n_pat, 0); + uint32_t n_feasible = 0; + for (uint32_t m = 1; m < n_pat; ++m) + activity[m] = activity[m & (m - 1u)] + row.coef[std::countr_zero(m)]; + for (uint32_t m = 0; m < n_pat; ++m) { + feasible[m] = activity[m] <= row.rhs ? 1 : 0; + n_feasible += feasible[m]; + } + if (n_feasible == 0 || n_feasible == n_pat) return false; + + for (uint32_t m = 0; m < n_pat; ++m) { + bool extremal = true; + if (feasible[m] != 0) { + for (int i = 0; i < k && extremal; ++i) + if ((m >> i & 1u) == 0u && feasible[m | (1u << i)] != 0) extremal = false; + if (extremal) out.maximal_feasible.push_back(m); + } else { + for (int i = 0; i < k && extremal; ++i) + if ((m >> i & 1u) != 0u && feasible[m ^ (1u << i)] == 0) extremal = false; + if (extremal) out.minimal_infeasible.push_back(m); + } + } + cuopt_assert(!out.maximal_feasible.empty() && !out.minimal_infeasible.empty(), + "a non-degenerate partition has at least one extremal point on each side"); + + for (int i = 0; i + 1 < k; ++i) { + if (row.coef[i] <= row.coef[i + 1]) continue; + const uint32_t lo_bit = 1u << i; + const uint32_t hi_bit = 1u << (i + 1); + for (uint32_t m = 0; m < n_pat; ++m) { + if ((m & lo_bit) != 0u || (m & hi_bit) == 0u) continue; + // coef[i] > coef[i+1], so moving the set bit down raises the activity: a feasible point whose + // swap is infeasible witnesses that the two variables are not interchangeable. + if (feasible[m] != 0 && feasible[(m ^ hi_bit) | lo_bit] == 0) { + out.strict[i] = true; + break; + } + } + } + for (int i = k - 1; i >= 0; --i) + out.suffix_strict[i] = out.suffix_strict[i + 1] + (out.strict[i] ? 1 : 0); + out.lemma36_bound = out.suffix_strict[0]; + cuopt_assert(out.lemma36_bound < k, "at most k-1 strict steps exist in a row of length k"); + return true; +} + +// BHW Theorem 2.5: (w, t) is equivalent to the row iff every maximal feasible point M satisfies +// sum_{i in M} w_i <= t and every minimal infeasible point satisfies sum_{i in .} w_i >= t + 1. +// Taking t as the feasible-side maximum makes the first block hold by construction, leaving hi < +// lo; integer weights make "> t" and ">= t+1" the same statement. +bool accepts(const partition_t& part, const int64_t* w, int64_t& bound) +{ + int64_t hi = std::numeric_limits::min(); + for (uint32_t m : part.maximal_feasible) + hi = std::max(hi, weight_activity(w, m)); + int64_t lo = std::numeric_limits::max(); + for (uint32_t m : part.minimal_infeasible) + lo = std::min(lo, weight_activity(w, m)); + bound = hi; + return hi < lo; +} + +// The rewritten row must not enlarge this row's LP relaxation: every x in [0,1]^k with w.x <= t has +// to satisfy a.x <= rhs. With a > 0 and w >= 0 after N1 that is a fractional knapsack, max a.x +// subject to w.x <= t, solved by taking the zero-weight entries for free and then filling capacity +// in decreasing a_i/w_i order. Greedy leaves at most one fractional entry, so the comparison closes +// exactly in rationals and no tolerance is needed. The 0-1 knapsack DP in cuts/cuts.hpp cannot be +// substituted here: its optimum is a lower bound on the continuous one, which would let a weakening +// row through. +bool lp_no_weakening(const norm_row_t& row, const int64_t* w, int64_t t) +{ + cuopt_assert(t >= 0, "acceptance implies the origin is feasible, so the bound is non-negative"); + + std::array order{}; + int n_items = 0; + __int128 value = 0; + for (int i = 0; i < row.k; ++i) { + if (w[i] == 0) + value += row.coef[i]; + else + order[n_items++] = i; + } + std::sort(order.begin(), order.begin() + n_items, [&](int x, int y) { + const __int128 dx = (__int128)row.coef[x] * w[y]; + const __int128 dy = (__int128)row.coef[y] * w[x]; + return dx != dy ? dx > dy : x < y; + }); + + int64_t capacity = t; + int fractional = -1; + for (int p = 0; p < n_items; ++p) { + const int i = order[p]; + if (w[i] <= capacity) { + value += row.coef[i]; + capacity -= w[i]; + } else { + if (capacity > 0) fractional = i; + break; + } + } + + if (fractional < 0) return value <= (__int128)row.rhs; + return value * w[fractional] + (__int128)capacity * row.coef[fractional] <= + (__int128)row.rhs * w[fractional]; +} + +// Debug companion to accepts(): checks the 0/1 partition over every point rather than the extremal +// ones the search relies on. +[[maybe_unused]] bool verify_equivalent(const norm_row_t& row, const int64_t* w, int64_t t) +{ + const uint32_t n_pat = 1u << row.k; + for (uint32_t m = 0; m < n_pat; ++m) { + int64_t a_activity = 0; + int64_t w_activity = 0; + for (int i = 0; i < row.k; ++i) { + if ((m >> i & 1u) == 0u) continue; + a_activity += row.coef[i]; + w_activity += w[i]; + } + if ((a_activity <= row.rhs) != (w_activity <= t)) return false; + } + return true; +} + +struct search_state_t { + const norm_row_t* row = nullptr; + const partition_t* part = nullptr; + std::array w{}; + std::array best_w{}; + int64_t best_bound = 0; + int best_nonzeros = 0; + int64_t best_sum = 0; + bool found = false; +}; + +// Enumerates the non-increasing weight vectors reachable from the prefix fixed so far. N1 makes the +// weights non-negative and N3 makes them non-increasing, so a candidate is just a non-increasing +// sequence bounded by w[0]; Lemma 3.6's strict steps both shrink the branching factor and bound how +// low a position may go while still leaving room for the steps beneath it. +void search_positions(search_state_t& st, int pos) +{ + const int k = st.row->k; + if (pos == k) { + int64_t bound = 0; + if (!accepts(*st.part, st.w.data(), bound)) return; + if (!lp_no_weakening(*st.row, st.w.data(), bound)) return; + + int nonzeros = 0; + int64_t sum = 0; + for (int i = 0; i < k; ++i) { + nonzeros += st.w[i] != 0 ? 1 : 0; + sum += st.w[i]; + } + // Same max|w| by construction at this depth, so prefer dropping variables, then smaller + // weights. + if (st.found && + (nonzeros > st.best_nonzeros || (nonzeros == st.best_nonzeros && sum >= st.best_sum))) + return; + st.found = true; + st.best_nonzeros = nonzeros; + st.best_sum = sum; + st.best_bound = bound; + st.best_w = st.w; + return; + } + + const int64_t upper = st.w[pos - 1] - (st.part->strict[pos - 1] ? 1 : 0); + const int64_t lower = st.part->suffix_strict[pos]; + for (int64_t v = upper; v >= lower; --v) { + st.w[pos] = v; + search_positions(st, pos + 1); + } +} + +// Fallback for the rows whose smallest equivalent magnitude exceeds BHW_EXACT_MAX_WEIGHT, where the +// exhaustive search gives up: w = round(a / min a), which was the optimal scale in 131 of 131 +// measured cases, and the all-ones clause form. N3 already sorted a, so both are non-increasing. +bool heuristic_reduce(const norm_row_t& row, + const partition_t& part, + std::vector& weights, + int64_t& bound) +{ + const int k = row.k; + const int64_t a_min = row.coef[k - 1]; + cuopt_assert(a_min > 0, "N1 leaves every coefficient positive"); + + // Only a strict gain is worth installing: smaller magnitude, or the same magnitude with a + // variable dropped. + int64_t best_max = row.coef[0]; + int best_nonzeros = k; + bool found = false; + + std::array candidate{}; + for (int variant = 0; variant < 2; ++variant) { + for (int i = 0; i < k; ++i) + candidate[i] = variant == 0 ? (row.coef[i] + a_min / 2) / a_min : 1; + + int64_t candidate_bound = 0; + if (!accepts(part, candidate.data(), candidate_bound)) continue; + if (!lp_no_weakening(row, candidate.data(), candidate_bound)) continue; + + int64_t candidate_max = 0; + int nonzeros = 0; + for (int i = 0; i < k; ++i) { + candidate_max = std::max(candidate_max, candidate[i]); + nonzeros += candidate[i] != 0 ? 1 : 0; + } + if (candidate_max > best_max || (candidate_max == best_max && nonzeros >= best_nonzeros)) + continue; + + found = true; + best_max = candidate_max; + best_nonzeros = nonzeros; + weights.assign(candidate.begin(), candidate.begin() + k); + bound = candidate_bound; + } + return found; +} + +// Reduce one normalized shape. The caller undoes N1/N3 on the result. +bool reduce_shape(const norm_row_t& row, std::vector& weights, int64_t& bound) +{ + partition_t part; + if (!build_partition(row, part)) return false; + + const int64_t current = row.coef[0]; // N3 puts the largest coefficient first + cuopt_assert(current >= 2, "rows already at magnitude one are rejected before normalization"); + // Lemma 3.6 bounds max|w| from below over every equivalent inequality, so this row is provably + // irreducible in magnitude and not worth searching. + if (part.lemma36_bound >= current) return false; + + search_state_t st; + st.row = &row; + st.part = ∂ + const int64_t m_high = std::min(BHW_EXACT_MAX_WEIGHT, current - 1); + for (int64_t m = std::max(part.lemma36_bound, 1); m <= m_high; ++m) { + st.found = false; + st.w[0] = m; + search_positions(st, 1); + if (!st.found) continue; + // First m with any acceptance, so this is the minimum achievable max|w|. + weights.assign(st.best_w.begin(), st.best_w.begin() + row.k); + bound = st.best_bound; + return true; + } + return heuristic_reduce(row, part, weights, bound); +} + +} // namespace + +template +bhw_row_rewrite_t bhw_reduce_row( + const f_t* coefficients, int len, f_t side, int direction, bhw_shape_cache_t* cache) +{ + cuopt_assert(direction == 1 || direction == -1, + "direction is the sign that orients the row to <="); + bhw_row_rewrite_t rewrite; + if (len < 2 || len > BHW_MAX_LEN) return rewrite; + if (!scaling_bound_finite(side)) return rewrite; + + // Integerize so the point partition is exact, then orient the row to a.x <= b. + const double scale = row_int_scale( + coefficients, len, side, std::numeric_limits::infinity(), BHW_MAX_LEN, BHW_INT_SCALE_MAX); + if (scale == 0.0) return rewrite; + + std::array integral{}; + int64_t largest = 0; + for (int j = 0; j < len; ++j) { + integral[j] = (int64_t)std::llround((double)coefficients[j] * scale) * direction; + largest = std::max(largest, std::abs(integral[j])); + } + // A row already at +/-1 has no magnitude to give back. The census puts most of the corpus here, + // so this test carries the screening cost. + if (largest <= 1) return rewrite; + + norm_row_t norm_row; + norm_row.k = len; + norm_row.rhs = (int64_t)std::llround((double)side * scale) * direction; + std::array order{}; + for (int j = 0; j < len; ++j) { + order[j] = j; + // N1: complementing x_j = 1 - y_j moves the negative coefficient onto the right-hand side. + if (integral[j] < 0) norm_row.rhs -= integral[j]; + } + // N3: descending by magnitude, ties broken by position so the shape key is deterministic. + std::sort(order.begin(), order.begin() + len, [&](int x, int y) { + const int64_t ax = std::abs(integral[x]); + const int64_t ay = std::abs(integral[y]); + return ax != ay ? ax > ay : x < y; + }); + for (int p = 0; p < len; ++p) { + const int j = order[p]; + norm_row.coef[p] = std::abs(integral[j]); + norm_row.slot[p] = j; + norm_row.flipped[p] = integral[j] < 0; + } + + bhw_shape_result_t computed; + const bhw_shape_result_t* result = nullptr; + if (cache != nullptr) { + std::vector key(norm_row.coef.begin(), norm_row.coef.begin() + len); + key.push_back(norm_row.rhs); + auto cached = cache->find(key); + if (cached == cache->end()) { + bhw_shape_result_t fresh; + fresh.accepted = reduce_shape(norm_row, fresh.weights, fresh.bound); + cached = cache->emplace(std::move(key), std::move(fresh)).first; + } + result = &cached->second; + } else { + computed.accepted = reduce_shape(norm_row, computed.weights, computed.bound); + result = &computed; + } + if (!result->accepted) return rewrite; + + cuopt_assert((int)result->weights.size() == len, "cached shape has the wrong length"); + cuopt_assert(*std::min_element(result->weights.begin(), result->weights.end()) >= 0, + "N1 leaves the reduced weights non-negative"); + cuopt_assert( + *std::max_element(result->weights.begin(), result->weights.end()) == result->weights[0], + "N3 leaves the reduced weights non-increasing"); + cuopt_assert(result->weights[0] < norm_row.coef[0] || + std::count(result->weights.begin(), result->weights.end(), (int64_t)0) > 0, + "an accepted rewrite must shrink the magnitude or drop a variable"); + cuopt_assert(verify_equivalent(norm_row, result->weights.data(), result->bound), + "BHW rewrite changed the 0/1 feasible set"); + + // Undo N3 and N1, then undo the orientation. Complementing back turns w_i y_i into w_i - w_i x_i, + // which flips the coefficient and moves w_i onto the bound. + rewrite.coefficients.assign(len, 0); + int64_t new_side = result->bound; + for (int p = 0; p < len; ++p) { + const int j = norm_row.slot[p]; + if (norm_row.flipped[p]) { + rewrite.coefficients[j] = -result->weights[p]; + new_side -= result->weights[p]; + } else { + rewrite.coefficients[j] = result->weights[p]; + } + } + for (int j = 0; j < len; ++j) + rewrite.coefficients[j] *= direction; + rewrite.side = new_side * direction; + rewrite.max_coef_before = norm_row.coef[0]; + rewrite.max_coef_after = result->weights[0]; + rewrite.accepted = true; + return rewrite; +} +template +papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& problem, + const papilo::ProblemUpdate& problemUpdate, + const papilo::Num& num, + papilo::Reductions& reductions, + const papilo::Timer& timer, + int& reason_of_infeasibility) +{ + const auto& constraint_matrix = problem.getConstraintMatrix(); + const auto& lhs_values = constraint_matrix.getLeftHandSides(); + const auto& rhs_values = constraint_matrix.getRightHandSides(); + const auto& row_flags = constraint_matrix.getRowFlags(); + const auto& domains = problem.getVariableDomains(); + const auto& col_flags = domains.flags; + const auto& lower_bounds = domains.lower_bounds; + const auto& upper_bounds = domains.upper_bounds; + + const int num_rows = constraint_matrix.getNRows(); + papilo::PresolveStatus status = papilo::PresolveStatus::kUnchanged; + int64_t coefficients_reduced = 0; + int64_t coefficients_dropped = 0; + std::vector row_shrinks; + + // Every eligible row is screened, not only problemUpdate.getChangedActivities(): that worklist is + // seeded in full once and afterwards fed only by activity changes, so a row whose side changes is + // never revisited. Rescreening is cheap because the shape cache absorbs the repetition and the + // reduction is idempotent. + for (int row = 0; row < num_rows; ++row) { + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + const int len = row_coefficients.getLength(); + if (len < 2 || len > BHW_MAX_LEN) continue; + + const auto& row_flag = row_flags[row]; + if (row_flag.test(papilo::RowFlag::kRedundant)) continue; + const bool lhs_infinite = row_flag.test(papilo::RowFlag::kLhsInf); + const bool rhs_infinite = row_flag.test(papilo::RowFlag::kRhsInf); + // Equal flags mean either a ranged row / equation (both sides finite) or a free row. + if (lhs_infinite == rhs_infinite) continue; + + const int* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); + bool all_binary = true; + for (int j = 0; j < len && all_binary; ++j) { + const int col = indices[j]; + all_binary = col_flags[col].test(papilo::ColFlag::kIntegral) && + !col_flags[col].test(papilo::ColFlag::kLbInf) && + !col_flags[col].test(papilo::ColFlag::kUbInf) && + !col_flags[col].test(papilo::ColFlag::kFixed) && num.isZero(lower_bounds[col]) && + num.isEq(upper_bounds[col], f_t{1}); + } + if (!all_binary) continue; + + const int direction = lhs_infinite ? 1 : -1; + const f_t side = lhs_infinite ? rhs_values[row] : lhs_values[row]; + const bhw_row_rewrite_t rewrite = + bhw_reduce_row(values, len, side, direction, &shape_cache_); + if (!rewrite.accepted) continue; + + cuopt_assert(rewrite.max_coef_after >= 1, + "an accepted rewrite keeps at least one nonzero weight"); + row_shrinks.push_back((double)rewrite.max_coef_before / (double)rewrite.max_coef_after); + + // Same shape as papilo's own sparsifier, SimplifyInequalities: lock the row, then the entries, + // then the side. Dropping an entry needs no column lock -- ProblemUpdate marks the column + // modified and derives the resulting singleton rows and empty columns itself. The zero has to + // be exact: SparseStorage::changeRowInplace compacts an entry out on newval == 0, not on + // num.isZero. + papilo::TransactionGuard guard{reductions}; + reductions.lockRow(row); + [[maybe_unused]] int emitted = 0; + for (int j = 0; j < len; ++j) { + if ((f_t)rewrite.coefficients[j] == values[j]) continue; + reductions.changeMatrixEntry(row, indices[j], (f_t)rewrite.coefficients[j]); + ++emitted; + ++coefficients_reduced; + coefficients_dropped += rewrite.coefficients[j] == 0; + } + if (direction == 1) { + if ((f_t)rewrite.side != rhs_values[row]) { + reductions.changeRowRHS(row, (f_t)rewrite.side); + ++emitted; + } + } else { + if ((f_t)rewrite.side != lhs_values[row]) { + reductions.changeRowLHS(row, (f_t)rewrite.side); + ++emitted; + } + } + // Reporting kReduced for a transaction that changes nothing would have papilo re-derive the + // same rewrite every round. + cuopt_assert(emitted > 0, "accepted rewrite emitted no reduction"); + status = papilo::PresolveStatus::kReduced; + } + + if (coefficients_reduced > 0) { + const size_t n_rows_rewritten = row_shrinks.size(); + cuopt_assert(n_rows_rewritten > 0, "a changed coefficient implies an accepted row"); + const double mean = + std::accumulate(row_shrinks.begin(), row_shrinks.end(), 0.0) / (double)n_rows_rewritten; + // Mean alone is outlier-driven here: one row collapsing from a large coefficient to 1 outweighs + // hundreds of modest rewrites, so report the median next to it. + const auto middle = row_shrinks.begin() + n_rows_rewritten / 2; + std::nth_element(row_shrinks.begin(), middle, row_shrinks.end()); + double median = *middle; + if (n_rows_rewritten % 2 == 0) + median = (median + *std::max_element(row_shrinks.begin(), middle)) / 2.0; + + CUOPT_LOG_DEBUG( + "BHW reduced %ld coefficients (%ld dropped) in %zu rows, " + "max|a| shrank %.1fx mean, %.1fx median", + coefficients_reduced, + coefficients_dropped, + n_rows_rewritten, + mean, + median); + } + + return status; +} + +#define INSTANTIATE(F_TYPE) \ + template class BHWCoeffReduce; \ + template bhw_row_rewrite_t bhw_reduce_row( \ + const F_TYPE*, int, F_TYPE, int, bhw_shape_cache_t*); + +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT +INSTANTIATE(float) +#endif + +#if MIP_INSTANTIATE_DOUBLE +INSTANTIATE(double) +#endif + +#undef INSTANTIATE + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp new file mode 100644 index 0000000000..f01892b652 --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp @@ -0,0 +1,101 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#if !defined(__clang__) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wstringop-overflow" // ignore boost error for pip wheel build +#pragma GCC diagnostic ignored "-Wnarrowing" +#endif +#include +#include +#include +#include +#if !defined(__clang__) +#pragma GCC diagnostic pop +#endif + +#include +#include +#include + +namespace cuopt::mathematical_optimization::mip { + +// Widest row we enumerate: building the point partition walks 2^BHW_MAX_LEN patterns. A census of +// the PaPILO-reduced MIPLIB corpus (6.2M all-binary one-sided rows) puts 98.0% of the exactly +// integralizable ones at nnz <= 12; raising this to 40 would add 2.0%. +static constexpr int BHW_MAX_LEN = 12; +// Largest max|w| the exhaustive search considers before falling back to the heuristic candidates. +// Accepted weights were at most 3 throughout the corpus study. +static constexpr int64_t BHW_EXACT_MAX_WEIGHT = 6; +// Largest per-row rational multiplier / denominator used to integerize a row (passed to +// row_int_scale as its maxdnom/maxfinal caps). +static constexpr int64_t BHW_INT_SCALE_MAX = 1000000; // 1e6 + +// Outcome for one canonical row shape, in BHW's normalized frame (coefficients complemented to be +// positive and sorted descending). Rejections are cached too, since re-deriving them is the bulk of +// the work on instances whose rows repeat. +struct bhw_shape_result_t { + std::vector weights; + int64_t bound = 0; + bool accepted = false; +}; + +// Keyed by the normalized coefficients followed by the normalized right-hand side. The reduction is +// a pure function of that key, so entries stay valid across presolve rounds and problems. +using bhw_shape_cache_t = std::map, bhw_shape_result_t>; + +struct bhw_row_rewrite_t { + std::vector coefficients; // one per input entry, in input order; 0 drops that entry + int64_t side = 0; // replaces the row's finite side + // Largest |coefficient| before and after reduction, both in the integerized frame. Comparable + // only there: the row is scaled on the way in, so the input coefficients sit in a different + // frame. + int64_t max_coef_before = 0; + int64_t max_coef_after = 0; + bool accepted = false; +}; + +// Rewrite one one-sided all-binary row with smaller integer coefficients spanning the same 0/1 +// feasible set. direction is +1 for "coefficients . x <= side" and -1 for ">= side"; side is the +// finite side of the row. Rejects the row (accepted = false) unless it integerizes exactly, admits +// a strictly smaller equivalent form, and that form does not enlarge the row's LP relaxation. cache +// may be null to skip memoization. +// +// The caller checks that every entry is a binary integer variable and that exactly one side of the +// row is finite. Exposed for testing: BHWCoeffReduce::execute only screens rows and emits the +// result, so this covers the whole reduction without any papilo types. +template +bhw_row_rewrite_t bhw_reduce_row( + const f_t* coefficients, int len, f_t side, int direction, bhw_shape_cache_t* cache); + +// Bradley-Hammer-Wolsey coefficient reduction: replace an all-binary row by an equivalent one with +// smaller integer coefficients. See bhw_coeff_reduce.cpp for the lineage. +template +class BHWCoeffReduce : public papilo::PresolveMethod { + public: + BHWCoeffReduce() : papilo::PresolveMethod() + { + this->setName("bhwcoeffreduce"); + this->setType(papilo::PresolverType::kIntegralCols); + this->setTiming(papilo::PresolverTiming::kMedium); + } + + papilo::PresolveStatus execute(const papilo::Problem& problem, + const papilo::ProblemUpdate& problemUpdate, + const papilo::Num& num, + papilo::Reductions& reductions, + const papilo::Timer& timer, + int& reason_of_infeasibility) override; + + private: + // Only touched from execute, which papilo runs one task at a time per presolver object. + bhw_shape_cache_t shape_cache_; +}; + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 3d5b165abb..10784f7021 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -40,6 +40,7 @@ #include #include #include +#include #include #include #include @@ -677,6 +678,7 @@ void set_presolve_methods( if (category == problem_category_t::MIP) { // cuOpt custom GF2 presolver maybe_add(uptr(new cuopt::mathematical_optimization::mip::GF2Presolve())); + maybe_add(uptr(new cuopt::mathematical_optimization::mip::BHWCoeffReduce())); } // fast presolvers maybe_add(uptr(new papilo::SingletonCols())); diff --git a/cpp/tests/internal/CMakeLists.txt b/cpp/tests/internal/CMakeLists.txt index 9b5cd2ecb7..99f26aa440 100644 --- a/cpp/tests/internal/CMakeLists.txt +++ b/cpp/tests/internal/CMakeLists.txt @@ -28,6 +28,7 @@ ConfigureTest(NUMOPT_INTERNAL_TEST ${CUOPT_TEST_DIR}/mip/empty_fixed_problems_test.cu ${CUOPT_TEST_DIR}/mip/presolve_test.cu ${CUOPT_TEST_DIR}/mip/block_bve_test.cu + ${CUOPT_TEST_DIR}/mip/bhw_coeff_reduce_test.cpp ${CUOPT_TEST_DIR}/mip/gf2_presolve_test.cpp ${CUOPT_TEST_DIR}/mip/termination_test.cu ${CUOPT_TEST_DIR}/mip/determinism_test.cu diff --git a/cpp/tests/mip/bhw_coeff_reduce_test.cpp b/cpp/tests/mip/bhw_coeff_reduce_test.cpp new file mode 100644 index 0000000000..432bf2d078 --- /dev/null +++ b/cpp/tests/mip/bhw_coeff_reduce_test.cpp @@ -0,0 +1,171 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include + +#include + +#include +#include +#include +#include + +namespace cuopt::mathematical_optimization::test { + +using mip::BHW_MAX_LEN; +using mip::bhw_reduce_row; +using mip::bhw_row_rewrite_t; +using mip::bhw_shape_cache_t; + +namespace { + +bhw_row_rewrite_t reduce(const std::vector& coefficients, + double side, + int direction = 1, + bhw_shape_cache_t* cache = nullptr) +{ + return bhw_reduce_row( + coefficients.data(), (int)coefficients.size(), side, direction, cache); +} + +// The whole point of the pass: the rewritten row must accept exactly the same 0/1 points as the +// original. Checked over every point, independently of the extremal-point test the search uses. +bool same_feasible_set(const std::vector& coefficients, + double side, + int direction, + const bhw_row_rewrite_t& rewrite) +{ + constexpr double tol = 1e-9; + const int k = (int)coefficients.size(); + for (uint32_t mask = 0; mask < (1u << k); ++mask) { + double original = 0.0; + int64_t rewritten = 0; + for (int i = 0; i < k; ++i) { + if ((mask >> i & 1u) == 0u) continue; + original += coefficients[i]; + rewritten += rewrite.coefficients[i]; + } + const bool original_ok = direction == 1 ? original <= side + tol : original >= side - tol; + const bool rewritten_ok = + direction == 1 ? rewritten <= rewrite.side : rewritten >= rewrite.side; + if (original_ok != rewritten_ok) return false; + } + return true; +} + +} // namespace + +// Bradley, Hammer and Wolsey (1974) open with this row and reduce it to 4,4,2,2,1,1,1,0 <= 5. That +// rewrite enlarges the LP relaxation -- it admits a fractional point of activity 91.75 against a +// right-hand side of 80 -- so the LP-strength check rejects it and no smaller equivalent form +// survives. +TEST(bhw_coeff_reduce, rejects_a_rewrite_that_weakens_the_relaxation) +{ + EXPECT_FALSE(reduce({65, 64, 41, 22, 13, 12, 8, 2}, 80).accepted); +} + +// The shape that motivated the pass: bnatt400 leaves rows whose coefficients are 9/10 and 1/3 +// against a right-hand side of 2/3. Integerizing and reducing lands them in int8 range. +TEST(bhw_coeff_reduce, rational_row_integerizes_and_reduces) +{ + const std::vector row{0.9, 1.0 / 3, 1.0 / 3, 1.0 / 3}; + const auto reduced = reduce(row, 2.0 / 3); + ASSERT_TRUE(reduced.accepted); + EXPECT_EQ(reduced.coefficients, std::vector({3, 1, 1, 1})); + EXPECT_EQ(reduced.side, 2); + EXPECT_TRUE(same_feasible_set(row, 2.0 / 3, 1, reduced)); +} + +// The >= orientation is normalized by negation, so the same row negated must come back negated. +TEST(bhw_coeff_reduce, greater_equal_row_keeps_its_orientation) +{ + const std::vector row{-0.9, -1.0 / 3, -1.0 / 3, -1.0 / 3}; + const auto reduced = reduce(row, -2.0 / 3, -1); + ASSERT_TRUE(reduced.accepted); + EXPECT_EQ(reduced.coefficients, std::vector({-3, -1, -1, -1})); + EXPECT_EQ(reduced.side, -2); + EXPECT_TRUE(same_feasible_set(row, -2.0 / 3, -1, reduced)); +} + +TEST(bhw_coeff_reduce, rejects_rows_with_nothing_to_give_back) +{ + // Already at unit magnitude. + EXPECT_FALSE(reduce({1, 1, 1, 1}, 2).accepted); + // Does not integerize within the rational cap. + EXPECT_FALSE(reduce({M_PI, 1, 1}, 2).accepted); + // Outside the enumerable width. + EXPECT_FALSE(reduce({5}, 2).accepted); + EXPECT_FALSE(reduce(std::vector(BHW_MAX_LEN + 1, 3.0), 5).accepted); + // Every point feasible, so there is nothing to separate. + EXPECT_FALSE(reduce({3, 2, 2}, 100).accepted); + // No point feasible. + EXPECT_FALSE(reduce({3, 2, 2}, -1).accepted); +} + +TEST(bhw_coeff_reduce, memoized_result_matches_the_uncached_one) +{ + bhw_shape_cache_t cache; + const std::vector> rows{ + {0.9, 1.0 / 3, 1.0 / 3, 1.0 / 3}, {6, 4, 3, 2}, {-6, 4, 3, -2}, {9, 7, 6, 6, 4}}; + for (const auto& row : rows) { + for (int repeat = 0; repeat < 2; ++repeat) { + const auto cached = reduce(row, 12, 1, &cache); + const auto uncached = reduce(row, 12, 1, nullptr); + EXPECT_EQ(cached.accepted, uncached.accepted); + EXPECT_EQ(cached.coefficients, uncached.coefficients); + EXPECT_EQ(cached.side, uncached.side); + } + } +} + +// The invariant that matters, over mixed signs, both orientations and rational coefficients: an +// accepted rewrite never changes which 0/1 points satisfy the row. +TEST(bhw_coeff_reduce, accepted_rewrites_preserve_the_feasible_set) +{ + std::mt19937_64 rng(20260805); + bhw_shape_cache_t cache; + const int denominators[] = {1, 2, 3, 4, 5, 6, 8, 10, 12, 16}; + int accepted = 0; + + for (int trial = 0; trial < 20000; ++trial) { + const int len = 2 + (int)(rng() % 7); + const int direction = (rng() & 1u) != 0u ? 1 : -1; + const int denominator = denominators[rng() % 10]; + + std::vector row(len); + double positive_sum = 0.0; + double negative_sum = 0.0; + for (int i = 0; i < len; ++i) { + const int64_t numerator = 1 + (int64_t)(rng() % 30); + row[i] = (double)numerator / denominator * ((rng() & 3u) == 0u ? -1.0 : 1.0); + if (row[i] > 0.0) + positive_sum += row[i]; + else + negative_sum += row[i]; + } + // Put the side inside the activity range so the row is not trivially satisfied or violated. + double side = negative_sum + (positive_sum - negative_sum) * (double)(rng() % 1001) / 1000.0; + side = std::round(side * denominator) / denominator; + if (direction == -1) { + for (double& value : row) + value = -value; + side = -side; + } + + const auto reduced = reduce(row, side, direction, &cache); + if (!reduced.accepted) continue; + ++accepted; + + ASSERT_EQ((int)reduced.coefficients.size(), len); + ASSERT_TRUE(same_feasible_set(row, side, direction, reduced)) + << "rewrite changed the 0/1 feasible set on trial " << trial; + } + // Guards against the generator drifting into a corner where nothing is ever reduced. + EXPECT_GT(accepted, 1000); +} + +} // namespace cuopt::mathematical_optimization::test From 880064c009ed22c43e7986c3e4016285028fb887 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Thu, 6 Aug 2026 08:19:10 -0700 Subject: [PATCH 2/8] cleanup --- cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp | 12 +++++++++++- cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp | 6 +++--- cpp/tests/mip/bhw_coeff_reduce_test.cpp | 7 +++++++ 3 files changed, 21 insertions(+), 4 deletions(-) diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp index 862a5b61cb..5685018d56 100644 --- a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp @@ -365,7 +365,8 @@ bhw_row_rewrite_t bhw_reduce_row( int64_t largest = 0; for (int j = 0; j < len; ++j) { integral[j] = (int64_t)std::llround((double)coefficients[j] * scale) * direction; - largest = std::max(largest, std::abs(integral[j])); + if (integral[j] == 0) return rewrite; + largest = std::max(largest, std::abs(integral[j])); } // A row already at +/-1 has no magnitude to give back. The census puts most of the corpus here, // so this test carries the screening cost. @@ -444,6 +445,8 @@ bhw_row_rewrite_t bhw_reduce_row( rewrite.accepted = true; return rewrite; } +static constexpr int BHW_INTERRUPT_CHECK_STRIDE = 256; + template papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& problem, const papilo::ProblemUpdate& problemUpdate, @@ -460,6 +463,7 @@ papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& const auto& col_flags = domains.flags; const auto& lower_bounds = domains.lower_bounds; const auto& upper_bounds = domains.upper_bounds; + const auto& presolve_options = problemUpdate.getPresolveOptions(); const int num_rows = constraint_matrix.getNRows(); papilo::PresolveStatus status = papilo::PresolveStatus::kUnchanged; @@ -472,6 +476,12 @@ papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& // never revisited. Rescreening is cheap because the shape cache absorbs the repetition and the // reduction is idempotent. for (int row = 0; row < num_rows; ++row) { + if (reductions.size() >= presolve_options.max_reduction_seq) break; + if (row % BHW_INTERRUPT_CHECK_STRIDE == 0 && + papilo::PresolveMethod::is_interrupted( + timer, presolve_options.tlim, presolve_options.early_exit_callback)) + break; + auto row_coefficients = constraint_matrix.getRowCoefficients(row); const int len = row_coefficients.getLength(); if (len < 2 || len > BHW_MAX_LEN) continue; diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp index f01892b652..d62696f629 100644 --- a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp @@ -63,9 +63,9 @@ struct bhw_row_rewrite_t { // Rewrite one one-sided all-binary row with smaller integer coefficients spanning the same 0/1 // feasible set. direction is +1 for "coefficients . x <= side" and -1 for ">= side"; side is the -// finite side of the row. Rejects the row (accepted = false) unless it integerizes exactly, admits -// a strictly smaller equivalent form, and that form does not enlarge the row's LP relaxation. cache -// may be null to skip memoization. +// finite side of the row. Rejects the row (accepted = false) unless it integerizes exactly to +// nonzero coefficients, admits a strictly smaller equivalent form, and that form does not enlarge +// the row's LP relaxation. cache may be null to skip memoization. // // The caller checks that every entry is a binary integer variable and that exactly one side of the // row is finite. Exposed for testing: BHWCoeffReduce::execute only screens rows and emits the diff --git a/cpp/tests/mip/bhw_coeff_reduce_test.cpp b/cpp/tests/mip/bhw_coeff_reduce_test.cpp index 432bf2d078..7dd6769b6d 100644 --- a/cpp/tests/mip/bhw_coeff_reduce_test.cpp +++ b/cpp/tests/mip/bhw_coeff_reduce_test.cpp @@ -106,6 +106,13 @@ TEST(bhw_coeff_reduce, rejects_rows_with_nothing_to_give_back) EXPECT_FALSE(reduce({3, 2, 2}, -1).accepted); } +TEST(bhw_coeff_reduce, rejects_rows_with_a_zero_coefficient) +{ + EXPECT_FALSE(reduce({65, 64, 41, 22, 13, 12, 8, 2, 0}, 80).accepted); + EXPECT_FALSE(reduce({0, 9, 7, 6, 6, 4}, 20).accepted); + EXPECT_FALSE(reduce({6, 0}, 5).accepted); +} + TEST(bhw_coeff_reduce, memoized_result_matches_the_uncached_one) { bhw_shape_cache_t cache; From ad0a45c6801b2563cea99964cf7c91249a9c0f9a Mon Sep 17 00:00:00 2001 From: yboucher Date: Wed, 19 Aug 2026 06:34:14 -0700 Subject: [PATCH 3/8] cleanup --- .../presolve/bhw_coeff_reduce.cpp | 124 +++++++++++------- .../presolve/bhw_coeff_reduce.hpp | 8 +- cpp/tests/mip/bhw_coeff_reduce_test.cpp | 8 +- 3 files changed, 84 insertions(+), 56 deletions(-) diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp index 5685018d56..4d6de0c8d4 100644 --- a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp @@ -31,29 +31,28 @@ // points, which are supersets of BHW's ceilings and roofs because we skip their ordering condition. // That costs nothing: with all coefficients positive the activity is monotone, so the maximum of // w.x over the feasible points is attained at a maximal one and the minimum over the infeasible -// points at a minimal one. Checked exhaustively over 9.4M (row, weight) pairs against full 2^k -// equivalence, in both directions, before landing. +// points at a minimal one. // // BHW's minimization is not implemented. Section 5 obtains the minimum equivalent inequality by LP // over the polytope. We instead search weight vectors by increasing max|w| up to a cap, which is // minimal by construction below the cap, and fall back to two heuristic candidates above it. That // search is affordable at every width we enumerate only because conditions N1 and N3 restrict // candidates to non-negative non-increasing sequences: at most sum_{m=1..6} C(11+m, m) = 18563 -// vectors for a 12-entry row, measured at 187us for the worst shape. Lemma 3.6 supplies the lower -// bound that seeds and prunes the search. +// vectors for a 12-entry row. Lemma 3.6 supplies the lower bound that seeds and prunes the search. // -// TODO: extend the fallback path to the row generation of Section 6. On the PaPILO-reduced corpus -// the LP reaches 1489 reducible rows where these heuristics match the optimum on roughly 60%, and -// unlocks int8 on 321 of them. The separation oracle in step 3 is a 0-1 knapsack, max w.x subject -// to a.x <= b, which solve_knapsack_problem in cuts/cuts.hpp already solves by DP. Note that -// routine is the wrong tool for the LP-strength gate below, which needs the continuous relaxation. -// Boyd (1993) Generating Fenchel Cutting Planes for Knapsack Polyhedra, SIAM J. Optim. -// 3(4):734-750, treats the same separation problem as cut generation. +// TODO: extend the fallback path to the row generation of Section 6, which reaches rows the two +// heuristics below leave at a larger magnitude than necessary. The separation oracle in step 3 is a +// 0-1 knapsack, max w.x subject to a.x <= b, which solve_knapsack_problem in cuts/cuts.hpp already +// solves by DP. Boyd (1993) Generating Fenchel Cutting Planes for Knapsack Polyhedra, SIAM J. +// Optim. 3(4):734-750, treats the same separation problem as cut generation. namespace cuopt::mathematical_optimization::mip { namespace { +// Rows between two interrupt checks in execute; screening one row is cheap next to the check. +constexpr int BHW_INTERRUPT_CHECK_STRIDE = 256; + // One row in BHW's normalized frame: condition N1 (every coefficient positive, negatives // complemented by x_i = 1 - y_i) and condition N3 (coefficients sorted descending). N1 is what // makes the row activity monotone in x, and it also makes the weights of any equivalent inequality @@ -271,8 +270,9 @@ void search_positions(search_state_t& st, int pos) } // Fallback for the rows whose smallest equivalent magnitude exceeds BHW_EXACT_MAX_WEIGHT, where the -// exhaustive search gives up: w = round(a / min a), which was the optimal scale in 131 of 131 -// measured cases, and the all-ones clause form. N3 already sorted a, so both are non-increasing. +// exhaustive search gives up: w = round(a / min a) and the all-ones clause form, both put through +// the same acceptance and LP-strength gates as a searched vector. N3 already sorted a, so both +// candidates are non-increasing. bool heuristic_reduce(const norm_row_t& row, const partition_t& part, std::vector& weights, @@ -364,17 +364,17 @@ bhw_row_rewrite_t bhw_reduce_row( std::array integral{}; int64_t largest = 0; for (int j = 0; j < len; ++j) { - integral[j] = (int64_t)std::llround((double)coefficients[j] * scale) * direction; + integral[j] = std::llround((double)coefficients[j] * scale) * direction; if (integral[j] == 0) return rewrite; largest = std::max(largest, std::abs(integral[j])); } - // A row already at +/-1 has no magnitude to give back. The census puts most of the corpus here, - // so this test carries the screening cost. + // A row already at +/-1 has no magnitude to give back; rejecting it here is what keeps screening + // cheap on the rows that dominate a model. if (largest <= 1) return rewrite; norm_row_t norm_row; norm_row.k = len; - norm_row.rhs = (int64_t)std::llround((double)side * scale) * direction; + norm_row.rhs = std::llround((double)side * scale) * direction; std::array order{}; for (int j = 0; j < len; ++j) { order[j] = j; @@ -419,7 +419,7 @@ bhw_row_rewrite_t bhw_reduce_row( *std::max_element(result->weights.begin(), result->weights.end()) == result->weights[0], "N3 leaves the reduced weights non-increasing"); cuopt_assert(result->weights[0] < norm_row.coef[0] || - std::count(result->weights.begin(), result->weights.end(), (int64_t)0) > 0, + std::count(result->weights.begin(), result->weights.end(), 0) > 0, "an accepted rewrite must shrink the magnitude or drop a variable"); cuopt_assert(verify_equivalent(norm_row, result->weights.data(), result->bound), "BHW rewrite changed the 0/1 feasible set"); @@ -445,7 +445,59 @@ bhw_row_rewrite_t bhw_reduce_row( rewrite.accepted = true; return rewrite; } -static constexpr int BHW_INTERRUPT_CHECK_STRIDE = 256; + +namespace { + +// Coefficient-shrink figures behind the DEBUG line. Both the accumulation and the summary compile +// away below DEBUG, so a release build carries neither the per-row vector nor the selection. +struct bhw_stats_t { +#if (CUOPT_LOG_ACTIVE_LEVEL <= RAPIDS_LOGGER_LOG_LEVEL_DEBUG) + int64_t coefficients_reduced = 0; + int64_t coefficients_dropped = 0; + std::vector row_shrinks; + + void changed_coefficient(int64_t new_coefficient) + { + ++coefficients_reduced; + coefficients_dropped += new_coefficient == 0; + } + + void rewrote_row(int64_t max_coef_before, int64_t max_coef_after) + { + row_shrinks.push_back((double)max_coef_before / max_coef_after); + } + + void report() + { + if (coefficients_reduced == 0) return; + const size_t n_rows_rewritten = row_shrinks.size(); + cuopt_assert(n_rows_rewritten > 0, "a changed coefficient implies an accepted row"); + const double mean = + std::accumulate(row_shrinks.begin(), row_shrinks.end(), 0.0) / n_rows_rewritten; + // One row collapsing to magnitude 1 dominates the mean, so report the median next to it. + const auto middle = row_shrinks.begin() + n_rows_rewritten / 2; + std::nth_element(row_shrinks.begin(), middle, row_shrinks.end()); + double median = *middle; + if (n_rows_rewritten % 2 == 0) + median = (median + *std::max_element(row_shrinks.begin(), middle)) / 2.0; + + CUOPT_LOG_DEBUG( + "BHW reduced %ld coefficients (%ld dropped) in %zu rows, " + "max|a| shrank %.1fx mean, %.1fx median", + coefficients_reduced, + coefficients_dropped, + n_rows_rewritten, + mean, + median); + } +#else + void changed_coefficient(int64_t) {} + void rewrote_row(int64_t, int64_t) {} + void report() {} +#endif +}; + +} // namespace template papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& problem, @@ -467,9 +519,7 @@ papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& const int num_rows = constraint_matrix.getNRows(); papilo::PresolveStatus status = papilo::PresolveStatus::kUnchanged; - int64_t coefficients_reduced = 0; - int64_t coefficients_dropped = 0; - std::vector row_shrinks; + bhw_stats_t stats; // Every eligible row is screened, not only problemUpdate.getChangedActivities(): that worklist is // seeded in full once and afterwards fed only by activity changes, so a row whose side changes is @@ -514,10 +564,10 @@ papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& cuopt_assert(rewrite.max_coef_after >= 1, "an accepted rewrite keeps at least one nonzero weight"); - row_shrinks.push_back((double)rewrite.max_coef_before / (double)rewrite.max_coef_after); + stats.rewrote_row(rewrite.max_coef_before, rewrite.max_coef_after); // Same shape as papilo's own sparsifier, SimplifyInequalities: lock the row, then the entries, - // then the side. Dropping an entry needs no column lock -- ProblemUpdate marks the column + // then the side. Dropping an entry needs no column lock: ProblemUpdate marks the column // modified and derives the resulting singleton rows and empty columns itself. The zero has to // be exact: SparseStorage::changeRowInplace compacts an entry out on newval == 0, not on // num.isZero. @@ -528,8 +578,7 @@ papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& if ((f_t)rewrite.coefficients[j] == values[j]) continue; reductions.changeMatrixEntry(row, indices[j], (f_t)rewrite.coefficients[j]); ++emitted; - ++coefficients_reduced; - coefficients_dropped += rewrite.coefficients[j] == 0; + stats.changed_coefficient(rewrite.coefficients[j]); } if (direction == 1) { if ((f_t)rewrite.side != rhs_values[row]) { @@ -548,28 +597,7 @@ papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& status = papilo::PresolveStatus::kReduced; } - if (coefficients_reduced > 0) { - const size_t n_rows_rewritten = row_shrinks.size(); - cuopt_assert(n_rows_rewritten > 0, "a changed coefficient implies an accepted row"); - const double mean = - std::accumulate(row_shrinks.begin(), row_shrinks.end(), 0.0) / (double)n_rows_rewritten; - // Mean alone is outlier-driven here: one row collapsing from a large coefficient to 1 outweighs - // hundreds of modest rewrites, so report the median next to it. - const auto middle = row_shrinks.begin() + n_rows_rewritten / 2; - std::nth_element(row_shrinks.begin(), middle, row_shrinks.end()); - double median = *middle; - if (n_rows_rewritten % 2 == 0) - median = (median + *std::max_element(row_shrinks.begin(), middle)) / 2.0; - - CUOPT_LOG_DEBUG( - "BHW reduced %ld coefficients (%ld dropped) in %zu rows, " - "max|a| shrank %.1fx mean, %.1fx median", - coefficients_reduced, - coefficients_dropped, - n_rows_rewritten, - mean, - median); - } + stats.report(); return status; } diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp index d62696f629..fee1d27334 100644 --- a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp @@ -26,12 +26,12 @@ namespace cuopt::mathematical_optimization::mip { -// Widest row we enumerate: building the point partition walks 2^BHW_MAX_LEN patterns. A census of -// the PaPILO-reduced MIPLIB corpus (6.2M all-binary one-sided rows) puts 98.0% of the exactly -// integralizable ones at nnz <= 12; raising this to 40 would add 2.0%. +// Widest row we enumerate: building the point partition walks 2^BHW_MAX_LEN patterns, so this caps +// the cost of screening one row. Wider rows are left alone; nearly every row that integerizes +// exactly is narrower than this. static constexpr int BHW_MAX_LEN = 12; // Largest max|w| the exhaustive search considers before falling back to the heuristic candidates. -// Accepted weights were at most 3 throughout the corpus study. +// Reductions needing a larger magnitude are rare enough not to pay for the extra search depth. static constexpr int64_t BHW_EXACT_MAX_WEIGHT = 6; // Largest per-row rational multiplier / denominator used to integerize a row (passed to // row_int_scale as its maxdnom/maxfinal caps). diff --git a/cpp/tests/mip/bhw_coeff_reduce_test.cpp b/cpp/tests/mip/bhw_coeff_reduce_test.cpp index 7dd6769b6d..20b8fdb0f6 100644 --- a/cpp/tests/mip/bhw_coeff_reduce_test.cpp +++ b/cpp/tests/mip/bhw_coeff_reduce_test.cpp @@ -60,16 +60,16 @@ bool same_feasible_set(const std::vector& coefficients, } // namespace // Bradley, Hammer and Wolsey (1974) open with this row and reduce it to 4,4,2,2,1,1,1,0 <= 5. That -// rewrite enlarges the LP relaxation -- it admits a fractional point of activity 91.75 against a -// right-hand side of 80 -- so the LP-strength check rejects it and no smaller equivalent form +// rewrite enlarges the LP relaxation (it admits a fractional point of activity 91.75 against a +// right-hand side of 80), so the LP-strength check rejects it and no smaller equivalent form // survives. TEST(bhw_coeff_reduce, rejects_a_rewrite_that_weakens_the_relaxation) { EXPECT_FALSE(reduce({65, 64, 41, 22, 13, 12, 8, 2}, 80).accepted); } -// The shape that motivated the pass: bnatt400 leaves rows whose coefficients are 9/10 and 1/3 -// against a right-hand side of 2/3. Integerizing and reducing lands them in int8 range. +// The shape that motivated the pass: a 9/10 coefficient against thirds, with a fractional +// right-hand side. Integerizing and reducing lands it in int8 range. TEST(bhw_coeff_reduce, rational_row_integerizes_and_reduces) { const std::vector row{0.9, 1.0 / 3, 1.0 / 3, 1.0 / 3}; From 03011ffc0db19fb6e1a36c16dcae8e6c266c5052 Mon Sep 17 00:00:00 2001 From: yboucher Date: Fri, 28 Aug 2026 00:43:44 -0700 Subject: [PATCH 4/8] fix numerical vulnerability --- .../presolve/bhw_coeff_reduce.cpp | 47 ++++++++++++++++--- cpp/tests/mip/bhw_coeff_reduce_test.cpp | 11 +++++ 2 files changed, 51 insertions(+), 7 deletions(-) diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp index 4d6de0c8d4..109ac182d3 100644 --- a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp @@ -50,9 +50,6 @@ namespace cuopt::mathematical_optimization::mip { namespace { -// Rows between two interrupt checks in execute; screening one row is cheap next to the check. -constexpr int BHW_INTERRUPT_CHECK_STRIDE = 256; - // One row in BHW's normalized frame: condition N1 (every coefficient positive, negatives // complemented by x_i = 1 - y_i) and condition N3 (coefficients sorted descending). N1 is what // makes the row activity monotone in x, and it also makes the weights of any equivalent inequality @@ -86,6 +83,36 @@ int64_t weight_activity(const int64_t* w, uint32_t mask) return sum; } +template +bool integerization_preserves_binary_feasible_set( + const f_t* coefficients, + int len, + f_t side, + int direction, + const std::array& integral, + int64_t integral_side) +{ + cuopt_assert(len >= 2 && len <= BHW_MAX_LEN, "row length outside the enumerable range"); + const uint32_t n_pat = 1u << len; + // software implemented, but addition is very cheap + std::vector<_Float128> original_activity(n_pat, 0.0L); + std::vector integral_activity(n_pat, 0); + for (uint32_t m = 1; m < n_pat; ++m) { + const uint32_t previous = m & (m - 1u); + const int j = std::countr_zero(m); + original_activity[m] = original_activity[previous] + coefficients[j]; + integral_activity[m] = integral_activity[previous] + integral[j]; + } + + for (uint32_t m = 0; m < n_pat; ++m) { + const bool original_feasible = direction == 1 ? original_activity[m] <= side + : original_activity[m] >= side; + const bool integral_feasible = integral_activity[m] <= integral_side; + if (original_feasible != integral_feasible) return false; + } + return true; +} + // Splits {0,1}^k and collects the maximal feasible and minimal infeasible points. Returns false for // a degenerate row (all points feasible or all infeasible), which is left to other presolvers. bool build_partition(const norm_row_t& row, partition_t& out) @@ -373,8 +400,9 @@ bhw_row_rewrite_t bhw_reduce_row( if (largest <= 1) return rewrite; norm_row_t norm_row; - norm_row.k = len; - norm_row.rhs = std::llround((double)side * scale) * direction; + norm_row.k = len; + const int64_t integral_side = std::llround((double)side * scale) * direction; + norm_row.rhs = integral_side; std::array order{}; for (int j = 0; j < len; ++j) { order[j] = j; @@ -412,6 +440,11 @@ bhw_row_rewrite_t bhw_reduce_row( } if (!result->accepted) return rewrite; + // Rows sharing an integerized cache key can have different original floating-point partitions. + if (!integerization_preserves_binary_feasible_set( + coefficients, len, side, direction, integral, integral_side)) + return rewrite; + cuopt_assert((int)result->weights.size() == len, "cached shape has the wrong length"); cuopt_assert(*std::min_element(result->weights.begin(), result->weights.end()) >= 0, "N1 leaves the reduced weights non-negative"); @@ -527,8 +560,8 @@ papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& // reduction is idempotent. for (int row = 0; row < num_rows; ++row) { if (reductions.size() >= presolve_options.max_reduction_seq) break; - if (row % BHW_INTERRUPT_CHECK_STRIDE == 0 && - papilo::PresolveMethod::is_interrupted( + // Screening one row is cheap next to the interrupt check. + if (papilo::PresolveMethod::is_interrupted( timer, presolve_options.tlim, presolve_options.early_exit_callback)) break; diff --git a/cpp/tests/mip/bhw_coeff_reduce_test.cpp b/cpp/tests/mip/bhw_coeff_reduce_test.cpp index 20b8fdb0f6..2bd7dbf4bf 100644 --- a/cpp/tests/mip/bhw_coeff_reduce_test.cpp +++ b/cpp/tests/mip/bhw_coeff_reduce_test.cpp @@ -80,6 +80,17 @@ TEST(bhw_coeff_reduce, rational_row_integerizes_and_reduces) EXPECT_TRUE(same_feasible_set(row, 2.0 / 3, 1, reduced)); } +TEST(bhw_coeff_reduce, rejects_approximate_integerization_that_changes_the_feasible_set) +{ + constexpr int num_variables = 12; + constexpr double perturbation = 9.9e-7; + constexpr double small_coefficient = 1000.0; + std::vector row(num_variables, small_coefficient + perturbation); + row[0] = 12000.0; + + EXPECT_FALSE(reduce(row, 11000.0 - perturbation).accepted); +} + // The >= orientation is normalized by negation, so the same row negated must come back negated. TEST(bhw_coeff_reduce, greater_equal_row_keeps_its_orientation) { From 5a8bd2df6fdfff699b8bb42f11ee7d3f2e609b84 Mon Sep 17 00:00:00 2001 From: yboucher Date: Fri, 28 Aug 2026 00:57:46 -0700 Subject: [PATCH 5/8] some tidying-up --- .../presolve/bhw_coeff_reduce.cpp | 130 +++++------------- .../presolve/bhw_coeff_reduce.hpp | 37 ++--- cpp/tests/mip/bhw_coeff_reduce_test.cpp | 22 +-- 3 files changed, 46 insertions(+), 143 deletions(-) diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp index 109ac182d3..55335a2d97 100644 --- a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp @@ -18,43 +18,18 @@ #include #include #include -#include #include #include -// Gordon H. Bradley, Peter L. Hammer, Laurence Wolsey (1974) Coefficient reduction for inequalities -// in 0-1 variables. Mathematical Programming 7:263-282. -// -// Over binaries many inequalities share the same 0/1 feasible set. Theorem 2.5 characterizes that -// set of "equivalent" inequalities by a system over the ceilings and roofs of the row. The test -// here is the same characterization taken over the plain maximal feasible and minimal infeasible -// points, which are supersets of BHW's ceilings and roofs because we skip their ordering condition. -// That costs nothing: with all coefficients positive the activity is monotone, so the maximum of -// w.x over the feasible points is attained at a maximal one and the minimum over the infeasible -// points at a minimal one. -// -// BHW's minimization is not implemented. Section 5 obtains the minimum equivalent inequality by LP -// over the polytope. We instead search weight vectors by increasing max|w| up to a cap, which is -// minimal by construction below the cap, and fall back to two heuristic candidates above it. That -// search is affordable at every width we enumerate only because conditions N1 and N3 restrict -// candidates to non-negative non-increasing sequences: at most sum_{m=1..6} C(11+m, m) = 18563 -// vectors for a 12-entry row. Lemma 3.6 supplies the lower bound that seeds and prunes the search. -// -// TODO: extend the fallback path to the row generation of Section 6, which reaches rows the two -// heuristics below leave at a larger magnitude than necessary. The separation oracle in step 3 is a -// 0-1 knapsack, max w.x subject to a.x <= b, which solve_knapsack_problem in cuts/cuts.hpp already -// solves by DP. Boyd (1993) Generating Fenchel Cutting Planes for Knapsack Polyhedra, SIAM J. -// Optim. 3(4):734-750, treats the same separation problem as cut generation. +// Bradley, Hammer and Wolsey (1974), "Coefficient reduction for inequalities in 0-1 variables." +// Theorem 2.5 separates maximal feasible from minimal infeasible points. Positive normalized +// coefficients make activity monotone, so these supersets of BHW's ceilings and roofs suffice. +// Search non-negative, non-increasing weights by increasing max|w|; Lemma 3.6 bounds and prunes it. namespace cuopt::mathematical_optimization::mip { -namespace { - -// One row in BHW's normalized frame: condition N1 (every coefficient positive, negatives -// complemented by x_i = 1 - y_i) and condition N3 (coefficients sorted descending). N1 is what -// makes the row activity monotone in x, and it also makes the weights of any equivalent inequality -// non-negative, since a variable that tightens the row must tighten every inequality with the same -// feasible set. Both normalizations are undone before any reduction is returned. +// N1 complements negative coefficients; N3 sorts them descending. N1 makes activity monotone and +// the weights of equivalent inequalities non-negative. struct norm_row_t { int k = 0; std::array coef{}; // descending, all > 0 @@ -73,7 +48,7 @@ struct partition_t { int lemma36_bound = 0; }; -int64_t weight_activity(const int64_t* w, uint32_t mask) +static int64_t weight_activity(const int64_t* w, uint32_t mask) { int64_t sum = 0; while (mask != 0u) { @@ -84,7 +59,7 @@ int64_t weight_activity(const int64_t* w, uint32_t mask) } template -bool integerization_preserves_binary_feasible_set( +static bool integerization_preserves_binary_feasible_set( const f_t* coefficients, int len, f_t side, @@ -94,7 +69,6 @@ bool integerization_preserves_binary_feasible_set( { cuopt_assert(len >= 2 && len <= BHW_MAX_LEN, "row length outside the enumerable range"); const uint32_t n_pat = 1u << len; - // software implemented, but addition is very cheap std::vector<_Float128> original_activity(n_pat, 0.0L); std::vector integral_activity(n_pat, 0); for (uint32_t m = 1; m < n_pat; ++m) { @@ -113,9 +87,8 @@ bool integerization_preserves_binary_feasible_set( return true; } -// Splits {0,1}^k and collects the maximal feasible and minimal infeasible points. Returns false for -// a degenerate row (all points feasible or all infeasible), which is left to other presolvers. -bool build_partition(const norm_row_t& row, partition_t& out) +// Collects maximal feasible and minimal infeasible points; rejects degenerate partitions. +static bool build_partition(const norm_row_t& row, partition_t& out) { const int k = row.k; const uint32_t n_pat = 1u << k; @@ -168,11 +141,9 @@ bool build_partition(const norm_row_t& row, partition_t& out) return true; } -// BHW Theorem 2.5: (w, t) is equivalent to the row iff every maximal feasible point M satisfies -// sum_{i in M} w_i <= t and every minimal infeasible point satisfies sum_{i in .} w_i >= t + 1. -// Taking t as the feasible-side maximum makes the first block hold by construction, leaving hi < -// lo; integer weights make "> t" and ">= t+1" the same statement. -bool accepts(const partition_t& part, const int64_t* w, int64_t& bound) +// BHW Theorem 2.5: integer w is equivalent iff its maximum over maximal feasible points is below +// its minimum over minimal infeasible points. +static bool accepts(const partition_t& part, const int64_t* w, int64_t& bound) { int64_t hi = std::numeric_limits::min(); for (uint32_t m : part.maximal_feasible) @@ -184,14 +155,9 @@ bool accepts(const partition_t& part, const int64_t* w, int64_t& bound) return hi < lo; } -// The rewritten row must not enlarge this row's LP relaxation: every x in [0,1]^k with w.x <= t has -// to satisfy a.x <= rhs. With a > 0 and w >= 0 after N1 that is a fractional knapsack, max a.x -// subject to w.x <= t, solved by taking the zero-weight entries for free and then filling capacity -// in decreasing a_i/w_i order. Greedy leaves at most one fractional entry, so the comparison closes -// exactly in rationals and no tolerance is needed. The 0-1 knapsack DP in cuts/cuts.hpp cannot be -// substituted here: its optimum is a lower bound on the continuous one, which would let a weakening -// row through. -bool lp_no_weakening(const norm_row_t& row, const int64_t* w, int64_t t) +// Maximizing a.x over [0,1]^k with w.x <= t is a fractional knapsack. Greedy leaves at most one +// fractional entry, allowing the containment check to close exactly in integer arithmetic. +static bool lp_no_weakening(const norm_row_t& row, const int64_t* w, int64_t t) { cuopt_assert(t >= 0, "acceptance implies the origin is feasible, so the bound is non-negative"); @@ -228,9 +194,9 @@ bool lp_no_weakening(const norm_row_t& row, const int64_t* w, int64_t t) (__int128)row.rhs * w[fractional]; } -// Debug companion to accepts(): checks the 0/1 partition over every point rather than the extremal -// ones the search relies on. -[[maybe_unused]] bool verify_equivalent(const norm_row_t& row, const int64_t* w, int64_t t) +[[maybe_unused]] static bool verify_equivalent(const norm_row_t& row, + const int64_t* w, + int64_t t) { const uint32_t n_pat = 1u << row.k; for (uint32_t m = 0; m < n_pat; ++m) { @@ -257,11 +223,8 @@ struct search_state_t { bool found = false; }; -// Enumerates the non-increasing weight vectors reachable from the prefix fixed so far. N1 makes the -// weights non-negative and N3 makes them non-increasing, so a candidate is just a non-increasing -// sequence bounded by w[0]; Lemma 3.6's strict steps both shrink the branching factor and bound how -// low a position may go while still leaving room for the steps beneath it. -void search_positions(search_state_t& st, int pos) +// Enumerates non-negative, non-increasing weights; Lemma 3.6 bounds each suffix. +static void search_positions(search_state_t& st, int pos) { const int k = st.row->k; if (pos == k) { @@ -275,8 +238,7 @@ void search_positions(search_state_t& st, int pos) nonzeros += st.w[i] != 0 ? 1 : 0; sum += st.w[i]; } - // Same max|w| by construction at this depth, so prefer dropping variables, then smaller - // weights. + // At fixed max|w|, prefer fewer nonzeros, then smaller sum. if (st.found && (nonzeros > st.best_nonzeros || (nonzeros == st.best_nonzeros && sum >= st.best_sum))) return; @@ -296,21 +258,16 @@ void search_positions(search_state_t& st, int pos) } } -// Fallback for the rows whose smallest equivalent magnitude exceeds BHW_EXACT_MAX_WEIGHT, where the -// exhaustive search gives up: w = round(a / min a) and the all-ones clause form, both put through -// the same acceptance and LP-strength gates as a searched vector. N3 already sorted a, so both -// candidates are non-increasing. -bool heuristic_reduce(const norm_row_t& row, - const partition_t& part, - std::vector& weights, - int64_t& bound) +// Above the exact-search cap, try round(a/min(a)) and the all-ones row through the same gates. +static bool heuristic_reduce(const norm_row_t& row, + const partition_t& part, + std::vector& weights, + int64_t& bound) { const int k = row.k; const int64_t a_min = row.coef[k - 1]; cuopt_assert(a_min > 0, "N1 leaves every coefficient positive"); - // Only a strict gain is worth installing: smaller magnitude, or the same magnitude with a - // variable dropped. int64_t best_max = row.coef[0]; int best_nonzeros = k; bool found = false; @@ -342,13 +299,12 @@ bool heuristic_reduce(const norm_row_t& row, return found; } -// Reduce one normalized shape. The caller undoes N1/N3 on the result. -bool reduce_shape(const norm_row_t& row, std::vector& weights, int64_t& bound) +static bool reduce_shape(const norm_row_t& row, std::vector& weights, int64_t& bound) { partition_t part; if (!build_partition(row, part)) return false; - const int64_t current = row.coef[0]; // N3 puts the largest coefficient first + const int64_t current = row.coef[0]; cuopt_assert(current >= 2, "rows already at magnitude one are rejected before normalization"); // Lemma 3.6 bounds max|w| from below over every equivalent inequality, so this row is provably // irreducible in magnitude and not worth searching. @@ -371,8 +327,6 @@ bool reduce_shape(const norm_row_t& row, std::vector& weights, int64_t& return heuristic_reduce(row, part, weights, bound); } -} // namespace - template bhw_row_rewrite_t bhw_reduce_row( const f_t* coefficients, int len, f_t side, int direction, bhw_shape_cache_t* cache) @@ -383,7 +337,6 @@ bhw_row_rewrite_t bhw_reduce_row( if (len < 2 || len > BHW_MAX_LEN) return rewrite; if (!scaling_bound_finite(side)) return rewrite; - // Integerize so the point partition is exact, then orient the row to a.x <= b. const double scale = row_int_scale( coefficients, len, side, std::numeric_limits::infinity(), BHW_MAX_LEN, BHW_INT_SCALE_MAX); if (scale == 0.0) return rewrite; @@ -395,8 +348,7 @@ bhw_row_rewrite_t bhw_reduce_row( if (integral[j] == 0) return rewrite; largest = std::max(largest, std::abs(integral[j])); } - // A row already at +/-1 has no magnitude to give back; rejecting it here is what keeps screening - // cheap on the rows that dominate a model. + // can't coefficient-reduce a unit-magnitude row any further if (largest <= 1) return rewrite; norm_row_t norm_row; @@ -440,7 +392,7 @@ bhw_row_rewrite_t bhw_reduce_row( } if (!result->accepted) return rewrite; - // Rows sharing an integerized cache key can have different original floating-point partitions. + // check the feasible set remains unchanged under floating point math if (!integerization_preserves_binary_feasible_set( coefficients, len, side, direction, integral, integral_side)) return rewrite; @@ -479,10 +431,6 @@ bhw_row_rewrite_t bhw_reduce_row( return rewrite; } -namespace { - -// Coefficient-shrink figures behind the DEBUG line. Both the accumulation and the summary compile -// away below DEBUG, so a release build carries neither the per-row vector nor the selection. struct bhw_stats_t { #if (CUOPT_LOG_ACTIVE_LEVEL <= RAPIDS_LOGGER_LOG_LEVEL_DEBUG) int64_t coefficients_reduced = 0; @@ -507,7 +455,6 @@ struct bhw_stats_t { cuopt_assert(n_rows_rewritten > 0, "a changed coefficient implies an accepted row"); const double mean = std::accumulate(row_shrinks.begin(), row_shrinks.end(), 0.0) / n_rows_rewritten; - // One row collapsing to magnitude 1 dominates the mean, so report the median next to it. const auto middle = row_shrinks.begin() + n_rows_rewritten / 2; std::nth_element(row_shrinks.begin(), middle, row_shrinks.end()); double median = *middle; @@ -530,8 +477,6 @@ struct bhw_stats_t { #endif }; -} // namespace - template papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& problem, const papilo::ProblemUpdate& problemUpdate, @@ -554,13 +499,9 @@ papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& papilo::PresolveStatus status = papilo::PresolveStatus::kUnchanged; bhw_stats_t stats; - // Every eligible row is screened, not only problemUpdate.getChangedActivities(): that worklist is - // seeded in full once and afterwards fed only by activity changes, so a row whose side changes is - // never revisited. Rescreening is cheap because the shape cache absorbs the repetition and the - // reduction is idempotent. + // getChangedActivities() omits side-only changes, so screen every row. cache hits amortize for (int row = 0; row < num_rows; ++row) { if (reductions.size() >= presolve_options.max_reduction_seq) break; - // Screening one row is cheap next to the interrupt check. if (papilo::PresolveMethod::is_interrupted( timer, presolve_options.tlim, presolve_options.early_exit_callback)) break; @@ -599,11 +540,6 @@ papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& "an accepted rewrite keeps at least one nonzero weight"); stats.rewrote_row(rewrite.max_coef_before, rewrite.max_coef_after); - // Same shape as papilo's own sparsifier, SimplifyInequalities: lock the row, then the entries, - // then the side. Dropping an entry needs no column lock: ProblemUpdate marks the column - // modified and derives the resulting singleton rows and empty columns itself. The zero has to - // be exact: SparseStorage::changeRowInplace compacts an entry out on newval == 0, not on - // num.isZero. papilo::TransactionGuard guard{reductions}; reductions.lockRow(row); [[maybe_unused]] int emitted = 0; @@ -624,8 +560,6 @@ papilo::PresolveStatus BHWCoeffReduce::execute(const papilo::Problem& ++emitted; } } - // Reporting kReduced for a transaction that changes nothing would have papilo re-derive the - // same rewrite every round. cuopt_assert(emitted > 0, "accepted rewrite emitted no reduction"); status = papilo::PresolveStatus::kReduced; } diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp index fee1d27334..1070498ccd 100644 --- a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp @@ -26,56 +26,36 @@ namespace cuopt::mathematical_optimization::mip { -// Widest row we enumerate: building the point partition walks 2^BHW_MAX_LEN patterns, so this caps -// the cost of screening one row. Wider rows are left alone; nearly every row that integerizes -// exactly is narrower than this. +// Building the point partition visits at most 2^BHW_MAX_LEN patterns per row. static constexpr int BHW_MAX_LEN = 12; // Largest max|w| the exhaustive search considers before falling back to the heuristic candidates. -// Reductions needing a larger magnitude are rare enough not to pay for the extra search depth. static constexpr int64_t BHW_EXACT_MAX_WEIGHT = 6; -// Largest per-row rational multiplier / denominator used to integerize a row (passed to -// row_int_scale as its maxdnom/maxfinal caps). -static constexpr int64_t BHW_INT_SCALE_MAX = 1000000; // 1e6 +// Passed to row_int_scale as its maxdnom and maxfinal caps. +static constexpr int64_t BHW_INT_SCALE_MAX = 1000000; -// Outcome for one canonical row shape, in BHW's normalized frame (coefficients complemented to be -// positive and sorted descending). Rejections are cached too, since re-deriving them is the bulk of -// the work on instances whose rows repeat. struct bhw_shape_result_t { std::vector weights; int64_t bound = 0; bool accepted = false; }; -// Keyed by the normalized coefficients followed by the normalized right-hand side. The reduction is -// a pure function of that key, so entries stay valid across presolve rounds and problems. using bhw_shape_cache_t = std::map, bhw_shape_result_t>; struct bhw_row_rewrite_t { - std::vector coefficients; // one per input entry, in input order; 0 drops that entry - int64_t side = 0; // replaces the row's finite side - // Largest |coefficient| before and after reduction, both in the integerized frame. Comparable - // only there: the row is scaled on the way in, so the input coefficients sit in a different - // frame. + std::vector coefficients; // 0 drops the entry + int64_t side = 0; int64_t max_coef_before = 0; int64_t max_coef_after = 0; bool accepted = false; }; -// Rewrite one one-sided all-binary row with smaller integer coefficients spanning the same 0/1 -// feasible set. direction is +1 for "coefficients . x <= side" and -1 for ">= side"; side is the -// finite side of the row. Rejects the row (accepted = false) unless it integerizes exactly to -// nonzero coefficients, admits a strictly smaller equivalent form, and that form does not enlarge -// the row's LP relaxation. cache may be null to skip memoization. -// -// The caller checks that every entry is a binary integer variable and that exactly one side of the -// row is finite. Exposed for testing: BHWCoeffReduce::execute only screens rows and emits the -// result, so this covers the whole reduction without any papilo types. +// Rewrites a one-sided all-binary row with smaller integer coefficients and the same 0/1 feasible +// set. direction is +1 for <= and -1 for >=. The caller guarantees nonfixed +// binary variables and exactly one finite side. template bhw_row_rewrite_t bhw_reduce_row( const f_t* coefficients, int len, f_t side, int direction, bhw_shape_cache_t* cache); -// Bradley-Hammer-Wolsey coefficient reduction: replace an all-binary row by an equivalent one with -// smaller integer coefficients. See bhw_coeff_reduce.cpp for the lineage. template class BHWCoeffReduce : public papilo::PresolveMethod { public: @@ -94,7 +74,6 @@ class BHWCoeffReduce : public papilo::PresolveMethod { int& reason_of_infeasibility) override; private: - // Only touched from execute, which papilo runs one task at a time per presolver object. bhw_shape_cache_t shape_cache_; }; diff --git a/cpp/tests/mip/bhw_coeff_reduce_test.cpp b/cpp/tests/mip/bhw_coeff_reduce_test.cpp index 2bd7dbf4bf..9d3801fd13 100644 --- a/cpp/tests/mip/bhw_coeff_reduce_test.cpp +++ b/cpp/tests/mip/bhw_coeff_reduce_test.cpp @@ -32,8 +32,7 @@ bhw_row_rewrite_t reduce(const std::vector& coefficients, coefficients.data(), (int)coefficients.size(), side, direction, cache); } -// The whole point of the pass: the rewritten row must accept exactly the same 0/1 points as the -// original. Checked over every point, independently of the extremal-point test the search uses. +// Full 0/1 check independent of the production extremal-point test. bool same_feasible_set(const std::vector& coefficients, double side, int direction, @@ -59,17 +58,12 @@ bool same_feasible_set(const std::vector& coefficients, } // namespace -// Bradley, Hammer and Wolsey (1974) open with this row and reduce it to 4,4,2,2,1,1,1,0 <= 5. That -// rewrite enlarges the LP relaxation (it admits a fractional point of activity 91.75 against a -// right-hand side of 80), so the LP-strength check rejects it and no smaller equivalent form -// survives. +// BHW's opening example has an equivalent reduction that weakens the LP relaxation. TEST(bhw_coeff_reduce, rejects_a_rewrite_that_weakens_the_relaxation) { EXPECT_FALSE(reduce({65, 64, 41, 22, 13, 12, 8, 2}, 80).accepted); } -// The shape that motivated the pass: a 9/10 coefficient against thirds, with a fractional -// right-hand side. Integerizing and reducing lands it in int8 range. TEST(bhw_coeff_reduce, rational_row_integerizes_and_reduces) { const std::vector row{0.9, 1.0 / 3, 1.0 / 3, 1.0 / 3}; @@ -91,7 +85,6 @@ TEST(bhw_coeff_reduce, rejects_approximate_integerization_that_changes_the_feasi EXPECT_FALSE(reduce(row, 11000.0 - perturbation).accepted); } -// The >= orientation is normalized by negation, so the same row negated must come back negated. TEST(bhw_coeff_reduce, greater_equal_row_keeps_its_orientation) { const std::vector row{-0.9, -1.0 / 3, -1.0 / 3, -1.0 / 3}; @@ -104,16 +97,15 @@ TEST(bhw_coeff_reduce, greater_equal_row_keeps_its_orientation) TEST(bhw_coeff_reduce, rejects_rows_with_nothing_to_give_back) { - // Already at unit magnitude. + // Unit magnitude. EXPECT_FALSE(reduce({1, 1, 1, 1}, 2).accepted); - // Does not integerize within the rational cap. + // Scaling cap. EXPECT_FALSE(reduce({M_PI, 1, 1}, 2).accepted); - // Outside the enumerable width. + // Unsupported width. EXPECT_FALSE(reduce({5}, 2).accepted); EXPECT_FALSE(reduce(std::vector(BHW_MAX_LEN + 1, 3.0), 5).accepted); - // Every point feasible, so there is nothing to separate. + // Degenerate partition. EXPECT_FALSE(reduce({3, 2, 2}, 100).accepted); - // No point feasible. EXPECT_FALSE(reduce({3, 2, 2}, -1).accepted); } @@ -140,8 +132,6 @@ TEST(bhw_coeff_reduce, memoized_result_matches_the_uncached_one) } } -// The invariant that matters, over mixed signs, both orientations and rational coefficients: an -// accepted rewrite never changes which 0/1 points satisfy the row. TEST(bhw_coeff_reduce, accepted_rewrites_preserve_the_feasible_set) { std::mt19937_64 rng(20260805); From 9043c70a77191a2e8799130669cfa380448e7bf0 Mon Sep 17 00:00:00 2001 From: yboucher Date: Fri, 28 Aug 2026 01:01:11 -0700 Subject: [PATCH 6/8] style --- cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp | 12 +++++------- cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp | 2 +- cpp/tests/mip/bhw_coeff_reduce_test.cpp | 4 ++-- 3 files changed, 8 insertions(+), 10 deletions(-) diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp index 55335a2d97..598918f67e 100644 --- a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp @@ -79,8 +79,8 @@ static bool integerization_preserves_binary_feasible_set( } for (uint32_t m = 0; m < n_pat; ++m) { - const bool original_feasible = direction == 1 ? original_activity[m] <= side - : original_activity[m] >= side; + const bool original_feasible = + direction == 1 ? original_activity[m] <= side : original_activity[m] >= side; const bool integral_feasible = integral_activity[m] <= integral_side; if (original_feasible != integral_feasible) return false; } @@ -194,9 +194,7 @@ static bool lp_no_weakening(const norm_row_t& row, const int64_t* w, int64_t t) (__int128)row.rhs * w[fractional]; } -[[maybe_unused]] static bool verify_equivalent(const norm_row_t& row, - const int64_t* w, - int64_t t) +[[maybe_unused]] static bool verify_equivalent(const norm_row_t& row, const int64_t* w, int64_t t) { const uint32_t n_pat = 1u << row.k; for (uint32_t m = 0; m < n_pat; ++m) { @@ -352,9 +350,9 @@ bhw_row_rewrite_t bhw_reduce_row( if (largest <= 1) return rewrite; norm_row_t norm_row; - norm_row.k = len; + norm_row.k = len; const int64_t integral_side = std::llround((double)side * scale) * direction; - norm_row.rhs = integral_side; + norm_row.rhs = integral_side; std::array order{}; for (int j = 0; j < len; ++j) { order[j] = j; diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp index 1070498ccd..73cc5d833e 100644 --- a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp @@ -43,7 +43,7 @@ using bhw_shape_cache_t = std::map, bhw_shape_result_t>; struct bhw_row_rewrite_t { std::vector coefficients; // 0 drops the entry - int64_t side = 0; + int64_t side = 0; int64_t max_coef_before = 0; int64_t max_coef_after = 0; bool accepted = false; diff --git a/cpp/tests/mip/bhw_coeff_reduce_test.cpp b/cpp/tests/mip/bhw_coeff_reduce_test.cpp index 9d3801fd13..285b6d8d79 100644 --- a/cpp/tests/mip/bhw_coeff_reduce_test.cpp +++ b/cpp/tests/mip/bhw_coeff_reduce_test.cpp @@ -76,8 +76,8 @@ TEST(bhw_coeff_reduce, rational_row_integerizes_and_reduces) TEST(bhw_coeff_reduce, rejects_approximate_integerization_that_changes_the_feasible_set) { - constexpr int num_variables = 12; - constexpr double perturbation = 9.9e-7; + constexpr int num_variables = 12; + constexpr double perturbation = 9.9e-7; constexpr double small_coefficient = 1000.0; std::vector row(num_variables, small_coefficient + perturbation); row[0] = 12000.0; From 2c476e9c954f98ac0e84abcfe3814e421b987ac7 Mon Sep 17 00:00:00 2001 From: yboucher Date: Fri, 28 Aug 2026 09:47:19 -0700 Subject: [PATCH 7/8] bhw runs as delayed --- cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp index 73cc5d833e..91b7f5feb2 100644 --- a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp @@ -64,6 +64,8 @@ class BHWCoeffReduce : public papilo::PresolveMethod { this->setName("bhwcoeffreduce"); this->setType(papilo::PresolverType::kIntegralCols); this->setTiming(papilo::PresolverTiming::kMedium); + // can interfere with some papilo reductions by causing them to miss their trigger condition + this->setDelayed(true); } papilo::PresolveStatus execute(const papilo::Problem& problem, From 7e1de69d9c7089ffcf791363e8616fce63c16dc2 Mon Sep 17 00:00:00 2001 From: yboucher Date: Fri, 28 Aug 2026 10:12:05 -0700 Subject: [PATCH 8/8] style --- cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp index 91b7f5feb2..d334bf399f 100644 --- a/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp @@ -65,7 +65,7 @@ class BHWCoeffReduce : public papilo::PresolveMethod { this->setType(papilo::PresolverType::kIntegralCols); this->setTiming(papilo::PresolverTiming::kMedium); // can interfere with some papilo reductions by causing them to miss their trigger condition - this->setDelayed(true); + this->setDelayed(true); } papilo::PresolveStatus execute(const papilo::Problem& problem,