Root heuristics

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

@@ -41,6 +41,11 @@ struct mip_submip_hyper_params_t {
   // number of simplex iteration from the parent B&B.
   f_t iteration_limit_ratio = 0.8;
 
+  // If there is not enough variables fixed or we already found an improving solution,
+  // perform a short DFS to quickly find a feasible solution. This setting controls
+  // the maximum number of nodes allow for backtracking.
+  i_t dfs_max_backtrack = 5;
+
   // Run CPU FJ over the sub-MIP
   bool enable_cpufj = true;
 };

cpp/src/branch_and_bound/branch_and_bound.hpp

@@ -34,12 +34,14 @@
 #include <cuopt/mathematical_optimization/pdlp/solver_settings.hpp>
 
 #include <mip_heuristics/presolve/third_party_presolve.hpp>
+#include <mip_heuristics/root_heuristics.hpp>
 
 #include <omp.h>
 
 #include <atomic>
 #include <functional>
 #include <future>
+#include <list>
 #include <memory>
 #include <vector>
 
@@ -368,24 +370,32 @@ class branch_and_bound_t {
 
   // Launch a new RINS worker
   bool launch_rins_worker(const std::vector<f_t>& sol);
-  void set_solution_from_submip(const std::vector<f_t>& solution,
+  void set_solution_from_submip(const simplex::lp_problem_t<i_t, f_t>& lp,
+                                const std::vector<f_t>& solution,
                                 const third_party_presolve_t<i_t, f_t>& presolver,
-                                f_t fixrate,
-                                f_t obj);
+                                f_t fixrate);
 
   // Solve the RINS sub-MIP
   void solve_submip(diving_worker_t<i_t, f_t>* worker,
                     const std::vector<f_t>& current_incumbent,
+                    const std::vector<simplex::variable_type_t>& var_type,
                     i_t num_var_fixed,
                     i_t num_integers,
                     i_t submip_level,
-                    std::string_view log_prefix);
+                    std::string_view log_prefix,
+                    bool is_root_heuristic);
 
   // Creates and solves the RINS sub-MIP
-  void rins(diving_worker_t<i_t, f_t>* rins_worker, const std::vector<f_t>& node_solution);
-
-  // Get the simplex settings for solving the LP of a single node
-  simplex::simplex_solver_settings_t<i_t, f_t> get_node_lp_settings();
+  void rins(diving_worker_t<i_t, f_t>* worker,
+            const std::vector<f_t>& current_incumbent,
+            const std::vector<simplex::variable_type_t>& var_types,
+            bool is_root_heuristic);
+
+  void launch_root_heuristics(const simplex::lp_problem_t<i_t, f_t>& lp,
+                              const std::vector<f_t>& sol,
+                              i_t cut_pass,
+                              std::list<root_heuristics_t<i_t, f_t>>& heuristics,
+                              omp_atomic_t<i_t>& worker_count);
 
   // Solve the LP relaxation of a leaf node
   simplex::dual_status_t solve_node_lp(mip_node_t<i_t, f_t>* node_ptr,

cpp/src/branch_and_bound/constants.hpp

@@ -27,15 +27,31 @@ enum class heuristics_origin_t {
 // doi: 10.1007/s10107-004-0518-7.
 enum class search_strategy_t : int {
   BEST_FIRST           = 0,  // Best-First + Plunging.
-  PSEUDOCOST_DIVING    = 1,  // Pseudocost diving (9.2.5)
-  LINE_SEARCH_DIVING   = 2,  // Line search diving (9.2.4)
-  GUIDED_DIVING        = 3,  // Guided diving (9.2.3).
-  COEFFICIENT_DIVING   = 4,  // Coefficient diving (9.2.1)
+  PSEUDOCOST_DIVING    = 1,  // Pseudocost diving [1, Section 9.2.5]
+  LINE_SEARCH_DIVING   = 2,  // Line search diving [1, Section 9.2.4]
+  GUIDED_DIVING        = 3,  // Guided diving. [1, Section 9.2.3]
+  COEFFICIENT_DIVING   = 4,  // Coefficient diving [1, Section 9.2.1]
   FARKAS_DIVING        = 5,  // Farkas Diving (see [2])
-  VECTOR_LENGTH_DIVING = 6,  // Vector Length Diving (9.2.6)
-  SUBMIP               = 7   // RINS (see [3])
+  VECTOR_LENGTH_DIVING = 6,  // Vector Length Diving [1, Section 9.2.6]
+  RINS                 = 7,  // RINS (see [3])
 };
 
 enum class branch_direction_t { NONE = -1, DOWN = 0, UP = 1 };
 
+inline const char* search_strategy_to_string(search_strategy_t search_strategy)
+{
+  switch (search_strategy) {
+    case search_strategy_t::BEST_FIRST: return "BEST_FIRST";
+    case search_strategy_t::PSEUDOCOST_DIVING: return "PSEUDOCOST_DIVING";
+    case search_strategy_t::LINE_SEARCH_DIVING: return "LINE_SEARCH_DIVING";
+    case search_strategy_t::GUIDED_DIVING: return "GUIDED_DIVING";
+    case search_strategy_t::COEFFICIENT_DIVING: return "COEFFICIENT_DIVING";
+    case search_strategy_t::FARKAS_DIVING: return "FARKAS_DIVING";
+    case search_strategy_t::VECTOR_LENGTH_DIVING: return "VECTOR_LENGTH_DIVING";
+    case search_strategy_t::RINS: return "RINS";
+  }
+
+  return "UNKNOWN";
+}
+
 }  // namespace cuopt::mathematical_optimization::mip

cpp/src/branch_and_bound/worker.hpp

@@ -243,6 +243,8 @@ class diving_worker_t : public branch_and_bound_worker_t<i_t, f_t> {
   // The best-first worker that is associated with this diving worker. Used for controlling the
   // number of active diving workers.
   bfs_worker_t<i_t, f_t>* bfs_worker{nullptr};
+
+  std::atomic<int> halt = false;
 };
 
 struct submip_stats_t {

cpp/src/dual_simplex/solve.cpp

@@ -104,6 +104,12 @@ f_t compute_user_objective(const lp_problem_t<i_t, f_t>& lp, f_t obj)
   return user_obj;
 }
 
+template <typename i_t, typename f_t>
+f_t compute_internal_objective(const lp_problem_t<i_t, f_t>& lp, f_t user_obj)
+{
+  return user_obj / lp.obj_scale - lp.obj_constant;
+}
+
 template <typename i_t, typename f_t>
 lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original_lp,
                                           const f_t start_time,
@@ -813,6 +819,8 @@ template double compute_user_objective<int, double>(const lp_problem_t<int, doub
 
 template double compute_user_objective(const lp_problem_t<int, double>& lp, double obj);
 
+template double compute_internal_objective(const lp_problem_t<int, double>& lp, double user_obj);
+
 template lp_status_t solve_linear_program_advanced(
   const lp_problem_t<int, double>& original_lp,
   const double start_time,

cpp/src/dual_simplex/solve.hpp

@@ -63,6 +63,9 @@ f_t compute_user_objective(const lp_problem_t<i_t, f_t>& lp, const std::vector<f
 template <typename i_t, typename f_t>
 f_t compute_user_objective(const lp_problem_t<i_t, f_t>& lp, f_t obj);
 
+template <typename i_t, typename f_t>
+f_t compute_internal_objective(const lp_problem_t<i_t, f_t>& lp, f_t user_obj);
+
 template <typename i_t, typename f_t>
 lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original_lp,
                                           const f_t start_time,

cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu

@@ -1814,13 +1814,20 @@ void fj_cpu_worker_t<i_t, f_t>::create_worker(
 }
 
 template <typename i_t, typename f_t>
-void fj_cpu_worker_t<i_t, f_t>::run_async(f_t time_limit, double work_unit_limit)
+void fj_cpu_worker_t<i_t, f_t>::run_async(f_t time_limit,
+                                          double work_unit_limit,
+                                          omp_atomic_t<i_t>* worker_count)
 {
   if (!fj_cpu) return;
 
-#pragma omp task shared(fj_cpu) firstprivate(time_limit, work_unit_limit) \
+  if (worker_count) ++(*worker_count);
+#pragma omp task shared(fj_cpu) firstprivate(time_limit, work_unit_limit, worker_count) \
   priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) depend(out : *fj_cpu)
-  cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit);
+  {
+    cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit);
+    fj_cpu.reset();
+    if (worker_count) --(*worker_count);
+  }
 }
 
 template <typename i_t, typename f_t>

cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh

@@ -10,6 +10,8 @@
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
 
+#include <utilities/omp_helpers.hpp>
+
 #include <atomic>
 #include <functional>
 #include <limits>
@@ -47,8 +49,9 @@ struct fj_cpu_worker_t {
 
   // Run the worker asynchronously (i.e., launch an openmp task and then continue the
   // execution). Call `stop()` for stopping the worker
-  void run_async(f_t time_limit         = std::numeric_limits<f_t>::infinity(),
-                 double work_unit_limit = std::numeric_limits<double>::infinity());
+  void run_async(f_t time_limit                  = std::numeric_limits<f_t>::infinity(),
+                 double work_unit_limit          = std::numeric_limits<double>::infinity(),
+                 omp_atomic_t<i_t>* worker_count = nullptr);
 
   // Run the CPU FJ synchronously (i.e., wait for it to finish before proceeding)
   void run_sync(f_t time_limit         = std::numeric_limits<f_t>::infinity(),

cpp/src/mip_heuristics/presolve/third_party_presolve.cpp

@@ -1290,9 +1290,9 @@ void third_party_presolve_t<i_t, f_t>::undo(std::vector<f_t>& primal_solution,
 
 template <typename i_t, typename f_t>
 void third_party_presolve_t<i_t, f_t>::uncrush_primal_solution(
-  const std::vector<f_t>& reduced_primal, std::vector<f_t>& full_primal) const
+  const std::vector<f_t>& reduced_primal, std::vector<f_t>& full_primal, bool check_postsolve) const
 {
-  if (presolver_ == cuopt::mathematical_optimization::presolver_t::PSLP) {
+  if (presolver_ == PSLP) {
     cuopt_expects(false,
                   error_type_t::RuntimeError,
                   "This code path should be never called, as this is meant for callbacks and they "
@@ -1308,7 +1308,7 @@ void third_party_presolve_t<i_t, f_t>::uncrush_primal_solution(
 
   bool is_optimal = false;
   auto status = post_solver.undo(reduced_sol, full_sol, *papilo_post_solve_storage_, is_optimal);
-  check_postsolve_status(status);
+  if (check_postsolve) check_postsolve_status(status);
   full_primal = std::move(full_sol.primal);
 }
 

cpp/src/mip_heuristics/presolve/third_party_presolve.hpp

@@ -141,7 +141,8 @@ class third_party_presolve_t {
             bool dual_postsolve);
 
   void uncrush_primal_solution(const std::vector<f_t>& reduced_primal,
-                               std::vector<f_t>& full_primal) const;
+                               std::vector<f_t>& full_primal,
+                               bool check_postsolve = true) const;
 
   void crush_primal_solution(const optimization_problem_t<i_t, f_t>& reduced_problem,
                              const std::vector<f_t>& original_primal,

cpp/src/mip_heuristics/root_heuristics.hpp

@@ -0,0 +1,60 @@
+/* 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 <branch_and_bound/worker.hpp>
+#include <dual_simplex/user_problem.hpp>
+#include "feasibility_jump/fj_cpu_worker.cuh"
+
+namespace cuopt::mathematical_optimization::mip {
+
+template <typename i_t, typename f_t>
+struct root_heuristics_t {
+  std::unique_ptr<diving_worker_t<i_t, f_t>> submip_worker_;
+  fj_cpu_worker_t<i_t, f_t> fj_cpu_worker_;
+  std::vector<simplex::variable_type_t> var_types_;
+  csr_matrix_t<i_t, f_t> Arow_;
+
+  root_heuristics_t(const csr_matrix_t<i_t, f_t>& Arow,
+                    const std::vector<simplex::variable_type_t>& var_types)
+    : submip_worker_(nullptr), var_types_(var_types), Arow_(Arow) {};
+
+  ~root_heuristics_t() { stop(); }
+
+  void stop()
+  {
+    fj_cpu_worker_.stop();
+    if (submip_worker_) {
+      submip_worker_->halt              = true;
+      diving_worker_t<i_t, f_t>* worker = submip_worker_.get();
+#pragma omp taskwait depend(in : *worker)
+    }
+  }
+
+  diving_worker_t<i_t, f_t>* create_submip_worker(
+    i_t id,
+    const simplex::lp_problem_t<i_t, f_t>& lp,
+    const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
+    f_t root_obj,
+    const std::vector<simplex::variable_status_t>& root_vstatus,
+    const std::vector<f_t>& sol)
+  {
+    submip_worker_ =
+      std::make_unique<diving_worker_t<i_t, f_t>>(id, lp, Arow_, var_types_, settings);
+    submip_worker_->start_node       = mip_node_t<i_t, f_t>(root_obj, root_vstatus);
+    submip_worker_->leaf_vstatus     = root_vstatus;
+    submip_worker_->leaf_solution.x  = sol;
+    submip_worker_->recompute_bounds = false;
+    submip_worker_->recompute_basis  = true;
+    submip_worker_->search_strategy  = search_strategy_t::RINS;
+    submip_worker_->set_active();
+
+    return submip_worker_.get();
+  }
+};
+}  // namespace cuopt::mathematical_optimization::mip