MIP symmetry detection and orbital fixing using dejavu

cpp/CMakeLists.txt

@@ -276,6 +276,20 @@ FetchContent_MakeAvailable(pslp)
 set(BUILD_SHARED_LIBS ${BUILD_SHARED_LIBS_SAVED})
 
 
+# dejavu - header-only graph automorphism library for MIP symmetry detection
+# https://github.com/markusa4/dejavu (header-only, skip its CMakeLists.txt)
+FetchContent_Declare(
+  dejavu
+  GIT_REPOSITORY "https://github.com/markusa4/dejavu.git"
+  GIT_TAG "v2.1"
+  GIT_PROGRESS TRUE
+  EXCLUDE_FROM_ALL
+  SYSTEM
+  SOURCE_SUBDIR _nonexistent
+)
+FetchContent_MakeAvailable(dejavu)
+message(STATUS "dejavu (graph automorphism): ${dejavu_SOURCE_DIR}")
+
 include(${rapids-cmake-dir}/cpm/rapids_logger.cmake)
 # generate logging macros
 rapids_cpm_rapids_logger(BUILD_EXPORT_SET cuopt-exports INSTALL_EXPORT_SET cuopt-exports)
@@ -471,7 +485,8 @@ target_include_directories(cuopt PRIVATE
 )
 
 target_include_directories(cuopt SYSTEM PRIVATE
-        "${pslp_SOURCE_DIR}/include"
+  "${pslp_SOURCE_DIR}/include"
+  "${dejavu_SOURCE_DIR}"
 )
 
 target_include_directories(cuopt

cpp/include/cuopt/linear_programming/constants.h

@@ -59,6 +59,7 @@
 #define CUOPT_MIP_HEURISTICS_ONLY                  "mip_heuristics_only"
 #define CUOPT_MIP_SCALING                          "mip_scaling"
 #define CUOPT_MIP_PRESOLVE                         "mip_presolve"
+#define CUOPT_MIP_SYMMETRY                         "mip_symmetry"
 #define CUOPT_MIP_RELIABILITY_BRANCHING            "mip_reliability_branching"
 #define CUOPT_MIP_CUT_PASSES                       "mip_cut_passes"
 #define CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS      "mip_mixed_integer_rounding_cuts"

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

@@ -90,6 +90,7 @@ class mip_solver_settings_t {
   bool heuristics_only          = false;
   i_t reliability_branching     = -1;
   i_t num_cpu_threads           = -1;  // -1 means use default number of threads in branch and bound
+  i_t symmetry                  = -1;
   i_t max_cut_passes            = 10;  // number of cut passes to make
   i_t mir_cuts                  = -1;
   i_t mixed_integer_gomory_cuts = -1;

cpp/src/branch_and_bound/branch_and_bound.cpp

@@ -8,6 +8,7 @@
 #include <branch_and_bound/branch_and_bound.hpp>
 #include <branch_and_bound/mip_node.hpp>
 #include <branch_and_bound/pseudo_costs.hpp>
+#include <branch_and_bound/symmetry.hpp>
 
 #include <cuts/cuts.hpp>
 #include <mip_heuristics/presolve/conflict_graph/clique_table.cuh>
@@ -248,11 +249,13 @@ branch_and_bound_t<i_t, f_t>::branch_and_bound_t(
   const simplex_solver_settings_t<i_t, f_t>& solver_settings,
   f_t start_time,
   const probing_implied_bound_t<i_t, f_t>& probing_implied_bound,
-  std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table)
+  std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table,
+  mip_symmetry_t<i_t, f_t>* symmetry)
   : original_problem_(user_problem),
     settings_(solver_settings),
     probing_implied_bound_(probing_implied_bound),
     clique_table_(std::move(clique_table)),
+    symmetry_(symmetry),
     original_lp_(user_problem.handle_ptr, 1, 1, 1),
     Arow_(1, 1, 0),
     incumbent_(1),
@@ -730,6 +733,20 @@ void branch_and_bound_t<i_t, f_t>::set_final_solution(mip_solution_t<i_t, f_t>&
   settings_.log.printf("Explored %d nodes in %.2fs.\n",
                        exploration_stats_.nodes_explored,
                        toc(exploration_stats_.start_time));
+  if (exploration_stats_.orbital_fixing_nodes.load() > 0 ||
+      exploration_stats_.orbital_conflict_nodes.load() > 0) {
+    settings_.log.printf("Orbital fixing applied at %lld nodes, %lld total variable fixings, "
+                         "%lld nodes with conflicting orbits\n",
+                         (long long)exploration_stats_.orbital_fixing_nodes.load(),
+                         (long long)exploration_stats_.orbital_fixings_applied.load(),
+                         (long long)exploration_stats_.orbital_conflict_nodes.load());
+  }
+  if (exploration_stats_.lexical_reduction_nodes.load() > 0) {
+    settings_.log.printf("Lexical reduction applied at %lld nodes, %lld total variable fixings, %lld nodes pruned\n",
+                         (long long)exploration_stats_.lexical_reduction_nodes.load(),
+                         (long long)exploration_stats_.lexical_reduction_fixings_applied.load(),
+                         (long long)exploration_stats_.lexical_reduction_pruned_nodes.load());
+  }
   settings_.log.printf("Absolute Gap %e Objective %.16e %s Bound %.16e\n",
                        gap,
                        obj,
@@ -1387,42 +1404,79 @@ dual::status_t branch_and_bound_t<i_t, f_t>::solve_node_lp(
   bool feasible            = worker->set_lp_variable_bounds(node_ptr, settings_);
   dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED;
   worker->leaf_edge_norms  = edge_norms_;
+  if (worker->recompute_bounds && worker->orbital_fixing &&
+      worker->search_strategy == BEST_FIRST) {
+    worker->orbital_fixing->reset(symmetry_, node_ptr);
+  }
 
   if (feasible) {
-    i_t node_iter     = 0;
-    f_t lp_start_time = tic();
-
-    lp_status = dual_phase2_with_advanced_basis(2,
-                                                0,
-                                                worker->recompute_basis,
-                                                lp_start_time,
-                                                worker->leaf_problem,
-                                                lp_settings,
-                                                leaf_vstatus,
-                                                worker->basis_factors,
-                                                worker->basic_list,
-                                                worker->nonbasic_list,
-                                                worker->leaf_solution,
-                                                node_iter,
-                                                worker->leaf_edge_norms);
 
-    if (lp_status == dual::status_t::NUMERICAL) {
-      log.debug("Numerical issue node %d. Resolving from scratch.\n", node_ptr->node_id);
-      lp_status_t second_status = solve_linear_program_with_advanced_basis(worker->leaf_problem,
-                                                                           lp_start_time,
-                                                                           lp_settings,
-                                                                           worker->leaf_solution,
-                                                                           worker->basis_factors,
-                                                                           worker->basic_list,
-                                                                           worker->nonbasic_list,
-                                                                           leaf_vstatus,
-                                                                           worker->leaf_edge_norms);
+    // Perform orbital fixing
+    auto* of = worker->orbital_fixing.get();
+    if (of != nullptr && !of->disabled()) {
+      i_t prev_fix = node_ptr->orbital_fix_zero.size() + node_ptr->orbital_fix_one.size();
+      i_t conflicts = of->orbital_fixing(symmetry_, settings_, node_ptr, worker->leaf_problem,
+                                         worker->start_lower, worker->start_upper);
+      i_t new_fix = node_ptr->orbital_fix_zero.size() + node_ptr->orbital_fix_one.size();
+      if (new_fix > prev_fix) {
+        ++stats.orbital_fixing_nodes;
+        stats.orbital_fixings_applied += (new_fix - prev_fix);
+      }
+      if (conflicts > 0) { ++stats.orbital_conflict_nodes; }
+    } else if (of != nullptr) {
+      of->propagate_cumulative_fixings(node_ptr);
+    }
 
-      lp_status = convert_lp_status_to_dual_status(second_status);
+    if (settings_.symmetry == 2 && worker->lexical_reduction != nullptr) {
+      i_t lexical_reductions_info =
+        worker->lexical_reduction->lexical_reduce(symmetry_, node_ptr, worker->leaf_problem);
+      if (lexical_reductions_info > 0) {
+        stats.lexical_reduction_nodes++;
+        stats.lexical_reduction_fixings_applied += lexical_reductions_info;
+      }
+      if (lexical_reductions_info == -1) {
+        feasible = false;
+        stats.lexical_reduction_pruned_nodes++;
+      }
     }
 
-    stats.total_lp_solve_time += toc(lp_start_time);
-    stats.total_lp_iters += node_iter;
+    if (feasible) {
+      i_t node_iter     = 0;
+      f_t lp_start_time = tic();
+
+      lp_status = dual_phase2_with_advanced_basis(2,
+                                                  0,
+                                                  worker->recompute_basis,
+                                                  lp_start_time,
+                                                  worker->leaf_problem,
+                                                  lp_settings,
+                                                  leaf_vstatus,
+                                                  worker->basis_factors,
+                                                  worker->basic_list,
+                                                  worker->nonbasic_list,
+                                                  worker->leaf_solution,
+                                                  node_iter,
+                                                  worker->leaf_edge_norms);
+
+      if (lp_status == dual::status_t::NUMERICAL) {
+        log.debug("Numerical issue node %d. Resolving from scratch.\n", node_ptr->node_id);
+        lp_status_t second_status =
+          solve_linear_program_with_advanced_basis(worker->leaf_problem,
+                                                   lp_start_time,
+                                                   lp_settings,
+                                                   worker->leaf_solution,
+                                                   worker->basis_factors,
+                                                   worker->basic_list,
+                                                   worker->nonbasic_list,
+                                                   leaf_vstatus,
+                                                   worker->leaf_edge_norms);
+
+        lp_status = convert_lp_status_to_dual_status(second_status);
+      }
+
+      stats.total_lp_solve_time += toc(lp_start_time);
+      stats.total_lp_iters += node_iter;
+    }
   }
 
 #ifdef LOG_NODE_SIMPLEX
@@ -1438,6 +1492,7 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(branch_and_bound_worker_t<i_t, f_
   stack.push_front(worker->start_node);
   worker->recompute_basis  = true;
   worker->recompute_bounds = true;
+  worker->ensure_orbital_fixing();
 
   f_t lower_bound = get_lower_bound();
   f_t upper_bound = upper_bound_;
@@ -1563,6 +1618,7 @@ template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::dive_with(branch_and_bound_worker_t<i_t, f_t>* worker)
 {
   raft::common::nvtx::range scope("BB::diving_thread");
+  if (worker->orbital_fixing) { worker->orbital_fixing->disable(); }
   logger_t log;
   log.log = false;
 
@@ -1659,7 +1715,7 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
   std::array<i_t, num_search_strategies> max_num_workers_per_type =
     get_max_workers(num_workers, strategies);
 
-  worker_pool_.init(num_workers, original_lp_, Arow_, var_types_, settings_);
+  worker_pool_.init(num_workers, original_lp_, Arow_, var_types_, symmetry_, settings_);
   active_workers_per_strategy_.fill(0);
 
 #ifdef CUOPT_LOG_DEBUG
@@ -1801,7 +1857,7 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
 template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::single_threaded_solve()
 {
-  worker_pool_.init(1, original_lp_, Arow_, var_types_, settings_);
+  worker_pool_.init(1, original_lp_, Arow_, var_types_, symmetry_, settings_);
   branch_and_bound_worker_t<i_t, f_t>* worker = worker_pool_.get_idle_worker();
 
   f_t lower_bound = get_lower_bound();
@@ -2524,6 +2580,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
                                basic_list,
                                nonbasic_list,
                                basis_update,
+                               symmetry_,
                                pc_);
   }
 
@@ -2593,6 +2650,17 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   node_queue_.push(search_tree_.root.get_down_child());
   node_queue_.push(search_tree_.root.get_up_child());
 
+  if (symmetry_ != nullptr) {
+    i_t removed = symmetry_->generators.template prune_by_bounds<f_t>(
+      original_lp_.lower, original_lp_.upper);
+    if (removed > 0) {
+      symmetry_->num_generators = static_cast<int>(symmetry_->generators.num_generators());
+      settings_.log.printf(
+        "Pruned %d generators invalidated by root-level bound tightening, %d remain\n",
+        removed, symmetry_->num_generators);
+    }
+  }
+
   settings_.log.printf("Exploring the B&B tree using %d threads\n\n", settings_.num_threads);
   node_concurrent_halt_ = 0;
 

cpp/src/branch_and_bound/branch_and_bound.hpp

@@ -47,6 +47,9 @@ struct clique_table_t;
 
 namespace cuopt::linear_programming::dual_simplex {
 
+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
@@ -80,7 +83,8 @@ class branch_and_bound_t {
                      const simplex_solver_settings_t<i_t, f_t>& solver_settings,
                      f_t start_time,
                      const probing_implied_bound_t<i_t, f_t>& probing_implied_bound,
-                     std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table = nullptr);
+                     std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table = nullptr,
+                     mip_symmetry_t<i_t, f_t>* symmetry = nullptr);
 
   // 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; }
@@ -162,6 +166,7 @@ class branch_and_bound_t {
   const simplex_solver_settings_t<i_t, f_t> settings_;
   const probing_implied_bound_t<i_t, f_t>& probing_implied_bound_;
   std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table_;
+  mip_symmetry_t<i_t, f_t>* symmetry_;
   std::future<std::shared_ptr<detail::clique_table_t<i_t, f_t>>> clique_table_future_;
   std::atomic<bool> signal_extend_cliques_{false};
 

cpp/src/branch_and_bound/branch_and_bound_worker.hpp

@@ -8,6 +8,7 @@
 #pragma once
 
 #include <branch_and_bound/mip_node.hpp>
+#include <branch_and_bound/symmetry.hpp>
 
 #include <dual_simplex/basis_updates.hpp>
 #include <dual_simplex/bounds_strengthening.hpp>
@@ -46,6 +47,12 @@ struct branch_and_bound_stats_t {
   omp_atomic_t<int64_t> total_lp_iters   = 0;
   omp_atomic_t<i_t> nodes_since_last_log = 0;
   omp_atomic_t<f_t> last_log             = 0.0;
+  omp_atomic_t<int64_t> orbital_fixing_nodes = 0;
+  omp_atomic_t<int64_t> orbital_fixings_applied = 0;
+  omp_atomic_t<int64_t> orbital_conflict_nodes = 0;
+  omp_atomic_t<int64_t> lexical_reduction_nodes = 0;
+  omp_atomic_t<int64_t> lexical_reduction_fixings_applied = 0;
+  omp_atomic_t<int64_t> lexical_reduction_pruned_nodes = 0;
 };
 
 template <typename i_t, typename f_t>
@@ -73,6 +80,20 @@ class branch_and_bound_worker_t {
 
   pcgenerator_t rng;
 
+  std::unique_ptr<orbital_fixing_t<i_t, f_t>> orbital_fixing;
+  std::unique_ptr<lexical_reduction_t<i_t, f_t>> lexical_reduction;
+  mip_symmetry_t<i_t, f_t>* symmetry_ptr = nullptr;
+
+  void ensure_orbital_fixing()
+  {
+    if (orbital_fixing == nullptr && symmetry_ptr != nullptr) {
+      orbital_fixing = std::make_unique<orbital_fixing_t<i_t, f_t>>(*symmetry_ptr);
+    }
+    if (lexical_reduction == nullptr && symmetry_ptr != nullptr) {
+      lexical_reduction = std::make_unique<lexical_reduction_t<i_t, f_t>>(symmetry_ptr->num_original_vars);
+    }
+  }
+
   bool recompute_basis  = true;
   bool recompute_bounds = true;
 
@@ -169,13 +190,15 @@ class branch_and_bound_worker_pool_t {
             const lp_problem_t<i_t, f_t>& original_lp,
             const csr_matrix_t<i_t, f_t>& Arow,
             const std::vector<variable_type_t>& var_type,
+            mip_symmetry_t<i_t, f_t>* symmetry,
             const simplex_solver_settings_t<i_t, f_t>& settings)
   {
     workers_.resize(num_workers);
     num_idle_workers_ = num_workers;
     for (i_t i = 0; i < num_workers; ++i) {
       workers_[i] = std::make_unique<branch_and_bound_worker_t<i_t, f_t>>(
         i, original_lp, Arow, var_type, settings);
+      workers_[i]->symmetry_ptr = symmetry;
       idle_workers_.push_front(i);
     }
 

cpp/src/branch_and_bound/mip_node.hpp

@@ -271,6 +271,8 @@ class mip_node_t {
     copy.children[1]        = nullptr;
     copy.status             = node_status_t::PENDING;
 
+    copy.orbital_fix_zero = orbital_fix_zero;
+    copy.orbital_fix_one  = orbital_fix_one;
     copy.origin_worker_id = origin_worker_id;
     copy.creation_seq     = creation_seq;
     return copy;
@@ -293,6 +295,12 @@ class mip_node_t {
 
   std::vector<variable_status_t> vstatus;
 
+  // Cumulative orbital fixing bound changes from root to this node.
+  // Stored so that when a child starts a new plunge, the parent's
+  // orbital fixings can be restored without re-derivation.
+  std::vector<i_t> orbital_fix_zero;
+  std::vector<i_t> orbital_fix_one;
+
   // Worker-local identification for deterministic ordering:
   // - origin_worker_id: which worker created this node
   // - creation_seq: sequence number within that worker (cumulative across horizons, serial)