diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index 5fa939058c..0776e5972e 100644 --- a/cpp/src/mip_heuristics/CMakeLists.txt +++ b/cpp/src/mip_heuristics/CMakeLists.txt @@ -14,6 +14,7 @@ set(MIP_LP_NECESSARY_FILES ${CMAKE_CURRENT_SOURCE_DIR}/presolve/third_party_presolve.cpp ${CMAKE_CURRENT_SOURCE_DIR}/presolve/single_lock_dual_aggregation.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..598918f67e --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.cpp @@ -0,0 +1,585 @@ +/* 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 + +// 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 { + +// 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 + 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; +}; + +static 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; +} + +template +static 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; + 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; +} + +// 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; + 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: 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) + 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; +} + +// 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"); + + 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]; +} + +[[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) { + 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 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) { + 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]; + } + // 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; + 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); + } +} + +// 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"); + + 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; +} + +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]; + 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); +} + +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; + + 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] = std::llround((double)coefficients[j] * scale) * direction; + if (integral[j] == 0) return rewrite; + largest = std::max(largest, std::abs(integral[j])); + } + // can't coefficient-reduce a unit-magnitude row any further + if (largest <= 1) return rewrite; + + norm_row_t norm_row; + 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; + // 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; + + // 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; + + 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(), 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; +} + +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; + 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 +}; + +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 auto& presolve_options = problemUpdate.getPresolveOptions(); + + const int num_rows = constraint_matrix.getNRows(); + papilo::PresolveStatus status = papilo::PresolveStatus::kUnchanged; + bhw_stats_t stats; + + // 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; + if (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; + + 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"); + stats.rewrote_row(rewrite.max_coef_before, rewrite.max_coef_after); + + 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; + stats.changed_coefficient(rewrite.coefficients[j]); + } + 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; + } + } + cuopt_assert(emitted > 0, "accepted rewrite emitted no reduction"); + status = papilo::PresolveStatus::kReduced; + } + + stats.report(); + + 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..d334bf399f --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/bhw_coeff_reduce.hpp @@ -0,0 +1,82 @@ +/* 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 { + +// 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. +static constexpr int64_t BHW_EXACT_MAX_WEIGHT = 6; +// Passed to row_int_scale as its maxdnom and maxfinal caps. +static constexpr int64_t BHW_INT_SCALE_MAX = 1000000; + +struct bhw_shape_result_t { + std::vector weights; + int64_t bound = 0; + bool accepted = false; +}; + +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 max_coef_before = 0; + int64_t max_coef_after = 0; + bool accepted = false; +}; + +// 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); + +template +class BHWCoeffReduce : public papilo::PresolveMethod { + public: + BHWCoeffReduce() : 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, + const papilo::ProblemUpdate& problemUpdate, + const papilo::Num& num, + papilo::Reductions& reductions, + const papilo::Timer& timer, + int& reason_of_infeasibility) override; + + private: + 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 b6976e430b..a8f32d3e62 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 @@ -678,6 +679,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 be78695698..b69fa08f46 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/single_lock_dual_aggregation_test.cpp ${CUOPT_TEST_DIR}/mip/termination_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..285b6d8d79 --- /dev/null +++ b/cpp/tests/mip/bhw_coeff_reduce_test.cpp @@ -0,0 +1,179 @@ +/* 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); +} + +// Full 0/1 check independent of the production extremal-point test. +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 + +// 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); +} + +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)); +} + +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); +} + +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) +{ + // Unit magnitude. + EXPECT_FALSE(reduce({1, 1, 1, 1}, 2).accepted); + // Scaling cap. + EXPECT_FALSE(reduce({M_PI, 1, 1}, 2).accepted); + // Unsupported width. + EXPECT_FALSE(reduce({5}, 2).accepted); + EXPECT_FALSE(reduce(std::vector(BHW_MAX_LEN + 1, 3.0), 5).accepted); + // Degenerate partition. + EXPECT_FALSE(reduce({3, 2, 2}, 100).accepted); + 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; + 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); + } + } +} + +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