Block Bounded Value Elimination presolving pass for binary problems

ci/validate_wheel.sh

@@ -22,7 +22,7 @@ PYDISTCHECK_ARGS=(
 if [[ "${package_dir}" == "python/libcuopt" ]]; then
     if [[ "${RAPIDS_CUDA_MAJOR}" == "12" ]]; then
         PYDISTCHECK_ARGS+=(
-            --max-allowed-size-compressed '690Mi'
+            --max-allowed-size-compressed '695Mi'
         )
     else
         PYDISTCHECK_ARGS+=(

cpp/include/cuopt/mathematical_optimization/constants.h

@@ -144,6 +144,9 @@
 #define CUOPT_MIP_HYPER_SUBMIP_ITERATION_LIMIT_RATIO "mip_hyper_submip_iteration_limit_ratio"
 #define CUOPT_MIP_HYPER_SUBMIP_ENABLE_CPUFJ          "mip_hyper_submip_enable_cpufj"
 
+/* @brief Block bounded-variable-elimination step of cuOpt's internal MIP presolve */
+#define CUOPT_MIP_HYPER_BLOCK_BVE "mip_hyper_block_bve"
+
 /* @brief QCQP (barrier) scaling hyper-parameters */
 #define CUOPT_QCQP_HYPER_RUIZ_EQUILIBRATION "qcqp_hyper_ruiz_equilibration"
 

cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp

@@ -160,6 +160,15 @@ class mip_solver_settings_t {
    * When this is `false`, probing is skipped even if presolve is otherwise on.
    */
   bool probing{true};
+  /**
+   * @brief Enable the block bounded-variable-elimination step of cuOpt's MIP presolve.
+   *
+   * Runs after trivial_presolve and eliminates blocks of functionally-determined binary auxiliary
+   * variables discovered via the probing-cache implication closure, re-encoding each block's
+   * projected relation as certified prime-implicate clauses. Requires the probing-cache step; a
+   * no-op when no certified reduction exists.
+   */
+  bool block_bve{true};
   /**
    * @brief Determinism mode for MIP solver.
    *

cpp/src/io/mps_writer.cpp

@@ -228,8 +228,8 @@ void mps_writer_t<i_t, f_t>::write(const std::string& mps_file_path)
   // save coefficients with full precision
   mps_file << std::setprecision(std::numeric_limits<f_t>::max_digits10);
 
-  // NAME section
-  mps_file << "NAME          " << problem_.get_problem_name() << "\n";
+  const std::string& pname = problem_.get_problem_name();
+  mps_file << "NAME          " << (pname.empty() ? "cuopt" : pname) << "\n";
 
   if (problem_.get_sense()) { mps_file << "OBJSENSE\n MAXIMIZE\n"; }
 

cpp/src/math_optimization/solver_settings.cu

@@ -211,6 +211,9 @@ solver_settings_t<i_t, f_t>::solver_settings_t() : pdlp_settings(), mip_settings
     {CUOPT_MIP_HYPER_DIVING_SHOW_TYPE, &mip_settings.diving_params.show_type, false, "log diving heuristic type when it finds a new incumbent"},
     // 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"},
+    // Kept a hyper-parameter while block-BVE bakes in: settable so a run can be bisected against it,
+    // but not yet a documented knob (no constant in the proto / server surfaces).
+    {CUOPT_MIP_HYPER_BLOCK_BVE, &mip_settings.block_bve, true, "eliminate blocks of binaries in cuOpt's MIP presolve (needs " CUOPT_MIP_PROBING ")"},
   };
   // String parameters
   string_parameters = {

cpp/src/mip_heuristics/CMakeLists.txt

@@ -33,6 +33,7 @@ set(MIP_NON_LP_FILES
   ${CMAKE_CURRENT_SOURCE_DIR}/local_search/rounding/simple_rounding.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/local_search/feasibility_pump/feasibility_pump.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/local_search/line_segment_search/line_segment_search.cu
+  ${CMAKE_CURRENT_SOURCE_DIR}/presolve/block_bve.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bounds_presolve.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bounds_update_data.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/presolve/semi_continuous.cu

cpp/src/mip_heuristics/diversity/diversity_manager.cu

@@ -11,13 +11,18 @@
 #include <mip_heuristics/mip_constants.hpp>
 #include <mip_heuristics/presolve/third_party_presolve.hpp>
 
+#include <mip_heuristics/presolve/block_bve.cuh>
 #include <mip_heuristics/presolve/conflict_graph/clique_table.cuh>
 #include <mip_heuristics/presolve/probing_cache.cuh>
 #include <mip_heuristics/presolve/trivial_presolve.cuh>
 #include <mip_heuristics/problem/problem_helpers.cuh>
 
 #include <pdlp/solve.cuh>
 
+#include <cuopt/mathematical_optimization/io/mps_data_model.hpp>
+#include <cuopt/mathematical_optimization/io/mps_writer.hpp>
+#include <cuopt/mathematical_optimization/optimization_problem_utils.hpp>
+#include <utilities/copy_helpers.hpp>
 #include <utilities/scope_guard.hpp>
 
 #include <memory>
@@ -39,6 +44,59 @@ size_t sub_mip_recombiner_config_t::max_n_of_vars_from_other =
 template <typename i_t, typename f_t>
 std::vector<recombiner_enum_t> recombiner_t<i_t, f_t>::enabled_recombiners;
 
+template <typename i_t, typename f_t>
+static cuopt::mathematical_optimization::io::mps_data_model_t<i_t, f_t> problem_to_mps_data_model(
+  const problem_t<i_t, f_t>& problem)
+{
+  auto stream = problem.handle_ptr->get_stream();
+  auto h_off  = cuopt::host_copy(problem.offsets, stream);
+  auto h_ind  = cuopt::host_copy(problem.variables, stream);
+  auto h_val  = cuopt::host_copy(problem.coefficients, stream);
+  auto h_clb  = cuopt::host_copy(problem.constraint_lower_bounds, stream);
+  auto h_cub  = cuopt::host_copy(problem.constraint_upper_bounds, stream);
+  auto h_obj  = cuopt::host_copy(problem.objective_coefficients, stream);
+  auto h_vb   = cuopt::host_copy(problem.variable_bounds, stream);
+  auto h_vt   = cuopt::host_copy(problem.variable_types, stream);
+  problem.handle_ptr->sync_stream();
+
+  const i_t n_vars = problem.n_variables;
+  std::vector<f_t> var_lower(n_vars), var_upper(n_vars);
+  for (i_t v = 0; v < n_vars; ++v) {
+    var_lower[v] = get_lower(h_vb[v]);
+    var_upper[v] = get_upper(h_vb[v]);
+  }
+  std::vector<char> var_types(n_vars);
+  for (i_t v = 0; v < n_vars; ++v)
+    var_types[v] = var_type_to_char(h_vt[v]);
+
+  cuopt::mathematical_optimization::io::mps_data_model_t<i_t, f_t> model;
+  model.set_maximize(false);
+  if (!h_off.empty()) {
+    model.set_csr_constraint_matrix(std::span<const f_t>{h_val.data(), h_val.size()},
+                                    std::span<const i_t>{h_ind.data(), h_ind.size()},
+                                    std::span<const i_t>{h_off.data(), h_off.size()});
+  }
+  if (problem.n_constraints != 0) {
+    model.set_constraint_lower_bounds(std::span<const f_t>{h_clb.data(), h_clb.size()});
+    model.set_constraint_upper_bounds(std::span<const f_t>{h_cub.data(), h_cub.size()});
+  }
+  if (n_vars != 0) {
+    model.set_objective_coefficients(std::span<const f_t>{h_obj.data(), h_obj.size()});
+    model.set_variable_lower_bounds(std::span<const f_t>{var_lower.data(), var_lower.size()});
+    model.set_variable_upper_bounds(std::span<const f_t>{var_upper.data(), var_upper.size()});
+    model.set_variable_types(var_types);
+  }
+  model.set_objective_scaling_factor(problem.presolve_data.objective_scaling_factor);
+  model.set_objective_offset(problem.presolve_data.objective_offset);
+  if (problem.original_problem_ptr != nullptr &&
+      !problem.original_problem_ptr->get_problem_name().empty()) {
+    model.set_problem_name(problem.original_problem_ptr->get_problem_name());
+  } else {
+    model.set_problem_name("cuopt");
+  }
+  return model;
+}
+
 template <typename i_t, typename f_t>
 diversity_manager_t<i_t, f_t>::diversity_manager_t(mip_solver_context_t<i_t, f_t>& context_)
   : context(context_),
@@ -279,6 +337,42 @@ void diversity_manager_t<i_t, f_t>::add_user_given_solutions(
   }
 }
 
+// Pin variables that a BVE projection table showed to have a single admissible value. Ids arrive in
+// the original frame and may repeat across blocks and rounds. Returns false when two blocks
+// disagree on a variable, which proves infeasibility since each fixing is a consequence of its
+// block alone.
+template <typename i_t, typename f_t>
+static bool apply_bve_fixings(problem_t<i_t, f_t>& problem,
+                              const std::vector<std::pair<i_t, bool>>& fixings,
+                              i_t& n_applied)
+{
+  n_applied = 0;
+  if (fixings.empty()) { return true; }
+  std::vector<std::pair<i_t, bool>> sorted(fixings);
+  std::sort(sorted.begin(), sorted.end());
+
+  const std::vector<i_t>& reverse_original_ids = problem.reverse_original_ids;
+  std::vector<i_t> var_indices;
+  std::vector<f_t> lb_values;
+  std::vector<f_t> ub_values;
+  for (size_t k = 0; k < sorted.size(); ++k) {
+    const auto [original_id, value] = sorted[k];
+    if (k > 0 && original_id == sorted[k - 1].first) {
+      if (value != sorted[k - 1].second) { return false; }
+      continue;
+    }
+    if (original_id < 0 || original_id >= (i_t)reverse_original_ids.size()) { continue; }
+    const i_t column = reverse_original_ids[original_id];
+    if (column < 0 || column >= problem.n_variables) { continue; }  // already eliminated
+    var_indices.push_back(column);
+    lb_values.push_back(value ? f_t(1) : f_t(0));
+    ub_values.push_back(value ? f_t(1) : f_t(0));
+  }
+  n_applied = (i_t)var_indices.size();
+  problem.update_variable_bounds(var_indices, lb_values, ub_values);
+  return true;
+}
+
 template <typename i_t, typename f_t>
 bool diversity_manager_t<i_t, f_t>::run_presolve(f_t time_limit, timer_t global_timer)
 {
@@ -304,19 +398,102 @@ bool diversity_manager_t<i_t, f_t>::run_presolve(f_t time_limit, timer_t global_
     CUOPT_LOG_INFO("Probing-cache step disabled via %s=false", CUOPT_MIP_PROBING);
     run_probing_cache = false;
   }
-  if (run_probing_cache) {
-    // Run probing cache before trivial presolve to discover variable implications
+  const bool remap_cache_ids           = true;
+  problem_ptr->related_vars_time_limit = context.settings.heuristic_params.related_vars_time_limit;
+
+  if (run_probing_cache && !global_timer.check_time_limit() && !presolve_timer.check_time_limit()) {
     const f_t max_time_on_probing = diversity_config.max_time_on_probing;
-    f_t time_for_probing_cache    = std::min(max_time_on_probing, time_limit);
+    f_t time_for_probing_cache =
+      std::min(max_time_on_probing, std::min(time_limit, (f_t)presolve_timer.remaining_time()));
     timer_t probing_timer{time_for_probing_cache};
-    // this function computes probing cache, finds singletons, substitutions and changes the problem
     bool problem_is_infeasible =
       compute_probing_cache(ls.constraint_prop.bounds_update, *problem_ptr, probing_timer);
     if (problem_is_infeasible) { return false; }
   }
-  const bool remap_cache_ids           = true;
-  problem_ptr->related_vars_time_limit = context.settings.heuristic_params.related_vars_time_limit;
+
   if (!global_timer.check_time_limit()) { trivial_presolve(*problem_ptr, remap_cache_ids); }
+
+  i_t max_bve_rounds          = 3;
+  const i_t n_vars_before_bve = problem_ptr->n_variables;
+  const i_t n_rows_before_bve = problem_ptr->n_constraints;
+
+  if (!run_probing_cache) max_bve_rounds = 0;
+  // Implications read off the projection tables, accumulated across rounds. They feed the next
+  // round's adjacency (pairs the cache never held) and are folded back into the cache afterwards.
+  probe_findings_t<i_t> bve_findings;
+  for (i_t bve_round = 0; bve_round < max_bve_rounds; ++bve_round) {
+    if (!context.settings.block_bve || problem_ptr->empty || global_timer.check_time_limit() ||
+        presolve_timer.check_time_limit()) {
+      break;
+    }
+
+    const i_t n_vars_before = problem_ptr->n_variables;
+    const i_t n_rows_before = problem_ptr->n_constraints;
+    auto impl_adj           = bve_build_impl_adj(ls.constraint_prop.bounds_update.probing_cache,
+                                       problem_ptr->reverse_original_ids,
+                                       problem_ptr->n_variables,
+                                       &bve_findings);
+    double bve_work_units   = 0.0;
+    timer_t bve_timer(global_timer.clamp_remaining_time(presolve_timer.remaining_time()));
+    const bool reduced =
+      block_bve_presolve(*problem_ptr, impl_adj, bve_timer, bve_work_units, &bve_findings);
+    CUOPT_LOG_DEBUG("Block-BVE outer round %d/%d: reduced=%d vars %d->%d rows %d->%d",
+                    bve_round + 1,
+                    max_bve_rounds,
+                    (int)reduced,
+                    n_vars_before,
+                    problem_ptr->n_variables,
+                    n_rows_before,
+                    problem_ptr->n_constraints);
+    if (!reduced) { break; }
+    if (problem_ptr->n_variables >= n_vars_before) { break; }
+  }
+
+  // Harvest the projections: tighten the cache in place, pin the variables the blocks left with a
+  // single value, then propagate.
+  ls.constraint_prop.bounds_update.probing_cache.merge_forcings(bve_findings.forcings,
+                                                                bve_findings.fixings);
+  i_t n_bve_fixings = 0;
+  if (!global_timer.check_time_limit()) {
+    if (!apply_bve_fixings(*problem_ptr, bve_findings.fixings, n_bve_fixings)) {
+      stats.presolve_time = timer.elapsed_time();
+      return false;
+    }
+    if (n_bve_fixings > 0) { trivial_presolve(*problem_ptr, remap_cache_ids); }
+  }
+  const bool bve_changed_model = problem_ptr->n_variables != n_vars_before_bve ||
+                                 problem_ptr->n_constraints != n_rows_before_bve ||
+                                 n_bve_fixings > 0;
+  if (bve_changed_model) {
+    CUOPT_LOG_DEBUG("Block-BVE projections fixed %d variables", n_bve_fixings);
+    // propagate fixings if any
+    if (!problem_ptr->empty && !global_timer.check_time_limit()) {
+      ls.constraint_prop.bounds_update.resize(*problem_ptr);
+      auto bve_term_crit = ls.constraint_prop.bounds_update.solve(*problem_ptr);
+      if (ls.constraint_prop.bounds_update.infeas_constraints_count > 0) {
+        stats.presolve_time = timer.elapsed_time();
+        return false;
+      }
+      if (termination_criterion_t::NO_UPDATE != bve_term_crit) {
+        ls.constraint_prop.bounds_update.set_updated_bounds(*problem_ptr);
+      }
+    }
+  }
+
+  if (const char* export_flag = std::getenv("CUOPT_EXPORT_GPU_PRESOLVED_PROBLEM");
+      export_flag != nullptr && std::atoi(export_flag) != 0) {
+    const std::string instance_name =
+      (problem_ptr->original_problem_ptr != nullptr &&
+       !problem_ptr->original_problem_ptr->get_problem_name().empty())
+        ? problem_ptr->original_problem_ptr->get_problem_name()
+        : std::string("cuopt");
+    const std::string mps_path = instance_name + "_gpupresolved.mps";
+    CUOPT_LOG_DEBUG("Exporting GPU-presolved problem to %s", mps_path.c_str());
+    auto model = problem_to_mps_data_model(*problem_ptr);
+    cuopt::mathematical_optimization::io::mps_writer_t<i_t, f_t> writer(model);
+    writer.write(mps_path);
+    exit(0);
+  }
   if (!problem_ptr->empty && !check_bounds_sanity(*problem_ptr)) { return false; }
   // if (!presolve_timer.check_time_limit() && !context.settings.heuristics_only &&
   //     !problem_ptr->empty) {

cpp/src/mip_heuristics/diversity/population.cu

@@ -233,6 +233,12 @@ std::vector<solution_t<i_t, f_t>> population_t<i_t, f_t>::get_external_solutions
           sol.compute_number_of_integers(),
           problem_ptr->n_integer_vars);
       }
+      if (std::abs(sol.get_objective() - h_entry.objective) > OBJECTIVE_EPSILON) {
+        CUOPT_LOG_DEBUG(
+          "External solution objective mismatch: sol.get_objective() = %g, h_entry.objective = %g",
+          sol.get_objective(),
+          h_entry.objective);
+      }
       sol.handle_ptr->sync_stream();
       return_vector.emplace_back(std::move(sol));
       counter++;