Recursive subMIP

cpp/src/branch_and_bound/CMakeLists.txt

@@ -4,11 +4,12 @@
 # cmake-format: on
 
 set(BRANCH_AND_BOUND_SRC_FILES
-  ${CMAKE_CURRENT_SOURCE_DIR}/branch_and_bound.cpp
-  ${CMAKE_CURRENT_SOURCE_DIR}/pseudo_costs.cpp
-  ${CMAKE_CURRENT_SOURCE_DIR}/diving_heuristics.cpp
-  )
+        ${CMAKE_CURRENT_SOURCE_DIR}/branch_and_bound.cpp
+        ${CMAKE_CURRENT_SOURCE_DIR}/pseudo_costs.cpp
+        ${CMAKE_CURRENT_SOURCE_DIR}/diving_heuristics.cpp
+        ${CMAKE_CURRENT_SOURCE_DIR}/presolve.cpp
+)
 
 
 set(CUOPT_SRC_FILES ${CUOPT_SRC_FILES}
-  ${BRANCH_AND_BOUND_SRC_FILES} PARENT_SCOPE)
+        ${BRANCH_AND_BOUND_SRC_FILES} PARENT_SCOPE)

cpp/src/branch_and_bound/branch_and_bound.hpp

@@ -12,6 +12,7 @@
 #include <branch_and_bound/mip_node.hpp>
 #include <branch_and_bound/node_queue.hpp>
 #include <branch_and_bound/pseudo_costs.hpp>
+#include <branch_and_bound/submip.hpp>
 #include <branch_and_bound/worker.hpp>
 #include <branch_and_bound/worker_pool.hpp>
 
@@ -41,26 +42,13 @@
 #include <vector>
 
 namespace cuopt::mathematical_optimization::mip {
+
 template <typename i_t, typename f_t>
 struct clique_table_t;
-}
-
-namespace cuopt::mathematical_optimization::mip {
 
 template <typename i_t, typename f_t>
 struct mip_symmetry_t;
 
-enum class mip_status_t {
-  OPTIMAL    = 0,  // The optimal integer solution was found
-  UNBOUNDED  = 1,  // The problem is unbounded
-  INFEASIBLE = 2,  // The problem is infeasible
-  TIME_LIMIT = 3,  // The solver reached a time limit
-  NODE_LIMIT = 4,  // The maximum number of nodes was reached (not implemented)
-  NUMERICAL  = 5,  // The solver encountered a numerical error
-  UNSET      = 6,  // The status is not set
-  WORK_LIMIT = 7,  // The solver reached a deterministic work limit
-};
-
 template <typename i_t, typename f_t>
 void upper_bound_callback(f_t upper_bound);
 
@@ -83,9 +71,19 @@ class branch_and_bound_t {
                      std::shared_ptr<mip::clique_table_t<i_t, f_t>> clique_table = nullptr,
                      mip_symmetry_t<i_t, f_t>* symmetry                          = nullptr);
 
+  branch_and_bound_t(branch_and_bound_t& other,
+                     const simplex::simplex_solver_settings_t<i_t, f_t>& solver_settings,
+                     const std::vector<f_t>& lower,
+                     const std::vector<f_t>& upper);
+
   // Set an initial guess based on the user_problem. This should be called before solve.
   void set_initial_guess(const std::vector<f_t>& user_guess) { guess_ = user_guess; }
 
+  void set_external_upper_bound(const omp_atomic_t<f_t>* external_upper_bound)
+  {
+    external_upper_bound_ = external_upper_bound;
+  }
+
   // Set the root solution found by PDLP
   void set_root_relaxation_solution(const std::vector<f_t>& primal,
                                     const std::vector<f_t>& dual,
@@ -111,6 +109,9 @@ class branch_and_bound_t {
   // Set a solution based on the user problem during the course of the solve
   void set_solution_from_heuristics(const std::vector<f_t>& solution);
 
+  // Set a solution based on the user problem during the course of the solve
+  void set_solution_from_submip(const std::vector<f_t>& solution, f_t fixrate);
+
   // This queues the solution to be processed at the correct work unit timestamp
   void queue_external_solution_deterministic(const std::vector<f_t>& solution, double work_unit_ts);
 
@@ -125,6 +126,10 @@ class branch_and_bound_t {
   // `bound` must be in B&B's internal objective space.
   void set_initial_upper_bound(f_t bound);
 
+  void warm_start(const pseudo_costs_t<i_t, f_t>& parent_pc,
+                  const std::vector<i_t>& reduced_to_original,
+                  i_t max_samples = 1);
+
   f_t get_upper_bound() const { return upper_bound_.load(); }
   bool has_solver_space_incumbent() const { return incumbent_.has_incumbent; }
 
@@ -151,6 +156,15 @@ class branch_and_bound_t {
                                 std::vector<f_t>& lower_bounds,
                                 std::vector<f_t>& upper_bounds);
 
+  // Whether obj should replace the stored incumbent. Must be called under mutex_upper_.
+  // Compares against the stored incumbent's objective, NOT against upper_bound_, because
+  // set_initial_upper_bound can set a tighter bound from an OG-space solution that has no
+  // corresponding solver-space incumbent (e.g. papilo can't crush it back).
+  bool improves_incumbent(f_t obj) const
+  {
+    return !incumbent_.has_incumbent || obj < incumbent_.objective;
+  }
+
   // The main entry routine. Returns the solver status and populates solution with the incumbent.
   mip_status_t solve(simplex::mip_solution_t<i_t, f_t>& solution);
 
@@ -199,18 +213,12 @@ class branch_and_bound_t {
   // original-space in the mip_solver_context_t), but does NOT imply incumbent_.has_incumbent.
   omp_atomic_t<f_t> upper_bound_;
 
+  // For submip
+  const omp_atomic_t<f_t>* external_upper_bound_{nullptr};
+
   // Solver-space incumbent tracked directly by B&B.
   simplex::mip_solution_t<i_t, f_t> incumbent_;
 
-  // Whether obj should replace the stored incumbent. Must be called under mutex_upper_.
-  // Compares against the stored incumbent's objective, NOT against upper_bound_, because
-  // set_initial_upper_bound can set a tighter bound from an OG-space solution that has no
-  // corresponding solver-space incumbent (e.g. papilo can't crush it back).
-  bool improves_incumbent(f_t obj) const
-  {
-    return !incumbent_.has_incumbent || obj < incumbent_.objective;
-  }
-
   // Structure with the general info of the solver.
   branch_and_bound_stats_t<i_t, f_t> exploration_stats_;
 
@@ -233,6 +241,7 @@ class branch_and_bound_t {
   std::atomic<int> root_concurrent_halt_{0};
   std::atomic<int> node_concurrent_halt_{0};
   bool is_root_solution_set{false};
+  bool root_warm_start_{false};
 
   // Pseudocosts
   pseudo_costs_t<i_t, f_t> pc_;
@@ -246,6 +255,9 @@ class branch_and_bound_t {
   // Worker pool dedicated to diving
   diving_worker_pool_t<i_t, f_t> diving_worker_pool_;
 
+  submip_worker_pool_t<i_t, f_t> submip_worker_pool_;
+  submip_stats_t submip_stats_;
+
   // Global status of the solver.
   omp_atomic_t<mip_status_t> solver_status_;
   omp_atomic_t<bool> is_running_{false};
@@ -262,7 +274,7 @@ class branch_and_bound_t {
   std::function<void(f_t)> user_bound_callback_;
 
   void print_table_header();
-  void report_heuristic(f_t obj);
+  void report_heuristic(f_t obj, char symbol);
   void report(char symbol,
               f_t obj,
               f_t lower_bound,
@@ -316,6 +328,14 @@ class branch_and_bound_t {
   // Launch a new diving worker from a given best-first worker.
   bool launch_diving_worker(bfs_worker_t<i_t, f_t>* bfs_worker);
 
+  // If the objective is integral or must move in steps than
+  // the lower bound will be different from the leaf objective.
+  // We use the leaf objective for RINS (on_optimal_callback)
+  // and if we are integer feasible (handle_integer_solution).
+  // We use the lower bound to decide if we should fathom the
+  // node or branch.
+  void apply_objective_step(mip_node_t<i_t, f_t>* node_ptr, f_t leaf_obj);
+
   // Launch a new best-first worker from a given bfs worker.
   void launch_bfs_worker(bfs_worker_t<i_t, f_t>* worker);
 
@@ -334,6 +354,15 @@ class branch_and_bound_t {
   // to find integer feasible solutions.
   void dive_with(diving_worker_t<i_t, f_t>* worker);
 
+  void launch_submip_worker(const std::vector<f_t>& sol);
+  void solve_submip(submip_worker_t<i_t, f_t>* worker,
+                    i_t num_var_fixed,
+                    i_t num_integers,
+                    i_t submip_level,
+                    std::string_view log_prefix);
+
+  void rins(submip_worker_t<i_t, f_t>* rins_worker, const std::vector<f_t>& node_solution);
+
   // Solve the LP relaxation of a leaf node
   simplex::dual_status_t solve_node_lp(mip_node_t<i_t, f_t>* node_ptr,
                                        branch_and_bound_worker_t<i_t, f_t>* worker,

cpp/src/branch_and_bound/constants.hpp

@@ -9,7 +9,7 @@
 
 namespace cuopt::mathematical_optimization::mip {
 
-constexpr int num_search_strategies = 7;
+constexpr int num_search_strategies = 8;
 
 // Indicate the search and variable selection algorithms used by each thread
 // in B&B (See [1]).
@@ -26,7 +26,8 @@ enum search_strategy_t : int {
   GUIDED_DIVING        = 3,  // Guided diving (9.2.3).
   COEFFICIENT_DIVING   = 4,  // Coefficient diving (9.2.1)
   FARKAS_DIVING        = 5,  // Farkas Diving (see [2])
-  VECTOR_LENGTH_DIVING = 6   // Vector Length Diving (9.2.6)
+  VECTOR_LENGTH_DIVING = 6,  // Vector Length Diving (9.2.6)
+  SUBMIP               = 7,  // RINS/RENS (akin to a guided diving, see HiGHS)
 };
 
 constexpr search_strategy_t search_strategies[] = {BEST_FIRST,
@@ -35,10 +36,23 @@ constexpr search_strategy_t search_strategies[] = {BEST_FIRST,
                                                    GUIDED_DIVING,
                                                    COEFFICIENT_DIVING,
                                                    FARKAS_DIVING,
-                                                   VECTOR_LENGTH_DIVING};
+                                                   VECTOR_LENGTH_DIVING,
+                                                   SUBMIP};
 
 enum class branch_direction_t { NONE = -1, DOWN = 0, UP = 1 };
 
 enum class branch_and_bound_mode_t { PARALLEL = 0, DETERMINISTIC = 1 };
 
+enum class mip_status_t {
+  OPTIMAL     = 0,  // The optimal integer solution was found
+  UNBOUNDED   = 1,  // The problem is unbounded
+  INFEASIBLE  = 2,  // The problem is infeasible
+  TIME_LIMIT  = 3,  // The solver reached a time limit
+  NODE_LIMIT  = 4,  // The maximum number of nodes was reached (not implemented)
+  NUMERICAL   = 5,  // The solver encountered a numerical error
+  UNSET       = 6,  // The status is not set
+  WORK_LIMIT  = 7,  // The solver reached a deterministic work limit
+  SUBMIP_HALT = 8   // Stop submip solve since it no longer valid
+};
+
 }  // namespace cuopt::mathematical_optimization::mip

cpp/src/branch_and_bound/diving_heuristics.hpp

@@ -23,6 +23,7 @@ namespace cuopt::mathematical_optimization::mip {
 inline char feasible_solution_symbol(search_strategy_t strategy, bool log_diving_type)
 {
   if (strategy == BEST_FIRST) return 'B';
+  if (strategy == SUBMIP) return 'S';
   if (!log_diving_type) { return 'D'; }
   switch (strategy) {
     case COEFFICIENT_DIVING: return 'C';
@@ -47,6 +48,7 @@ bool is_search_strategy_enabled(search_strategy_t strategy,
     case COEFFICIENT_DIVING: return settings.coefficient_diving != 0;
     case FARKAS_DIVING: return settings.farkas_diving != 0;
     case VECTOR_LENGTH_DIVING: return settings.vector_length_diving != 0;
+    case SUBMIP: return false;
   }
 
   return false;

cpp/src/branch_and_bound/presolve.cpp

@@ -0,0 +1,72 @@
+/* clang-format off */
+/*
+ * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-License-Identifier: Apache-2.0
+ */
+/* clang-format on */
+
+#include <branch_and_bound/presolve.hpp>
+
+#include <algorithm>
+#include <limits>
+
+namespace cuopt::mathematical_optimization::mip {
+
+mip_status_t presolve_status_to_mip_status(third_party_presolve_status_t status)
+{
+  switch (status) {
+    case third_party_presolve_status_t::OPTIMAL: return mip_status_t::OPTIMAL;
+    case third_party_presolve_status_t::INFEASIBLE: return mip_status_t::INFEASIBLE;
+    case third_party_presolve_status_t::UNBOUNDED: return mip_status_t::UNBOUNDED;
+    case third_party_presolve_status_t::UNBNDORINFEAS: return mip_status_t::INFEASIBLE;
+    case third_party_presolve_status_t::REDUCED: return mip_status_t::UNSET;
+    case third_party_presolve_status_t::UNCHANGED: return mip_status_t::UNSET;
+  }
+  return mip_status_t::UNSET;
+}
+
+template <typename i_t, typename f_t>
+mip_status_t presolver_t<i_t, f_t>::apply(
+  simplex::user_problem_t<i_t, f_t>& problem,
+  const simplex::simplex_solver_settings_t<i_t, f_t>& settings)
+{
+  f_t presolve_time_limit = std::min(0.1 * settings.time_limit, 60.0);
+  third_party_presolve_status_t status =
+    third_party_presolver_.apply(problem, settings, presolve_time_limit, 1);
+  return presolve_status_to_mip_status(status);
+}
+
+template <typename i_t, typename f_t>
+void presolver_t<i_t, f_t>::uncrush(const std::vector<f_t>& reduced_primal,
+                                    std::vector<f_t>& full_primal) const
+{
+  third_party_presolver_.uncrush_primal_solution(reduced_primal, full_primal);
+}
+template <typename i_t, typename f_t>
+void presolver_t<i_t, f_t>::crush(const std::vector<f_t>& full_primal,
+                                  std::vector<f_t>& reduced_primal) const
+{
+  third_party_presolver_.crush_primal_solution(full_primal, reduced_primal);
+}
+
+template <typename i_t, typename f_t>
+const std::vector<i_t>& presolver_t<i_t, f_t>::reduced_to_original_map() const
+{
+  return third_party_presolver_.get_reduced_to_original_map();
+}
+
+template <typename i_t, typename f_t>
+const std::vector<i_t>& presolver_t<i_t, f_t>::original_to_reduced_map() const
+{
+  return third_party_presolver_.get_original_to_reduced_map();
+}
+
+#if MIP_INSTANTIATE_FLOAT
+template class presolver_t<int, float>;
+#endif
+
+#if MIP_INSTANTIATE_DOUBLE
+template class presolver_t<int, double>;
+#endif
+
+}  // namespace cuopt::mathematical_optimization::mip

cpp/src/branch_and_bound/presolve.hpp

@@ -0,0 +1,44 @@
+/* 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
+
+#include <vector>
+
+#include <dual_simplex/simplex_solver_settings.hpp>
+#include <dual_simplex/user_problem.hpp>
+#include <mip_heuristics/mip_constants.hpp>
+#include <mip_heuristics/presolve/third_party_presolve.hpp>
+
+#include "constants.hpp"
+
+namespace cuopt::mathematical_optimization::mip {
+
+// Thin owner of a PaPILO presolver scoped to a single sub-MIP solve. apply() reduces the
+// problem in place; the retained column maps let a reduced-space solution be mapped back to
+// the original column space via uncrush().
+template <typename i_t, typename f_t>
+class presolver_t {
+ public:
+  // Presolve `problem` in place using PaPILO. Returns the presolve status; on
+  // INFEASIBLE/UNBOUNDED the problem is left untouched.
+  mip_status_t apply(simplex::user_problem_t<i_t, f_t>& problem,
+                     const simplex::simplex_solver_settings_t<i_t, f_t>& settings);
+
+  // Map a reduced-space primal solution back to the original column space.
+  void uncrush(const std::vector<f_t>& reduced_primal, std::vector<f_t>& full_primal) const;
+
+  void crush(const std::vector<f_t>& full_primal, std::vector<f_t>& reduced_primal) const;
+
+  const std::vector<i_t>& reduced_to_original_map() const;
+  const std::vector<i_t>& original_to_reduced_map() const;
+
+ private:
+  third_party_presolve_t<i_t, f_t> third_party_presolver_;
+};
+
+}  // namespace cuopt::mathematical_optimization::mip