More numerically robust CPUFJ feasibility checking

cpp/src/mip/feasibility_jump/feasibility_jump_impl_common.cuh

@@ -19,6 +19,24 @@
 
 namespace cuopt::linear_programming::detail {
 
+template <typename i_t, typename f_t, typename Iterator>
+HDI f_t fj_kahan_babushka_neumaier_sum(Iterator begin, Iterator end)
+{
+  f_t sum = 0;
+  f_t c   = 0;
+  for (Iterator it = begin; it != end; ++it) {
+    f_t delta = *it;
+    f_t t     = sum + delta;
+    if (fabs(sum) > fabs(delta)) {
+      c += (sum - t) + delta;
+    } else {
+      c += (delta - t) + sum;
+    }
+    sum = t;
+  }
+  return sum + c;
+}
+
 // Returns the current slack, and the variable delta that would nullify this slack ("tighten" it)
 template <typename i_t, typename f_t>
 HDI thrust::tuple<f_t, f_t> get_mtm_for_bound(

cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu

@@ -35,24 +35,6 @@ namespace cg = cooperative_groups;
 
 namespace cuopt::linear_programming::detail {
 
-template <typename i_t, typename f_t, typename Iterator>
-static DI f_t KahanBabushkaNeumaierSum(Iterator begin, Iterator end)
-{
-  f_t sum = 0;
-  f_t c   = 0;
-  for (Iterator it = begin; it != end; ++it) {
-    f_t delta = *it;
-    f_t t     = sum + delta;
-    if (fabs(sum) > fabs(delta)) {
-      c += (sum - t) + delta;
-    } else {
-      c += (delta - t) + sum;
-    }
-    sum = t;
-  }
-  return sum + c;
-}
-
 template <typename i_t, typename f_t>
 DI thrust::pair<f_t, f_t> move_objective_score(
   const typename fj_t<i_t, f_t>::climber_data_t::view_t& fj, i_t var_idx, f_t delta)
@@ -176,7 +158,7 @@ __global__ void init_lhs_and_violation(typename fj_t<i_t, f_t>::climber_data_t::
         return fj.pb.coefficients[j] * fj.incumbent_assignment[fj.pb.variables[j]];
       });
     fj.incumbent_lhs[cstr_idx] =
-      KahanBabushkaNeumaierSum<i_t, f_t>(delta_it + offset_begin, delta_it + offset_end);
+      fj_kahan_babushka_neumaier_sum<i_t, f_t>(delta_it + offset_begin, delta_it + offset_end);
     fj.incumbent_lhs_sumcomp[cstr_idx] = 0;
 
     f_t th_violation       = fj.excess_score(cstr_idx, fj.incumbent_lhs[cstr_idx]);
@@ -457,7 +439,8 @@ DI bool check_feasibility(const typename fj_t<i_t, f_t>::climber_data_t::view_t&
         return fj.pb.coefficients[j] * fj.incumbent_assignment[fj.pb.variables[j]];
       });
 
-    f_t lhs = KahanBabushkaNeumaierSum<i_t, f_t>(delta_it + offset_begin, delta_it + offset_end);
+    f_t lhs =
+      fj_kahan_babushka_neumaier_sum<i_t, f_t>(delta_it + offset_begin, delta_it + offset_end);
     cuopt_assert(fj.cstr_satisfied(cIdx, lhs), "constraint violated");
   }
   cuopt_func_call((check_variable_feasibility<i_t, f_t>(fj, check_integer)));

cpp/src/mip/feasibility_jump/fj_cpu.cu

@@ -124,6 +124,22 @@ static inline bool tabu_check(fj_cpu_climber_t<i_t, f_t>& fj_cpu,
   }
 }
 
+template <typename i_t, typename f_t>
+static bool check_variable_feasibility(const typename fj_t<i_t, f_t>::climber_data_t::view_t& fj,
+                                       bool check_integer = true)
+{
+  for (i_t var_idx = 0; var_idx < fj.pb.n_variables; var_idx += 1) {
+    auto val      = fj.incumbent_assignment[var_idx];
+    bool feasible = fj.pb.check_variable_within_bounds(var_idx, val);
+
+    if (!feasible) return false;
+    if (check_integer && fj.pb.is_integer_var(var_idx) &&
+        !fj.pb.is_integer(fj.incumbent_assignment[var_idx]))
+      return false;
+  }
+  return true;
+}
+
 template <typename i_t, typename f_t>
 static inline fj_staged_score_t compute_score(fj_cpu_climber_t<i_t, f_t>& fj_cpu,
                                               i_t var_idx,
@@ -321,17 +337,22 @@ static void apply_move(fj_cpu_climber_t<i_t, f_t>& fj_cpu,
   fj_cpu.h_incumbent_objective += fj_cpu.h_obj_coeffs[var_idx] * delta;
   if (fj_cpu.h_incumbent_objective < fj_cpu.h_best_objective &&
       fj_cpu.violated_constraints.empty()) {
-    cuopt_assert(fj_cpu.satisfied_constraints.size() == fj_cpu.view.pb.n_constraints, "");
-    fj_cpu.h_best_objective =
-      fj_cpu.h_incumbent_objective - fj_cpu.settings.parameters.breakthrough_move_epsilon;
-    fj_cpu.h_best_assignment = fj_cpu.h_assignment;
-    CUOPT_LOG_TRACE("%sCPUFJ: new best objective: %g\n",
-                    fj_cpu.log_prefix.c_str(),
-                    fj_cpu.pb_ptr->get_user_obj_from_solver_obj(fj_cpu.h_best_objective));
-    if (fj_cpu.improvement_callback) {
-      fj_cpu.improvement_callback(fj_cpu.h_best_objective, fj_cpu.h_assignment);
+    // recompute the LHS values to cancel out accumulation errors, then check if feasibility remains
+    recompute_lhs(fj_cpu);
+
+    if (fj_cpu.violated_constraints.empty() && check_variable_feasibility<i_t, f_t>(fj_cpu.view)) {
+      cuopt_assert(fj_cpu.satisfied_constraints.size() == fj_cpu.view.pb.n_constraints, "");
+      fj_cpu.h_best_objective =
+        fj_cpu.h_incumbent_objective - fj_cpu.settings.parameters.breakthrough_move_epsilon;
+      fj_cpu.h_best_assignment = fj_cpu.h_assignment;
+      CUOPT_LOG_TRACE("%sCPUFJ: new best objective: %g\n",
+                      fj_cpu.log_prefix.c_str(),
+                      fj_cpu.pb_ptr->get_user_obj_from_solver_obj(fj_cpu.h_best_objective));
+      if (fj_cpu.improvement_callback) {
+        fj_cpu.improvement_callback(fj_cpu.h_best_objective, fj_cpu.h_assignment);
+      }
+      fj_cpu.feasible_found = true;
     }
-    fj_cpu.feasible_found = true;
   }
 
   i_t tabu_tenure = fj_cpu.settings.parameters.tabu_tenure_min +
@@ -537,23 +558,24 @@ static thrust::tuple<fj_move_t, fj_staged_score_t> find_mtm_move_sat(
 }
 
 template <typename i_t, typename f_t>
-static void init_lhs(fj_cpu_climber_t<i_t, f_t>& fj_cpu)
+static void recompute_lhs(fj_cpu_climber_t<i_t, f_t>& fj_cpu)
 {
   cuopt_assert(fj_cpu.h_lhs.size() == fj_cpu.view.pb.n_constraints, "h_lhs size mismatch");
 
   fj_cpu.violated_constraints.clear();
   fj_cpu.satisfied_constraints.clear();
   for (i_t cstr_idx = 0; cstr_idx < fj_cpu.view.pb.n_constraints; ++cstr_idx) {
     auto [offset_begin, offset_end] = fj_cpu.view.pb.range_for_constraint(cstr_idx);
-    f_t lhs                         = 0;
-    for (i_t i = offset_begin; i < offset_end; ++i) {
-      lhs += fj_cpu.h_coefficients[i] * fj_cpu.h_assignment[fj_cpu.h_variables[i]];
-    }
-
-    fj_cpu.h_lhs[cstr_idx] = lhs;
+    auto delta_it =
+      thrust::make_transform_iterator(thrust::make_counting_iterator(0), [fj = fj_cpu.view](i_t j) {
+        return fj.pb.coefficients[j] * fj.incumbent_assignment[fj.pb.variables[j]];
+      });
+    fj_cpu.h_lhs[cstr_idx] =
+      fj_kahan_babushka_neumaier_sum<i_t, f_t>(delta_it + offset_begin, delta_it + offset_end);
+    fj_cpu.h_lhs_sumcomp[cstr_idx] = 0;
 
     f_t cstr_tolerance = fj_cpu.view.get_corrected_tolerance(cstr_idx);
-    f_t new_cost       = fj_cpu.view.excess_score(cstr_idx, lhs);
+    f_t new_cost       = fj_cpu.view.excess_score(cstr_idx, fj_cpu.h_lhs[cstr_idx]);
     if (new_cost < -cstr_tolerance) {
       fj_cpu.violated_constraints.insert(cstr_idx);
     } else {
@@ -695,7 +717,7 @@ static void perturb(fj_cpu_climber_t<i_t, f_t>& fj_cpu)
     fj_cpu.h_assignment[var_idx] = val;
   }
 
-  init_lhs(fj_cpu);
+  recompute_lhs(fj_cpu);
 }
 
 template <typename i_t, typename f_t>
@@ -759,6 +781,8 @@ static void init_fj_cpu(fj_cpu_climber_t<i_t, f_t>& fj_cpu,
   fj_cpu.view.incumbent_assignment =
     raft::device_span<f_t>(fj_cpu.h_assignment.data(), fj_cpu.h_assignment.size());
   fj_cpu.view.incumbent_lhs = raft::device_span<f_t>(fj_cpu.h_lhs.data(), fj_cpu.h_lhs.size());
+  fj_cpu.view.incumbent_lhs_sumcomp =
+    raft::device_span<f_t>(fj_cpu.h_lhs_sumcomp.data(), fj_cpu.h_lhs_sumcomp.size());
   fj_cpu.view.tabu_nodec_until =
     raft::device_span<i_t>(fj_cpu.h_tabu_nodec_until.data(), fj_cpu.h_tabu_nodec_until.size());
   fj_cpu.view.tabu_noinc_until =
@@ -824,7 +848,7 @@ static void init_fj_cpu(fj_cpu_climber_t<i_t, f_t>& fj_cpu,
   fj_cpu.var_bitmap.resize(fj_cpu.view.pb.n_variables, false);
   fj_cpu.iter_mtm_vars.reserve(fj_cpu.view.pb.n_variables);
 
-  init_lhs(fj_cpu);
+  recompute_lhs(fj_cpu);
 }
 
 template <typename i_t, typename f_t>
@@ -914,6 +938,14 @@ bool fj_t<i_t, f_t>::cpu_solve(fj_cpu_climber_t<i_t, f_t>& fj_cpu, f_t in_time_l
       break;
     }
 
+    // periodically recompute the LHS and violation scores
+    // to correct any accumulated numerical errors
+    cuopt_assert(fj_cpu.settings.parameters.lhs_refresh_period > 0,
+                 "lhs_refresh_period should be positive");
+    if (fj_cpu.iterations % fj_cpu.settings.parameters.lhs_refresh_period == 0) {
+      recompute_lhs(fj_cpu);
+    }
+
     fj_move_t move          = fj_move_t{-1, 0};
     fj_staged_score_t score = fj_staged_score_t::invalid();
     // Perform lift moves

cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp

@@ -60,7 +60,7 @@ TEST_P(TimeLimitTestFixture, time_limit)
 
   // Dual simplex is spending some time for factorizing the basis, and this computation does not
   // check for time limit
-  double excess_allowed_time = 2.0;
+  double excess_allowed_time = 3.0;
   EXPECT_NEAR(solve_time, target_solve_time, excess_allowed_time);
 }
 INSTANTIATE_TEST_SUITE_P(
@@ -71,7 +71,7 @@ INSTANTIATE_TEST_SUITE_P(
                     5,
                     CUOPT_METHOD_DUAL_SIMPLEX),  // LP, Dual Simplex
     std::make_tuple("/linear_programming/square41/square41.mps", 5, CUOPT_METHOD_PDLP),  // LP, PDLP
-    std::make_tuple("/mip/supportcase22.mps", 5, CUOPT_METHOD_DUAL_SIMPLEX)              // MIP
+    std::make_tuple("/mip/supportcase22.mps", 15, CUOPT_METHOD_DUAL_SIMPLEX)             // MIP
     ));
 
 TEST(c_api, iteration_limit)