diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 01afb9f98f..8fc0784899 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -158,6 +158,10 @@ /* @brief Block bounded-variable-elimination step of cuOpt's internal MIP presolve */ #define CUOPT_MIP_HYPER_BLOCK_BVE "mip_hyper_block_bve" +/* @brief Indicator-strengthening step that runs before Papilo presolve on MIPs */ +#define CUOPT_MIP_HYPER_PRESOLVE_INDICATOR_STRENGTHENING \ + "mip_hyper_presolve_indicator_strengthening" + /* @brief QCQP (barrier) scaling hyper-parameters */ #define CUOPT_QCQP_HYPER_RUIZ_EQUILIBRATION "qcqp_hyper_ruiz_equilibration" diff --git a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp index f3dc5fe340..07204846e6 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp @@ -176,6 +176,13 @@ class mip_solver_settings_t { * no-op when no certified reduction exists. */ bool block_bve{true}; + /** + * @brief Enable the indicator-strengthening step of presolve (MIP only). + * + * Runs before Papilo and only when the higher-level presolve is enabled. It appends implied + * indicator rows and lifts capacity rows by the indicator that bounds all of their members. + */ + bool indicator_strengthening{true}; /** * @brief Determinism mode for MIP solver. * diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 49a3e4fd2b..dff8f4c6cb 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -264,6 +264,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings // Recursive sub-MIP (RINS) hyper-parameters (hidden from default --help: name contains "hyper_") {CUOPT_MIP_HYPER_SUBMIP_ENABLE_CPUFJ, &mip_settings.submip_params.enable_cpufj, true, "run CPU FJ over the sub-MIP"}, {CUOPT_MIP_HYPER_BLOCK_BVE, &mip_settings.block_bve, true, "eliminate blocks of binaries in cuOpt's MIP presolve (needs " CUOPT_MIP_PROBING ")"}, + {CUOPT_MIP_HYPER_PRESOLVE_INDICATOR_STRENGTHENING, &mip_settings.indicator_strengthening, true, "append implied indicator rows and lift capacity rows before Papilo presolve"}, // PDLP scaling hyper-parameter (hidden from default --help: name contains "hyper_") {CUOPT_PDLP_HYPER_ENABLE_CURTIS_REID_SCALING, &pdlp_settings.hyper_params.do_curtis_reid_scaling, true, "Curtis-Reid prescaling, run before Ruiz/Pock-Chambolle scaling"}, }; diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index 6974fd7d35..db77c595d8 100644 --- a/cpp/src/mip_heuristics/CMakeLists.txt +++ b/cpp/src/mip_heuristics/CMakeLists.txt @@ -15,6 +15,7 @@ set(MIP_LP_NECESSARY_FILES ${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}/presolve/indicator_strengthening.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/indicator_strengthening.cpp b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp new file mode 100644 index 0000000000..d3fa383ec4 --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp @@ -0,0 +1,425 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "indicator_strengthening.hpp" + +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace cuopt::mathematical_optimization::mip { + +namespace { + +// Indicator variables $z$ are binaries that must be paid for before anything they own may be used, +// while member variables $x_j$ are the variables owned by a given indicator. +// They are linked via the following constraints: +// Link -> x_j - z <= 0 +// Disjunction -> y - \sum_{j \in \mathcal{S}} x_j <= 0 +// Capacity -> \sum_{i \in \mathcal{S}} x_i - s <= K +template +class indicator_strengthening_t { + public: + indicator_strengthening_t(const papilo::Problem& problem); + i_t add_implied_indicator_rows(papilo::Vec>& implied_entries) const; + i_t find_lift_capacity_rows(papilo::Vec>& lifted_entries) const; + i_t num_variable_upper_bounds() const { return variable_upper_bound_indicators.size(); } + i_t row_orientation(i_t row) const { return orientation[row]; } + + private: + const papilo::Problem& problem; + + std::vector is_binary; + + // +1 when the stored row reads A[i, :]^T x <= rhs, -1 when it reads A[i, :]^T x >= lhs, 0 for + // equations, ranges and free rows. + std::vector orientation; + + // CSR storing the implication graph for binary variables of the form x <= z, where x is a + // member variable and z, an indicator. Each row of the CSR corresponds to a member variable + // x and each column the indicator z that owns x. Here the coefficients of the CSR are irrelevant, + // so we just store the offsets (`variable_upper_bound_offsets`) and the nonzero columns + // (`variable_upper_bound_indicators`). + std::vector variable_upper_bound_offsets; + std::vector variable_upper_bound_indicators; + + i_t max_row_size = 0; +}; + +template +indicator_strengthening_t::indicator_strengthening_t(const papilo::Problem& problem) + : problem(problem) +{ + 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& row_sizes = constraint_matrix.getRowSizes(); + 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 i_t num_rows = constraint_matrix.getNRows(); + const i_t num_cols = problem.getNCols(); + + // A binary is an integral column with finite bounds [0, 1]. + is_binary.resize(num_cols); + for (i_t col = 0; col < num_cols; ++col) { + is_binary[col] = col_flags[col].test(papilo::ColFlag::kIntegral) && + !col_flags[col].test(papilo::ColFlag::kLbInf, papilo::ColFlag::kUbInf) && + lower_bounds[col] == 0.0 && upper_bounds[col] == 1.0; + } + + // Record the orientation of every one-sided row, and collect each two-variable row x - z <= 0 + // between binaries as a variable upper bound (member x, indicator z). + orientation.assign(num_rows, 0); + std::vector> variable_upper_bounds; + variable_upper_bounds.reserve(std::count(row_sizes.begin(), row_sizes.end(), 2)); + for (i_t row = 0; row < num_rows; ++row) { + max_row_size = std::max(max_row_size, row_sizes[row]); + const bool lhs_infinite = row_flags[row].test(papilo::RowFlag::kLhsInf); + const bool rhs_infinite = row_flags[row].test(papilo::RowFlag::kRhsInf); + if (lhs_infinite == rhs_infinite) { continue; } + const i_t direction = lhs_infinite ? 1 : -1; + orientation[row] = direction; + if (row_sizes[row] != 2) { continue; } + const f_t side = direction == 1 ? rhs_values[row] : lhs_values[row]; + if (side != 0.0) { continue; } + + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + const i_t* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); + if (!is_binary[indices[0]] || !is_binary[indices[1]]) { continue; } + const f_t v0 = direction * values[0]; + const f_t v1 = direction * values[1]; + if (v0 == 1.0 && v1 == -1.0) { + variable_upper_bounds.emplace_back(indices[0], indices[1]); + } else if (v0 == -1.0 && v1 == 1.0) { + variable_upper_bounds.emplace_back(indices[1], indices[0]); + } + } + + // Build the CSR from the (member, indicator) pairs: count the indicators of each member, + // prefix-sum the counts into offsets, then scatter the indicators into place. + variable_upper_bound_offsets.assign(num_cols + 1, 0); + for (const auto& variable_upper_bound : variable_upper_bounds) { + ++variable_upper_bound_offsets[variable_upper_bound.first + 1]; + } + for (i_t col = 0; col < num_cols; ++col) { + variable_upper_bound_offsets[col + 1] += variable_upper_bound_offsets[col]; + } + variable_upper_bound_indicators.resize(variable_upper_bounds.size()); + std::vector next(variable_upper_bound_offsets.begin(), + variable_upper_bound_offsets.end() - 1); + for (const auto& [member, indicator] : variable_upper_bounds) { + variable_upper_bound_indicators[next[member]++] = indicator; + } +} + +template +i_t indicator_strengthening_t::add_implied_indicator_rows( + papilo::Vec>& implied_entries) const +{ + const auto& constraint_matrix = problem.getConstraintMatrix(); + const auto& lhs_values = constraint_matrix.getLeftHandSides(); + const auto& rhs_values = constraint_matrix.getRightHandSides(); + const auto& row_sizes = constraint_matrix.getRowSizes(); + + const i_t num_rows = constraint_matrix.getNRows(); + const i_t num_cols = problem.getNCols(); + + i_t num_implied = 0; + + std::vector indicators; + indicators.reserve(max_row_size); + std::vector mark(num_cols, -1); + + // Loop over the constraints and identify the disjunctions (y - \sum_{j \in \mathcal{S}} x_j <= + // 0), + // $|\mathcal{S}| >= 2$ then construct a new constraint (y - \sum_{i \in \mathcal{D}} z_i <= 0), + // where + // $\mathcal{D}$ contains the first indicator associated with a member variable x_j. + for (i_t row = 0; row < num_rows; ++row) { + const i_t direction = orientation[row]; + if (direction == 0 || row_sizes[row] < 3) { continue; } + const f_t side = direction == 1 ? rhs_values[row] : lhs_values[row]; + if (side != 0.0) { continue; } + + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + const i_t len = row_coefficients.getLength(); + const i_t* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); + + // Every entry must be binary: the single +1 is the head y, and each -1 member is replaced by + // its indicator, deduplicated with `mark`. The row is rejected if no two members share an + // indicator, since the new row would then be implied by the LP relaxation. + const size_t num_members = len - 1; + i_t head = -1; + bool usable = true; + indicators.clear(); + for (i_t p = 0; p < len && usable; ++p) { + const i_t col = indices[p]; + const f_t v = direction * values[p]; + if (!is_binary[col]) { + usable = false; + } else if (v == 1.0) { + usable = head < 0; + head = col; + } else if (v == -1.0) { + const i_t start = variable_upper_bound_offsets[col]; + const i_t end = variable_upper_bound_offsets[col + 1]; + + // Take the first indicator that owns this variable, or the member itself if it has none. + const i_t z = start < end ? variable_upper_bound_indicators[start] : col; + if (mark[z] != row) { + mark[z] = row; + indicators.push_back(z); + usable = indicators.size() < num_members; + } + } else { + usable = false; + } + } + // Skip rows without a head, or whose head is also one of the indicators (a trivial row). + if (!usable || head < 0 || mark[head] == row) { continue; } + + // In the previous loop we insert the indicators in a random order, however, Papilo + // expects the triplets (row, col, val) to be sorted by row and then by column. + std::sort(indicators.begin(), indicators.end()); + const i_t implied_row = num_rows + num_implied; + const auto split = std::lower_bound(indicators.begin(), indicators.end(), head); + for (auto it = indicators.begin(); it != split; ++it) { + implied_entries.emplace_back(implied_row, *it, f_t{-1}); + } + implied_entries.emplace_back(implied_row, head, f_t{1}); + for (auto it = split; it != indicators.end(); ++it) { + implied_entries.emplace_back(implied_row, *it, f_t{-1}); + } + ++num_implied; + } + + return num_implied; +} + +template +i_t indicator_strengthening_t::find_lift_capacity_rows( + papilo::Vec>& lifted_entries) const +{ + const auto& constraint_matrix = problem.getConstraintMatrix(); + const auto& lhs_values = constraint_matrix.getLeftHandSides(); + const auto& rhs_values = constraint_matrix.getRightHandSides(); + const auto& domains = problem.getVariableDomains(); + const auto& col_flags = domains.flags; + const auto& lower_bounds = domains.lower_bounds; + + const i_t num_rows = constraint_matrix.getNRows(); + + i_t num_lifted = 0; + + std::vector members; + members.reserve(max_row_size); + + // Loop over the constraints and identify the capacity row \sum_{i \in \mathcal{S}} x_i - s <= K, + // and then tighten it as \sum_{i \in \mathcal{S}} x_i - s <= Kz, for a indicator z with + // x_i <= z + for (i_t row = 0; row < num_rows; ++row) { + // A capacity row is a one-sided row with a positive right-hand side K. + const i_t direction = orientation[row]; + if (direction == 0) { continue; } + const f_t capacity = direction * (direction == 1 ? rhs_values[row] : lhs_values[row]); + if (capacity <= 0.0) { continue; } + + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + const i_t len = row_coefficients.getLength(); + const i_t* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); + + // Every integral entry must be a +1 binary member with at least one indicator, and every + // continuous entry a -s slack with s >= 0. The pivot is the member with the fewest indicators. + members.clear(); + i_t pivot = -1; + bool usable = true; + for (i_t p = 0; p < len && usable; ++p) { + const i_t col = indices[p]; + const f_t v = direction * values[p]; + if (col_flags[col].test(papilo::ColFlag::kIntegral)) { + const i_t num_indicators = + variable_upper_bound_offsets[col + 1] - variable_upper_bound_offsets[col]; + usable = is_binary[col] && v == 1.0 && num_indicators > 0; + if (!usable) { continue; } + members.push_back(col); + + const i_t num_indicators_pivot = + pivot >= 0 ? variable_upper_bound_offsets[pivot + 1] - variable_upper_bound_offsets[pivot] + : 0; + if (pivot < 0 || num_indicators < num_indicators_pivot) { pivot = col; } + } else { + usable = + v < 0.0 && !col_flags[col].test(papilo::ColFlag::kLbInf) && lower_bounds[col] >= 0.0; + } + } + // Lifting needs at least two members, and only tightens the row when K < |S|. + if (!usable || members.size() < 2) { continue; } + const f_t n_members = members.size(); + if (capacity >= n_members) { continue; } + + // A shared indicator must be in the pivot's list, so only its candidates are tried; each is + // searched for in the lists of the other members, and the first one found in all of them wins. + i_t indicator = -1; + i_t pivot_start = variable_upper_bound_offsets[pivot]; + i_t pivot_end = variable_upper_bound_offsets[pivot + 1]; + + for (i_t p = pivot_start; p < pivot_end && indicator < 0; ++p) { + const i_t z = variable_upper_bound_indicators[p]; + bool shared = true; + for (size_t k = 0; k < members.size() && shared; ++k) { + if (members[k] == pivot) { continue; } + i_t start = variable_upper_bound_offsets[members[k]]; + i_t end = variable_upper_bound_offsets[members[k] + 1]; + std::span indicator_span(variable_upper_bound_indicators.begin() + start, end - start); + shared = std::find(indicator_span.begin(), indicator_span.end(), z) != indicator_span.end(); + } + if (shared) { indicator = z; } + } + if (indicator < 0) { continue; } + + // Add the -K z entry to the row; `strengthen_indicators` then sets its right-hand side to 0. + lifted_entries.emplace_back(row, indicator, direction * -capacity); + ++num_lifted; + } + + return num_lifted; +} + +} // namespace + +template +void strengthen_indicators(papilo::Problem& problem) +{ + const indicator_strengthening_t strengthening(problem); + if (strengthening.num_variable_upper_bounds() == 0) { return; } + + 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& row_sizes = constraint_matrix.getRowSizes(); + + const i_t num_rows = constraint_matrix.getNRows(); + const i_t num_cols = problem.getNCols(); + + // Bound the number of implied entries and lifted rows from above to reserve the buffers. + i_t max_implied_entries = 0; + i_t max_lifted_rows = 0; + for (i_t row = 0; row < num_rows; ++row) { + const i_t direction = strengthening.row_orientation(row); + if (direction == 0) { continue; } + const f_t capacity = direction * (direction == 1 ? rhs_values[row] : lhs_values[row]); + if (capacity == 0.0 && row_sizes[row] >= 3) { + max_implied_entries += row_sizes[row] - 1; + } else if (capacity > 0.0) { + ++max_lifted_rows; + } + } + + papilo::Vec> implied_entries; + papilo::Vec> lifted_entries; + implied_entries.reserve(max_implied_entries); + lifted_entries.reserve(max_lifted_rows); + const i_t num_implied = strengthening.add_implied_indicator_rows(implied_entries); + const i_t num_lifted = strengthening.find_lift_capacity_rows(lifted_entries); + if (num_implied == 0 && num_lifted == 0) { return; } + + // Rebuild the matrix as row-sorted triplets, inserting each lifted entry at its column position + // within its row, then append the implied rows after the original ones. + papilo::Vec> entries; + entries.reserve(constraint_matrix.getNnz() + lifted_entries.size() + implied_entries.size()); + auto lifted = lifted_entries.begin(); + for (i_t row = 0; row < num_rows; ++row) { + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + const i_t len = row_coefficients.getLength(); + const i_t* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); + i_t p = 0; + if (lifted != lifted_entries.end() && std::get<0>(*lifted) == row) { + for (; p < len && indices[p] < std::get<1>(*lifted); ++p) { + entries.emplace_back(row, indices[p], values[p]); + } + entries.push_back(*lifted); + ++lifted; + } + for (; p < len; ++p) { + entries.emplace_back(row, indices[p], values[p]); + } + } + entries.insert(entries.end(), implied_entries.begin(), implied_entries.end()); + + // Lifted rows get a zero right-hand side, and the implied rows read y - \sum_{z} z <= 0. + papilo::Vec lhs; + papilo::Vec rhs; + papilo::Vec flags; + lhs.reserve(num_rows + num_implied); + rhs.reserve(num_rows + num_implied); + flags.reserve(num_rows + num_implied); + lhs.assign(lhs_values.begin(), lhs_values.end()); + rhs.assign(rhs_values.begin(), rhs_values.end()); + flags.assign(row_flags.begin(), row_flags.end()); + for (const auto& entry : lifted_entries) { + const i_t row = std::get<0>(entry); + if (strengthening.row_orientation(row) == 1) { + rhs[row] = 0.0; + } else { + lhs[row] = 0.0; + } + } + papilo::RowFlags implied_flags; + implied_flags.set(papilo::RowFlag::kLhsInf); + lhs.resize(num_rows + num_implied, 0.0); + rhs.resize(num_rows + num_implied, 0.0); + flags.resize(num_rows + num_implied, implied_flags); + + // Name the implied rows if the problem carries constraint names. + const auto& constraint_names = problem.getConstraintNames(); + if (!constraint_names.empty()) { + papilo::Vec names; + names.reserve(constraint_names.size() + num_implied); + names.assign(constraint_names.begin(), constraint_names.end()); + for (i_t k = 0; k < num_implied; ++k) { + names.push_back("implied_indicator_" + std::to_string(k)); + } + problem.setConstraintNames(std::move(names)); + } + + // Replace the constraint matrix with the strengthened one. + papilo::SparseStorage storage( + std::move(entries), num_rows + num_implied, num_cols, true, 4.0, 30); + problem.setConstraintMatrix(std::move(storage), std::move(lhs), std::move(rhs), std::move(flags)); + + CUOPT_LOG_DEBUG( + "Indicator strengthening: %d implied indicator rows added, %d capacity rows lifted over %d " + "variable upper bounds", + num_implied, + num_lifted, + strengthening.num_variable_upper_bounds()); +} + +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT +template void strengthen_indicators(papilo::Problem&); +#endif + +#if MIP_INSTANTIATE_DOUBLE +template void strengthen_indicators(papilo::Problem&); +#endif + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/indicator_strengthening.hpp b/cpp/src/mip_heuristics/presolve/indicator_strengthening.hpp new file mode 100644 index 0000000000..203447f2df --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/indicator_strengthening.hpp @@ -0,0 +1,30 @@ +/* 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 +#if !defined(__clang__) +#pragma GCC diagnostic pop +#endif + +namespace cuopt::mathematical_optimization::mip { + +// Adds an implied indicator row y <= sum_{g in D} z_g for every implication row +// y <= sum_{j in S} x_j whose members are bounded by indicators x_j <= z_g, and lifts every +// capacity row sum_{i in S} x_i - s <= K whose members share an indicator z into sum_{i in S} x_i - +// s <= K z. Both need an integral indicator, so this is for MIPs only. +template +void strengthen_indicators(papilo::Problem& problem); + +} // 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 a28711a03d..81189c8e09 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -42,6 +42,7 @@ #include #include #include +#include #include #include #include @@ -747,8 +748,8 @@ 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())); + maybe_add(uptr(new GF2Presolve())); + maybe_add(uptr(new BHWCoeffReduce())); } // fast presolvers maybe_add(uptr(new papilo::SingletonCols())); @@ -927,14 +928,21 @@ third_party_presolve_status_t third_party_presolve_t::apply_papilo( // Capture original dimensions before papilo.apply() mutates papilo_problem // in place into its reduced form. - const i_t original_n_vars = static_cast(papilo_problem.getNCols()); - const i_t original_n_cons = static_cast(papilo_problem.getNRows()); - const i_t original_nnz = static_cast(papilo_problem.getConstraintMatrix().getNnz()); + const i_t original_n_vars = papilo_problem.getNCols(); + const i_t original_n_cons = papilo_problem.getNRows(); + const i_t original_nnz = papilo_problem.getConstraintMatrix().getNnz(); CUOPT_LOG_DEBUG("Original problem: %d constraints, %d variables, %d nonzeros", original_n_cons, original_n_vars, original_nnz); + + if (category == problem_category_t::MIP && indicator_strengthening_ && + (!reduction_allowlist_.has_value() || + reduction_allowlist_->count("indicatorstrengthening") > 0)) { + strengthen_indicators(papilo_problem); + } + CUOPT_LOG_INFO("\nRunning Papilo presolve (git hash %s)", PAPILO_GITHASH); if (category == problem_category_t::MIP) { dual_postsolve = false; } papilo::Presolve papilo_presolver; diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp index f1426749c3..8dfb69fe3f 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp @@ -115,6 +115,8 @@ class third_party_presolve_t { reduction_allowlist_ = std::move(allowlist); } + void set_indicator_strengthening(bool enabled) { indicator_strengthening_ = enabled; } + // Apply the presolve on an simplex::user_problem in-place. Used in sub MIP and (in the future) // restarts. third_party_presolve_status_t apply_to_subproblem( @@ -225,6 +227,7 @@ class third_party_presolve_t { f_t original_objective_scaling_factor_{1}; std::optional> reduction_allowlist_{}; + bool indicator_strengthening_{true}; }; // Just for testing the conversion: user_problem -> Papilo problem -> user_problem. diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 263122bd03..e6b6243d96 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -637,7 +637,8 @@ mip_solution_t solve_mip_helper( ? std::numeric_limits::infinity() : timer.remaining_time(); - presolver = std::make_unique>(); + presolver = std::make_unique>(); + presolver->set_indicator_strengthening(settings.indicator_strengthening); auto result = presolver->apply_presolve_from_op_problem( op_problem, cuopt::mathematical_optimization::problem_category_t::MIP, diff --git a/cpp/tests/internal/CMakeLists.txt b/cpp/tests/internal/CMakeLists.txt index 3cd1235717..3a5293a6d2 100644 --- a/cpp/tests/internal/CMakeLists.txt +++ b/cpp/tests/internal/CMakeLists.txt @@ -32,6 +32,7 @@ ConfigureTest(NUMOPT_INTERNAL_TEST ${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/indicator_strengthening_test.cpp ${CUOPT_TEST_DIR}/mip/arc_flow_test.cu ${CUOPT_TEST_DIR}/mip/single_lock_dual_aggregation_test.cpp ${CUOPT_TEST_DIR}/mip/termination_test.cu diff --git a/cpp/tests/mip/indicator_strengthening_test.cpp b/cpp/tests/mip/indicator_strengthening_test.cpp new file mode 100644 index 0000000000..9a90df418e --- /dev/null +++ b/cpp/tests/mip/indicator_strengthening_test.cpp @@ -0,0 +1,326 @@ +/* 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::mip::test { + +namespace { + +constexpr double kInf = std::numeric_limits::infinity(); +constexpr double kTol = 1e-9; + +using entries_t = std::vector>; + +struct row_t { + entries_t entries; + double lhs; + double rhs; +}; + +// Every column has a lower bound of 0. A row side of +-kInf is left unbounded. +papilo::Problem build_problem(const std::vector& upper_bounds, + const std::vector& integral, + const std::vector& rows) +{ + const int num_cols = upper_bounds.size(); + const int num_rows = rows.size(); + + papilo::ProblemBuilder builder; + builder.setNumCols(num_cols); + builder.setNumRows(num_rows); + for (int col = 0; col < num_cols; ++col) { + builder.setObj(col, 0.0); + builder.setColLb(col, 0.0); + builder.setColUb(col, upper_bounds[col]); + builder.setColIntegral(col, integral[col]); + } + for (int row = 0; row < num_rows; ++row) { + for (const auto& [col, value] : rows[row].entries) { + builder.addEntry(row, col, value); + } + builder.setRowLhsInf(row, rows[row].lhs == -kInf); + builder.setRowRhsInf(row, rows[row].rhs == kInf); + if (rows[row].lhs != -kInf) { builder.setRowLhs(row, rows[row].lhs); } + if (rows[row].rhs != kInf) { builder.setRowRhs(row, rows[row].rhs); } + } + return builder.build(); +} + +// The stored (column, value) pairs of a row, in storage order. +entries_t row_entries(const papilo::Problem& problem, int row) +{ + const auto coefficients = problem.getConstraintMatrix().getRowCoefficients(row); + const int len = coefficients.getLength(); + entries_t entries; + entries.reserve(len); + for (int p = 0; p < len; ++p) { + entries.emplace_back(coefficients.getIndices()[p], coefficients.getValues()[p]); + } + return entries; +} + +double row_activity(const papilo::Problem& problem, + int row, + const std::vector& point) +{ + double activity = 0.0; + for (const auto& [col, value] : row_entries(problem, row)) { + activity += value * point[col]; + } + return activity; +} + +bool is_feasible(const papilo::Problem& problem, const std::vector& point) +{ + const auto& matrix = problem.getConstraintMatrix(); + const auto& flags = matrix.getRowFlags(); + for (int row = 0; row < matrix.getNRows(); ++row) { + const double activity = row_activity(problem, row, point); + if (!flags[row].test(papilo::RowFlag::kLhsInf) && + activity < matrix.getLeftHandSides()[row] - kTol) { + return false; + } + if (!flags[row].test(papilo::RowFlag::kRhsInf) && + activity > matrix.getRightHandSides()[row] + kTol) { + return false; + } + } + return true; +} + +// Every integer point in the box [0, upper_bounds]. The continuous columns in these tests only meet +// integral row sides, so the integer grid reaches every vertex of their feasible intervals. +std::vector> integer_points(const std::vector& upper_bounds) +{ + const int num_cols = upper_bounds.size(); + size_t num_points = 1; + for (int upper_bound : upper_bounds) { + num_points *= upper_bound + 1; + } + std::vector> points(num_points, std::vector(num_cols)); + for (size_t k = 0; k < num_points; ++k) { + size_t rest = k; + for (int col = 0; col < num_cols; ++col) { + points[k][col] = rest % (upper_bounds[col] + 1); + rest /= upper_bounds[col] + 1; + } + } + return points; +} + +} // namespace + +// Reduction 1 alone. Two groups, z0 owning {x0, x1} and z1 owning {x2, x3}. +// +// y0 <= x0 + x1 + x2 members span {z0, z1} => y0 <= z0 + z1 +// x2 + x3 - y1 >= 0 members span {z1} => y1 <= z1 (stored as >=) +// y2 <= x0 + x2 one indicator per member, nothing to aggregate +TEST(IndicatorStrengthening, ImpliedIndicatorRows) +{ + constexpr int z0 = 0; + constexpr int z1 = 1; + constexpr int x0 = 2; + constexpr int x1 = 3; + constexpr int x2 = 4; + constexpr int x3 = 5; + constexpr int y0 = 6; + constexpr int y1 = 7; + constexpr int y2 = 8; + + const std::vector upper_bounds(9, 1); + const std::vector integral(9, true); + const std::vector rows{ + {{{x0, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x1, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x2, 1.0}, {z1, -1.0}}, -kInf, 0.0}, + {{{x3, 1.0}, {z1, -1.0}}, -kInf, 0.0}, + {{{y0, 1.0}, {x0, -1.0}, {x1, -1.0}, {x2, -1.0}}, -kInf, 0.0}, + {{{x2, 1.0}, {x3, 1.0}, {y1, -1.0}}, 0.0, kInf}, + {{{y2, 1.0}, {x0, -1.0}, {x2, -1.0}}, -kInf, 0.0}, + }; + const int num_rows = rows.size(); + + const auto original = build_problem(upper_bounds, integral, rows); + auto strengthened = build_problem(upper_bounds, integral, rows); + strengthen_indicators(strengthened); + + const auto& matrix = strengthened.getConstraintMatrix(); + ASSERT_EQ(matrix.getNRows(), num_rows + 2); + for (int row = 0; row < num_rows; ++row) { + EXPECT_EQ(row_entries(strengthened, row), row_entries(original, row)) << "row " << row; + } + + const int implied_y0 = num_rows; + const int implied_y1 = num_rows + 1; + EXPECT_EQ(row_entries(strengthened, implied_y0), (entries_t{{z0, -1.0}, {z1, -1.0}, {y0, 1.0}})); + EXPECT_EQ(row_entries(strengthened, implied_y1), (entries_t{{z1, -1.0}, {y1, 1.0}})); + for (int row : {implied_y0, implied_y1}) { + EXPECT_TRUE(matrix.getRowFlags()[row].test(papilo::RowFlag::kLhsInf)) << "row " << row; + EXPECT_FALSE(matrix.getRowFlags()[row].test(papilo::RowFlag::kRhsInf)) << "row " << row; + EXPECT_EQ(matrix.getRightHandSides()[row], 0.0) << "row " << row; + } + + for (const auto& point : integer_points(upper_bounds)) { + EXPECT_EQ(is_feasible(strengthened, point), is_feasible(original, point)); + } + + // LP point: y0 = 1 is covered by x0 + x1 + x2 = 1, but z0 + z1 = 2/3. + std::vector fractional(9, 0.0); + fractional[z0] = fractional[z1] = 1.0 / 3.0; + fractional[x0] = fractional[x1] = fractional[x2] = 1.0 / 3.0; + fractional[y0] = 1.0; + EXPECT_TRUE(is_feasible(original, fractional)); + EXPECT_GT(row_activity(strengthened, implied_y0, fractional), kTol); +} + +// Reduction 2 alone. z0 owns {x0, x1, x2}, z1 owns {x3}, s is a continuous slack in [0, 2]. +// +// x0 + x1 + x2 - s <= 1 => x0 + x1 + x2 - s <= z0 +// -x0 - x1 - x2 >= -2 => -x0 - x1 - x2 + 2 z0 >= 0 +// x0 + x3 <= 1 members share no indicator, unchanged +// x0 + x1 <= 2 capacity is not binding, unchanged +TEST(IndicatorStrengthening, LiftedCapacityRows) +{ + constexpr int z0 = 0; + constexpr int z1 = 1; + constexpr int x0 = 2; + constexpr int x1 = 3; + constexpr int x2 = 4; + constexpr int x3 = 5; + constexpr int s = 6; + + const std::vector upper_bounds{1, 1, 1, 1, 1, 1, 2}; + const std::vector integral{true, true, true, true, true, true, false}; + const std::vector rows{ + {{{x0, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x1, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x2, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x3, 1.0}, {z1, -1.0}}, -kInf, 0.0}, + {{{x0, 1.0}, {x1, 1.0}, {x2, 1.0}, {s, -1.0}}, -kInf, 1.0}, + {{{x0, -1.0}, {x1, -1.0}, {x2, -1.0}}, -2.0, kInf}, + {{{x0, 1.0}, {x3, 1.0}}, -kInf, 1.0}, + {{{x0, 1.0}, {x1, 1.0}}, -kInf, 2.0}, + }; + const int num_rows = rows.size(); + const int lifted_le = 4; + const int lifted_ge = 5; + + const auto original = build_problem(upper_bounds, integral, rows); + auto strengthened = build_problem(upper_bounds, integral, rows); + strengthen_indicators(strengthened); + + const auto& matrix = strengthened.getConstraintMatrix(); + ASSERT_EQ(matrix.getNRows(), num_rows); + for (int row = 0; row < num_rows; ++row) { + if (row == lifted_le || row == lifted_ge) { continue; } + EXPECT_EQ(row_entries(strengthened, row), row_entries(original, row)) << "row " << row; + EXPECT_EQ(matrix.getRightHandSides()[row], + original.getConstraintMatrix().getRightHandSides()[row]) + << "row " << row; + } + + EXPECT_EQ(row_entries(strengthened, lifted_le), + (entries_t{{z0, -1.0}, {x0, 1.0}, {x1, 1.0}, {x2, 1.0}, {s, -1.0}})); + EXPECT_TRUE(matrix.getRowFlags()[lifted_le].test(papilo::RowFlag::kLhsInf)); + EXPECT_EQ(matrix.getRightHandSides()[lifted_le], 0.0); + + EXPECT_EQ(row_entries(strengthened, lifted_ge), + (entries_t{{z0, 2.0}, {x0, -1.0}, {x1, -1.0}, {x2, -1.0}})); + EXPECT_TRUE(matrix.getRowFlags()[lifted_ge].test(papilo::RowFlag::kRhsInf)); + EXPECT_EQ(matrix.getLeftHandSides()[lifted_ge], 0.0); + + for (const auto& point : integer_points(upper_bounds)) { + EXPECT_EQ(is_feasible(strengthened, point), is_feasible(original, point)); + } + + // LP point: x0 + x1 + x2 = 1 fits the capacity, but z0 = 1/3 does not pay for it. + std::vector fractional(7, 0.0); + fractional[z0] = 1.0 / 3.0; + fractional[x0] = fractional[x1] = fractional[x2] = 1.0 / 3.0; + EXPECT_TRUE(is_feasible(original, fractional)); + EXPECT_GT(row_activity(strengthened, lifted_le, fractional), kTol); + EXPECT_LT(row_activity(strengthened, lifted_ge, fractional), -kTol); +} + +// Both reductions on one set-cover-like model: z0 owns {x0, x1}, z1 owns {x2, x3}, y must be +// covered by one of the members, and each group has a capacity row. +// +// y <= x0 + x1 + x2 + x3 => appended y <= z0 + z1 +// x0 + x1 - s <= 1 => x0 + x1 - s <= z0 +// x2 + x3 <= 1 => x2 + x3 <= z1 +TEST(IndicatorStrengthening, ImpliedAndLiftedRows) +{ + constexpr int z0 = 0; + constexpr int z1 = 1; + constexpr int x0 = 2; + constexpr int x1 = 3; + constexpr int x2 = 4; + constexpr int x3 = 5; + constexpr int y = 6; + constexpr int s = 7; + + const std::vector upper_bounds{1, 1, 1, 1, 1, 1, 1, 2}; + const std::vector integral{true, true, true, true, true, true, true, false}; + const std::vector rows{ + {{{x0, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x1, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x2, 1.0}, {z1, -1.0}}, -kInf, 0.0}, + {{{x3, 1.0}, {z1, -1.0}}, -kInf, 0.0}, + {{{y, 1.0}, {x0, -1.0}, {x1, -1.0}, {x2, -1.0}, {x3, -1.0}}, -kInf, 0.0}, + {{{x0, 1.0}, {x1, 1.0}, {s, -1.0}}, -kInf, 1.0}, + {{{x2, 1.0}, {x3, 1.0}}, -kInf, 1.0}, + }; + const int num_rows = rows.size(); + const int lifted_z0 = 5; + const int lifted_z1 = 6; + const int implied_y = num_rows; + + const auto original = build_problem(upper_bounds, integral, rows); + auto strengthened = build_problem(upper_bounds, integral, rows); + strengthen_indicators(strengthened); + + const auto& matrix = strengthened.getConstraintMatrix(); + ASSERT_EQ(matrix.getNRows(), num_rows + 1); + for (int row = 0; row < lifted_z0; ++row) { + EXPECT_EQ(row_entries(strengthened, row), row_entries(original, row)) << "row " << row; + } + + EXPECT_EQ(row_entries(strengthened, lifted_z0), + (entries_t{{z0, -1.0}, {x0, 1.0}, {x1, 1.0}, {s, -1.0}})); + EXPECT_EQ(row_entries(strengthened, lifted_z1), (entries_t{{z1, -1.0}, {x2, 1.0}, {x3, 1.0}})); + EXPECT_EQ(row_entries(strengthened, implied_y), (entries_t{{z0, -1.0}, {z1, -1.0}, {y, 1.0}})); + for (int row : {lifted_z0, lifted_z1, implied_y}) { + EXPECT_TRUE(matrix.getRowFlags()[row].test(papilo::RowFlag::kLhsInf)) << "row " << row; + EXPECT_FALSE(matrix.getRowFlags()[row].test(papilo::RowFlag::kRhsInf)) << "row " << row; + EXPECT_EQ(matrix.getRightHandSides()[row], 0.0) << "row " << row; + } + + for (const auto& point : integer_points(upper_bounds)) { + EXPECT_EQ(is_feasible(strengthened, point), is_feasible(original, point)); + } + + // LP point: every original row holds, and each of the three new rows is violated by 1/3. + std::vector fractional(8, 0.0); + fractional[z0] = fractional[z1] = 1.0 / 3.0; + fractional[x0] = fractional[x1] = fractional[x2] = fractional[x3] = 1.0 / 3.0; + fractional[y] = 1.0; + EXPECT_TRUE(is_feasible(original, fractional)); + for (int row : {lifted_z0, lifted_z1, implied_y}) { + EXPECT_GT(row_activity(strengthened, row, fractional), kTol) << "row " << row; + } +} + +} // namespace cuopt::mathematical_optimization::mip::test