[BP 3] Objective propagation

cpp/src/branch_and_bound/branch_and_bound.cpp

@@ -1213,7 +1213,6 @@ struct deterministic_bfs_policy_t
                          branch_direction_t dir) override
   {
     switch (status) {
-      case node_status_t::INFEASIBLE: this->worker.record_infeasible(node); break;
       case node_status_t::FATHOMED: this->worker.record_fathomed(node, node->lower_bound); break;
       case node_status_t::INTEGER_FEASIBLE:
         this->worker.record_integer_solution(node, node->lower_bound);
@@ -1365,6 +1364,23 @@ void branch_and_bound_t<i_t, f_t>::snap_to_lattice(const lp_problem_t<i_t, f_t>&
   }
 }
 
+template <typename i_t, typename f_t>
+f_t branch_and_bound_t<i_t, f_t>::compute_node_cutoff(const lp_problem_t<i_t, f_t>& lp,
+                                                      f_t upper_bound) const
+{
+  if (lp.objective_step.has_step()) {
+    f_t step = lp.objective_step.step_size;
+    f_t bias = lp.objective_step.bias;
+    // Any improving feasible solution lies on the lattice strictly below upper_bound.
+    f_t k = std::ceil((upper_bound - bias) / step - settings_.integer_tol);
+    return (k - 1) * step + bias;
+  } else if (lp.objective_is_integral) {
+    // Any improving feasible solution has an integer objective strictly below upper_bound.
+    return std::ceil(upper_bound - settings_.integer_tol) - 1;
+  }
+  return upper_bound;
+}
+
 template <typename i_t, typename f_t>
 template <typename WorkerT, typename Policy>
 std::pair<node_status_t, branch_direction_t> branch_and_bound_t<i_t, f_t>::update_tree_impl(
@@ -1384,15 +1400,14 @@ std::pair<node_status_t, branch_direction_t> branch_and_bound_t<i_t, f_t>::updat
   worker->recompute_basis = true;
 
   if (lp_status == dual_status_t::DUAL_UNBOUNDED) {
-    node_ptr->lower_bound = inf;
-    policy.graphviz(search_tree, node_ptr, "infeasible", 0.0);
-    search_tree.update(node_ptr, node_status_t::INFEASIBLE);
-    status = node_status_t::INFEASIBLE;
+    node_ptr->lower_bound = upper_bound;
+    policy.graphviz(search_tree, node_ptr, "infeasible", upper_bound);
+    search_tree.update(node_ptr, node_status_t::FATHOMED);
+    status = node_status_t::FATHOMED;
 
   } else if (lp_status == dual_status_t::CUTOFF) {
-    f_t leaf_obj          = compute_objective(leaf_problem, leaf_solution.x);
     node_ptr->lower_bound = upper_bound;
-    policy.graphviz(search_tree, node_ptr, "cut off", leaf_obj);
+    policy.graphviz(search_tree, node_ptr, "cut off", upper_bound);
     search_tree.update(node_ptr, node_status_t::FATHOMED);
     status = node_status_t::FATHOMED;
 
@@ -1584,21 +1599,8 @@ dual_status_t branch_and_bound_t<i_t, f_t>::solve_node_lp(
     lp_settings.concurrent_halt = &node_concurrent_halt_;
   }
   lp_settings.set_log(false);
-  f_t cutoff = upper_bound_.load();
-  if (worker->leaf_problem.objective_step.has_step()) {
-    f_t step = worker->leaf_problem.objective_step.step_size;
-    f_t bias = worker->leaf_problem.objective_step.bias;
-    // Any improving feasible solution must have objective <= cutoff - step.
-    f_t k               = std::floor((cutoff - bias) / step + settings_.integer_tol);
-    lp_settings.cut_off = (k - 1) * step + bias + settings_.dual_tol;
-  } else if (worker->leaf_problem.objective_is_integral) {
-    // If the objective is integral, any feasible solution should produce an upper bound that is
-    // (approximately) integral. We add a small tolerance and floor this value to get an integer,
-    // we then subtract 1, to stop simplex on problems that cannot improve the primal objective.
-    lp_settings.cut_off = std::floor(cutoff + settings_.integer_tol) - 1 + settings_.dual_tol;
-  } else {
-    lp_settings.cut_off = cutoff + settings_.dual_tol;
-  }
+  const f_t cutoff       = compute_node_cutoff(worker->leaf_problem, upper_bound_.load());
+  lp_settings.cut_off    = cutoff + settings_.dual_tol;
   lp_settings.inside_mip = 2;
   lp_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time);
   if (lp_settings.time_limit <= 0.0) { return dual_status_t::TIME_LIMIT; }
@@ -1642,6 +1644,13 @@ dual_status_t branch_and_bound_t<i_t, f_t>::solve_node_lp(
       lp_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time);
       if (lp_settings.time_limit <= 0.0) { return dual_status_t::TIME_LIMIT; }
 
+      if (worker->search_strategy == search_strategy_t::BEST_FIRST) {
+        if (!worker->domain.propagate_objective(
+              worker->Arow, worker->var_types, settings_, worker->leaf_problem, cutoff)) {
+          return dual_status_t::CUTOFF;
+        }
+      }
+
       i_t node_iter          = 0;
       f_t lp_start_time      = tic();
       f_t node_work_estimate = 0.0;
@@ -1729,7 +1738,6 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(bfs_worker_t<i_t, f_t>* worker,
 
   worker->recompute_basis = true;
   worker->ensure_orbital_fixing();
-
   // The branching of the start node itself is applied when its LP is solved.
   if (start_node->parent != nullptr) {
     start_node->parent->rebuild_variable_bounds(worker->domain);
@@ -4260,7 +4268,6 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
     return mip_status_t::OPTIMAL;
   }
 
-  is_running_ = true;
   lp_settings.concurrent_halt =
     settings_.concurrent_halt ? settings_.concurrent_halt : &node_concurrent_halt_;
   reduced_cost_bounds_t<i_t, f_t> reduced_cost_bounds(original_lp_.num_cols);
@@ -4355,6 +4362,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   }
 
   if (num_fractional != 0 && settings_.max_cut_passes > 0) { print_table_header(); }
+  is_running_ = true;
 
   cut_pool_t<i_t, f_t> cut_pool(original_lp_.num_cols, settings_);
   cut_generation_t<i_t, f_t> cut_generation(cut_pool,
@@ -5269,7 +5277,8 @@ node_status_t branch_and_bound_t<i_t, f_t>::solve_node_deterministic(
   simplex_solver_settings_t<i_t, f_t> lp_settings = settings_;
   lp_settings.set_log(false);
 
-  lp_settings.cut_off       = worker.local_upper_bound + settings_.dual_tol;
+  const f_t cutoff          = compute_node_cutoff(worker.leaf_problem, worker.local_upper_bound);
+  lp_settings.cut_off       = cutoff + settings_.dual_tol;
   lp_settings.inside_mip    = 2;
   lp_settings.time_limit    = remaining_time;
   lp_settings.scale_columns = false;
@@ -5284,19 +5293,32 @@ node_status_t branch_and_bound_t<i_t, f_t>::solve_node_deterministic(
     feasible = worker.rebuild_bounds_from_stack(original_lp_, settings_);
   }
 
+  const bool is_cutoff =
+    feasible && !worker.domain.propagate_objective(
+                  worker.Arow, worker.var_types, settings_, worker.leaf_problem, cutoff);
+
   if (settings_.deterministic) {
     // TEMP APPROXIMATION;
     worker.work_context.record_work_sync_on_horizon(worker.domain.last_nnz_processed / 1e8);
   }
 
   if (!feasible) {
-    node_ptr->lower_bound = std::numeric_limits<f_t>::infinity();
-    search_tree.update(node_ptr, node_status_t::INFEASIBLE);
-    worker.record_infeasible(node_ptr);
+    node_ptr->lower_bound = worker.local_upper_bound;
+    search_tree.update(node_ptr, node_status_t::FATHOMED);
+    worker.record_fathomed(node_ptr, worker.local_upper_bound);
     --exploration_stats_.nodes_unexplored;
     ++exploration_stats_.nodes_explored;
     worker.recompute_basis = true;
-    return node_status_t::INFEASIBLE;
+    return node_status_t::FATHOMED;
+  }
+
+  if (is_cutoff) {
+    ++exploration_stats_.nodes_explored;
+    --exploration_stats_.nodes_unexplored;
+    deterministic_bfs_policy_t<i_t, f_t> policy{*this, worker};
+    auto [status, round_dir] =
+      update_tree_impl(node_ptr, search_tree, &worker, dual_status_t::CUTOFF, policy);
+    return status;
   }
 
   // Solve LP relaxation
@@ -5894,12 +5916,15 @@ void branch_and_bound_t<i_t, f_t>::deterministic_dive(
     // Setup LP settings
     simplex_solver_settings_t<i_t, f_t> lp_settings = settings_;
     lp_settings.set_log(false);
-    lp_settings.cut_off       = worker.local_upper_bound + settings_.dual_tol;
+    const f_t cutoff          = compute_node_cutoff(worker.leaf_problem, worker.local_upper_bound);
+    lp_settings.cut_off       = cutoff + settings_.dual_tol;
     lp_settings.inside_mip    = 2;
     lp_settings.time_limit    = remaining_time;
     lp_settings.scale_columns = false;
 
-    feasible = worker.update_variable_bounds(node_ptr, settings_);
+    feasible = worker.update_variable_bounds(node_ptr, settings_) &&
+               worker.domain.propagate_objective(
+                 worker.Arow, worker.var_types, settings_, worker.leaf_problem, cutoff);
 
     if (settings_.deterministic) {
       // TEMP APPROXIMATION;

cpp/src/branch_and_bound/branch_and_bound.hpp

@@ -356,6 +356,10 @@ class branch_and_bound_t {
                        mip_node_t<i_t, f_t>* node_ptr,
                        f_t leaf_obj);
 
+  // Largest objective value of a solution that improves on upper_bound. This is the value just
+  // below upper_bound on the objective lattice if the objective is integral or moves in steps.
+  f_t compute_node_cutoff(const simplex::lp_problem_t<i_t, f_t>& lp, f_t upper_bound) const;
+
   // Launch a new best-first worker from a given bfs worker.
   void launch_bfs_worker(bfs_worker_t<i_t, f_t>* worker);
 

cpp/src/branch_and_bound/deterministic_workers.hpp

@@ -261,15 +261,6 @@ class deterministic_bfs_worker_t
     ++total_nodes_pruned;
   }
 
-  void record_infeasible(mip_node_t<i_t, f_t>* node)
-  {
-    record_event(
-      bb_event_t<i_t, f_t>::make_infeasible(this->clock, this->worker_id, node->creation_seq));
-    ++nodes_processed_this_horizon;
-    ++this->total_nodes_processed;
-    ++total_nodes_infeasible;
-  }
-
   void record_numerical(mip_node_t<i_t, f_t>* node)
   {
     record_event(

cpp/src/branch_and_bound/mip_node.hpp

@@ -27,16 +27,15 @@ namespace cuopt::mathematical_optimization::mip {
 enum class node_status_t : int {
   PENDING          = 0,  // Node is still in the tree, waiting to be solved
   INTEGER_FEASIBLE = 1,  // Node has an integer feasible solution
-  INFEASIBLE       = 2,  // Node is infeasible
-  FATHOMED         = 3,  // Node objective is greater than the upper bound
-  HAS_CHILDREN     = 4,  // Node has children to explore
-  NUMERICAL        = 5   // Encountered numerical issue when solving the LP relaxation
+  FATHOMED         = 2,  // Node is infeasible or has objective is greater than the upper bound
+  HAS_CHILDREN     = 3,  // Node has children to explore
+  NUMERICAL        = 4   // Encountered numerical issue when solving the LP relaxation
 };
 
 inline bool inactive_status(node_status_t status)
 {
   return (status == node_status_t::FATHOMED || status == node_status_t::INTEGER_FEASIBLE ||
-          status == node_status_t::INFEASIBLE || status == node_status_t::NUMERICAL);
+          status == node_status_t::NUMERICAL);
 }
 
 template <typename i_t, typename f_t>

cpp/src/dual_simplex/domain.cpp

@@ -85,6 +85,12 @@ void domain_t<i_t, f_t>::compute_activities(const csr_matrix_t<i_t, f_t>& Arow,
   for (i_t i = 0; i < Arow.m; ++i) {
     compute_row_activity(i, Arow, lp.lower, lp.upper);
   }
+
+  objective_vars.clear();
+  objective_vars.reserve(lp.num_cols);
+  for (i_t j = 0; j < lp.num_cols; ++j) {
+    if (lp.objective[j] != 0) { objective_vars.push_back(j); }
+  }
 }
 
 template <typename i_t, typename f_t>
@@ -422,6 +428,95 @@ bool domain_t<i_t, f_t>::apply_and_propagate(const csr_matrix_t<i_t, f_t>& Arow,
   return run_bound_propagation(Arow, var_types, settings, lp);
 }
 
+template <typename i_t, typename f_t>
+bool domain_t<i_t, f_t>::propagate_objective(const csr_matrix_t<i_t, f_t>& Arow,
+                                             const std::vector<variable_type_t>& var_types,
+                                             const simplex_solver_settings_t<i_t, f_t>& settings,
+                                             lp_problem_t<i_t, f_t>& lp,
+                                             f_t cutoff)
+{
+  if (!std::isfinite(cutoff) || objective_vars.empty()) { return true; }
+
+  // The work is added to the one of the preceding propagation, so the callers record both.
+  const size_t prior_nnz = last_nnz_processed;
+  size_t work            = 0;
+  const f_t rhs          = cutoff + settings.primal_tol * std::max<f_t>(1.0, std::abs(cutoff));
+  bool feasible          = true;
+
+  for (i_t round = 0; round < params.max_objective_rounds; ++round) {
+    // Minimum objective over the bounds in lp, with the infinite contributions flagged as in
+    // compute_row_activity.
+    f_t min     = 0;
+    f_t min_err = 0;
+    i_t min_inf = 0;
+    for (i_t j : objective_vars) {
+      const f_t c_j       = lp.objective[j];
+      const f_t bound_min = c_j < 0 ? lp.upper[j] : lp.lower[j];
+      const f_t alpha_min = c_j * bound_min;
+      if (std::isfinite(alpha_min) && std::abs(alpha_min) < params.huge_value) {
+        dot2_add(c_j, bound_min, min, min_err);
+      } else {
+        ++min_inf;
+      }
+    }
+    work += objective_vars.size();
+
+    if (min_inf == 0 && min + min_err > rhs) {
+      feasible = false;
+      break;
+    }
+    if (min_inf > 1) { break; }
+
+    const size_t num_changes = bound_changes.size();
+    for (i_t j : objective_vars) {
+      const f_t lb = lp.lower[j];
+      const f_t ub = lp.upper[j];
+      if (lb == ub) { continue; }
+
+      const f_t c_j       = lp.objective[j];
+      const f_t bound_min = c_j < 0 ? ub : lb;
+      const f_t alpha_min = c_j * bound_min;
+      const bool alpha_min_inf =
+        !std::isfinite(alpha_min) || std::abs(alpha_min) >= params.huge_value;
+
+      // With a single infinite contribution, only its variable can be bounded.
+      if (min_inf == 1 && !alpha_min_inf) { continue; }
+
+      // c_j x_j <= rhs - (minimum objective of the other variables): an upper bound if c_j > 0, a
+      // lower bound otherwise. The slack is formed with Dot2, as in the row propagation.
+      f_t slack     = rhs;
+      f_t slack_err = -min_err;
+      dot2_add(-1.0, min, slack, slack_err);
+      if (min_inf == 0) { dot2_add(c_j, bound_min, slack, slack_err); }
+      const f_t gamma = (slack + slack_err) / c_j;
+      if (std::abs(gamma) >= params.max_derived_bound) { continue; }
+
+      const f_t new_lower = c_j > 0 ? lb : gamma;
+      const f_t new_upper = c_j > 0 ? gamma : ub;
+      if (!tighten_bounds(
+            j, new_lower, new_upper, bound_change_origin_t::OBJECTIVE, var_types, settings, lp)) {
+        feasible = false;
+        break;
+      }
+    }
+    work += objective_vars.size();
+
+    if (!feasible || bound_changes.size() == num_changes) { break; }
+
+    // Propagate the rows of the tightened variables, which may raise the minimum objective.
+    feasible = run_bound_propagation(Arow, var_types, settings, lp);
+    work += last_nnz_processed;
+    if (!feasible) { break; }
+  }
+
+  if (!feasible) { clear_queues(); }
+
+  work += nnz_processed;
+  nnz_processed      = 0;
+  last_nnz_processed = prior_nnz + work;
+  return feasible;
+}
+
 template <typename i_t, typename f_t>
 void domain_t<i_t, f_t>::backtrack_to_parent(lp_problem_t<i_t, f_t>& lp)
 {

cpp/src/dual_simplex/domain.hpp

@@ -59,6 +59,8 @@ struct domain_params {
   double min_improvement_factor   = 1e3;
   // Derived bounds larger than this in magnitude are discarded
   double max_derived_bound = 1e8;
+  // Maximum number of rounds of objective propagation followed by row propagation
+  int max_objective_rounds = 3;
 };
 
 // The finite parts of the activities are accumulated with compensated (Dot2) summation: the
@@ -164,7 +166,18 @@ class domain_t {
                            bound_change_t<i_t, f_t> bound_change,
                            lp_problem_t<i_t, f_t>& lp);
 
+  // Propagates c^T x <= cutoff, followed by the rows affected by its tightenings, for up to
+  // max_objective_rounds rounds. Returns false if no solution within the bounds in lp can reach
+  // the cutoff.
+  bool propagate_objective(const csr_matrix_t<i_t, f_t>& Arow,
+                           const std::vector<variable_type_t>& var_types,
+                           const simplex_solver_settings_t<i_t, f_t>& settings,
+                           lp_problem_t<i_t, f_t>& lp,
+                           f_t cutoff);
+
   size_t last_nnz_processed{0};
+  size_t objective_fixings{0};
+  size_t objective_cutoffs{0};
 
   // Implications applied when an integer variable is fixed to 0 or 1. Either may be null, and the
   // clique table is skipped until it is ready.
@@ -181,6 +194,9 @@ class domain_t {
   std::vector<uint8_t> var_queued;
   circular_deque_t<i_t> var_queue;
 
+  // Variables with a nonzero objective coefficient.
+  std::vector<i_t> objective_vars;
+
   std::vector<bound_change_t<i_t, f_t>> bound_changes;
 
   size_t nnz_processed{0};