Fix issue with infinite lower bounds and try to bound free variables in barrier

cpp/cuopt_cli.cpp

@@ -13,6 +13,7 @@
 #include <cuopt/linear_programming/solve.hpp>
 #include <mps_parser/parser.hpp>
 #include <utilities/logger.hpp>
+#include <utilities/timer.hpp>
 
 #include <raft/core/device_setter.hpp>
 #include <raft/core/handle.hpp>
@@ -92,11 +93,15 @@ int run_single_file(const std::string& file_path,
                     bool solve_relaxation,
                     cuopt::linear_programming::solver_settings_t<int, double>& settings)
 {
+  cuopt::init_logger_t log(settings.get_parameter<std::string>(CUOPT_LOG_FILE),
+                           settings.get_parameter<bool>(CUOPT_LOG_TO_CONSOLE));
+
   std::string base_filename = file_path.substr(file_path.find_last_of("/\\") + 1);
 
   constexpr bool input_mps_strict = false;
   cuopt::mps_parser::mps_data_model_t<int, double> mps_data_model;
   bool parsing_failed = false;
+  auto timer          = cuopt::timer_t(settings.get_parameter<double>(CUOPT_TIME_LIMIT));
   {
     CUOPT_LOG_INFO("Reading file %s", base_filename.c_str());
     try {
@@ -111,6 +116,7 @@ int run_single_file(const std::string& file_path,
     CUOPT_LOG_ERROR("Parsing MPS failed. Exiting!");
     return -1;
   }
+  CUOPT_LOG_INFO("Read file %s in %.2f seconds", base_filename.c_str(), timer.elapsed_time());
 
   // Determine memory backend and create problem using interface
   // Create handle only for GPU memory backend (avoid CUDA init on CPU-only hosts)

cpp/src/barrier/barrier.cu

@@ -2151,10 +2151,9 @@ void barrier_solver_t<i_t, f_t>::gpu_compute_residual_norms(const rmm::device_uv
     std::max(device_vector_norm_inf<i_t, f_t>(data.d_primal_residual_, stream_view_),
              device_vector_norm_inf<i_t, f_t>(data.d_bound_residual_, stream_view_));
   dual_residual_norm = device_vector_norm_inf<i_t, f_t>(data.d_dual_residual_, stream_view_);
-  // TODO: CMM understand why rhs and not residual
   complementarity_residual_norm =
-    std::max(device_vector_norm_inf<i_t, f_t>(data.d_complementarity_xz_rhs_, stream_view_),
-             device_vector_norm_inf<i_t, f_t>(data.d_complementarity_wv_rhs_, stream_view_));
+    std::max(device_vector_norm_inf<i_t, f_t>(data.d_complementarity_xz_residual_, stream_view_),
+             device_vector_norm_inf<i_t, f_t>(data.d_complementarity_wv_residual_, stream_view_));
 }
 
 template <typename i_t, typename f_t>
@@ -3492,7 +3491,9 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time,
     f_t relative_primal_residual = primal_residual_norm / (1.0 + norm_b);
     f_t relative_dual_residual   = dual_residual_norm / (1.0 + norm_c);
     f_t relative_complementarity_residual =
-      complementarity_residual_norm / (1.0 + std::abs(primal_objective));
+      complementarity_residual_norm /
+      (1.0 + std::min(std::abs(compute_user_objective(lp, primal_objective)),
+                      std::abs(primal_objective)));
 
     dense_vector_t<i_t, f_t> upper(lp.upper);
     data.gather_upper_bounds(upper, data.restrict_u_);
@@ -3508,11 +3509,11 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time,
     float64_t elapsed_time = toc(start_time);
     settings.log.printf("%3d   %+.12e %+.12e %.2e %.2e %.2e %.1f\n",
                         iter,
-                        primal_objective,
-                        dual_objective,
-                        primal_residual_norm,
-                        dual_residual_norm,
-                        complementarity_residual_norm,
+                        compute_user_objective(lp, primal_objective),
+                        compute_user_objective(lp, dual_objective),
+                        relative_primal_residual,
+                        relative_dual_residual,
+                        relative_complementarity_residual,
                         elapsed_time);
 
     bool converged = primal_residual_norm < settings.barrier_relative_feasibility_tol &&
@@ -3654,7 +3655,9 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time,
       relative_primal_residual = primal_residual_norm / (1.0 + norm_b);
       relative_dual_residual   = dual_residual_norm / (1.0 + norm_c);
       relative_complementarity_residual =
-        complementarity_residual_norm / (1.0 + std::abs(primal_objective));
+        complementarity_residual_norm /
+        (1.0 + std::min(std::abs(compute_user_objective(lp, primal_objective)),
+                        std::abs(primal_objective)));
 
       if (relative_primal_residual < settings.barrier_relaxed_feasibility_tol &&
           relative_dual_residual < settings.barrier_relaxed_optimality_tol &&

cpp/src/dual_simplex/presolve.cpp

@@ -821,6 +821,168 @@ i_t presolve(const lp_problem_t<i_t, f_t>& original,
 {
   problem = original;
   std::vector<char> row_sense(problem.num_rows, '=');
+  // Check for free variables
+  i_t free_variables = 0;
+  for (i_t j = 0; j < problem.num_cols; j++) {
+    if (problem.lower[j] == -inf && problem.upper[j] == inf) { free_variables++; }
+  }
+
+  if (settings.barrier_presolve && free_variables > 0) {
+    // Try to remove free variables
+    std::vector<i_t> constraints_to_check;
+    std::vector<i_t> current_free_variables;
+    std::vector<i_t> row_marked(problem.num_rows, 0);
+    current_free_variables.reserve(problem.num_cols);
+    constraints_to_check.reserve(problem.num_rows);
+    for (i_t j = 0; j < problem.num_cols; j++) {
+      if (problem.lower[j] == -inf && problem.upper[j] == inf) {
+        current_free_variables.push_back(j);
+        const i_t col_start = problem.A.col_start[j];
+        const i_t col_end   = problem.A.col_start[j + 1];
+        for (i_t p = col_start; p < col_end; p++) {
+          const i_t i = problem.A.i[p];
+          if (row_marked[i] == 0) {
+            row_marked[i] = 1;
+            constraints_to_check.push_back(i);
+          }
+        }
+      }
+    }
+
+    i_t removed_free_variables = 0;
+
+    if (constraints_to_check.size() > 0) {
+      // Check if the constraints are feasible
+      csr_matrix_t<i_t, f_t> Arow(0, 0, 0);
+      problem.A.to_compressed_row(Arow);
+
+      // The constraints are in the form:
+      // sum_j a_j x_j = beta
+      for (i_t i : constraints_to_check) {
+        const i_t row_start   = Arow.row_start[i];
+        const i_t row_end     = Arow.row_start[i + 1];
+        f_t lower_activity_i  = 0.0;
+        f_t upper_activity_i  = 0.0;
+        i_t lower_inf_i       = 0;
+        i_t upper_inf_i       = 0;
+        i_t last_free_i       = -1;
+        f_t last_free_coeff_i = 0.0;
+        for (i_t p = row_start; p < row_end; p++) {
+          const i_t j       = Arow.j[p];
+          const f_t aij     = Arow.x[p];
+          const f_t lower_j = problem.lower[j];
+          const f_t upper_j = problem.upper[j];
+          if (lower_j == -inf && upper_j == inf) {
+            last_free_i       = j;
+            last_free_coeff_i = aij;
+          }
+          if (aij > 0) {
+            if (lower_j > -inf) {
+              lower_activity_i += aij * lower_j;
+            } else {
+              lower_inf_i++;
+            }
+            if (upper_j < inf) {
+              upper_activity_i += aij * upper_j;
+            } else {
+              upper_inf_i++;
+            }
+          } else {
+            if (upper_j < inf) {
+              lower_activity_i += aij * upper_j;
+            } else {
+              lower_inf_i++;
+            }
+            if (lower_j > -inf) {
+              upper_activity_i += aij * lower_j;
+            } else {
+              upper_inf_i++;
+            }
+          }
+        }
+
+        if (last_free_i == -1) { continue; }
+
+        // sum_j a_ij x_j == beta
+
+        const f_t rhs = problem.rhs[i];
+        // sum_{k != j} a_ik x_k + a_ij x_j == rhs
+        // Suppose that -inf < x_j < inf  and all other variables x_k with k != j are bounded
+        // a_ij x_j == rhs - sum_{k != j} a_ik x_k
+        // So if a_ij > 0, we have
+        //  x_j == 1/a_ij * (rhs - sum_{k != j} a_ik x_k)
+        // We can derive two bounds from  this:
+        // x_j <= 1/a_ij * (rhs - lower_activity_i) and
+        // x_j >= 1/a_ij * (rhs - upper_activity_i)
+
+        // If a_ij < 0, we have
+        // x_j == 1/a_ij * (rhs - sum_{k != j} a_ik x_k
+        // And we can derive two bounds from this:
+        // x_j >= 1/a_ij * (rhs - lower_activity_i)
+        // x_j <= 1/a_ij * (rhs - upper_activity_i)
+        const i_t j         = last_free_i;
+        const f_t a_ij      = last_free_coeff_i;
+        const f_t max_bound = 1e10;
+        bool bounded        = false;
+        if (a_ij > 0) {
+          if (lower_inf_i == 1) {
+            const f_t new_upper = 1.0 / a_ij * (rhs - lower_activity_i);
+            if (new_upper < max_bound) {
+              problem.upper[j] = new_upper;
+              bounded          = true;
+            }
+          }
+          if (upper_inf_i == 1) {
+            const f_t new_lower = 1.0 / a_ij * (rhs - upper_activity_i);
+            if (new_lower > -max_bound) {
+              problem.lower[j] = new_lower;
+              bounded          = true;
+            }
+          }
+        } else if (a_ij < 0) {
+          if (lower_inf_i == 1) {
+            const f_t new_lower = 1.0 / a_ij * (rhs - lower_activity_i);
+            if (new_lower > -max_bound) {
+              problem.lower[j] = new_lower;
+              bounded          = true;
+            }
+          }
+          if (upper_inf_i == 1) {
+            const f_t new_upper = 1.0 / a_ij * (rhs - upper_activity_i);
+            if (new_upper < max_bound) {
+              problem.upper[j] = new_upper;
+              bounded          = true;
+            }
+          }
+        }
+
+        if (bounded) { removed_free_variables++; }
+      }
+    }
+
+    for (i_t j : current_free_variables) {
+      if (problem.lower[j] > -inf && problem.upper[j] < inf) {
+        // We don't need two bounds. Pick the smallest one.
+        if (std::abs(problem.lower[j]) < std::abs(problem.upper[j])) {
+          // Restore the inf in the upper bound. Barrier will not require an additional w variable
+          problem.upper[j] = inf;
+        } else {
+          // Restores the -inf in the lower bound. Barrier will require an additional w variable
+          problem.lower[j] = -inf;
+        }
+      }
+    }
+
+    i_t new_free_variables = 0;
+    for (i_t j = 0; j < problem.num_cols; j++) {
+      if (problem.lower[j] == -inf && problem.upper[j] == inf) { new_free_variables++; }
+    }
+    if (removed_free_variables != 0) {
+      settings.log.printf("Bounded %d free variables\n", removed_free_variables);
+    }
+    assert(new_free_variables == free_variables - removed_free_variables);
+    free_variables = new_free_variables;
+  }
 
   // The original problem may have a variable without a lower bound
   // but a finite upper bound
@@ -834,7 +996,43 @@ i_t presolve(const lp_problem_t<i_t, f_t>& original,
     settings.log.printf("%d variables with no lower bound\n", no_lower_bound);
   }
 
-  // FIXME:: handle no lower bound case for barrier presolve
+  // Handle -inf < x_j <= u_j by substituting x'_j = -x_j, giving -u_j <= x'_j < inf
+  if (settings.barrier_presolve && no_lower_bound > 0) {
+    presolve_info.negated_variables.reserve(no_lower_bound);
+    for (i_t j = 0; j < problem.num_cols; j++) {
+      if (problem.lower[j] == -inf && problem.upper[j] < inf) {
+        presolve_info.negated_variables.push_back(j);
+
+        problem.lower[j] = -problem.upper[j];
+        problem.upper[j] = inf;
+        problem.objective[j] *= -1;
+
+        const i_t col_start = problem.A.col_start[j];
+        const i_t col_end   = problem.A.col_start[j + 1];
+        for (i_t p = col_start; p < col_end; p++) {
+          problem.A.x[p] *= -1.0;
+        }
+      }
+    }
+
+    // (1/2) x^T Q x with x = D x' (D_ii = -1 for negated columns) is (1/2) x'^T D Q D x'.
+    // One pass: Q'_{ik} = D_{ii} D_{kk} Q_{ik} — flip iff exactly one of {i,k} is negated.
+    if (problem.Q.n > 0 && !presolve_info.negated_variables.empty()) {
+      std::vector<bool> is_negated(static_cast<size_t>(problem.num_cols), false);
+      for (i_t const j : presolve_info.negated_variables) {
+        is_negated[static_cast<size_t>(j)] = true;
+      }
+      for (i_t row = 0; row < problem.Q.m; ++row) {
+        const i_t q_start         = problem.Q.row_start[row];
+        const i_t q_end           = problem.Q.row_start[row + 1];
+        const bool is_negated_row = is_negated[static_cast<size_t>(row)];
+        for (i_t p = q_start; p < q_end; ++p) {
+          const i_t col = problem.Q.j[p];
+          if (is_negated_row != is_negated[static_cast<size_t>(col)]) { problem.Q.x[p] *= -1.0; }
+        }
+      }
+    }
+  }
 
   // The original problem may have nonzero lower bounds
   // 0 != l_j <= x_j <= u_j
@@ -939,12 +1137,6 @@ i_t presolve(const lp_problem_t<i_t, f_t>& original,
     remove_empty_cols(problem, num_empty_cols, presolve_info);
   }
 
-  // Check for free variables
-  i_t free_variables = 0;
-  for (i_t j = 0; j < problem.num_cols; j++) {
-    if (problem.lower[j] == -inf && problem.upper[j] == inf) { free_variables++; }
-  }
-
   problem.Q.check_matrix("Before free variable expansion");
 
   if (settings.barrier_presolve && free_variables > 0) {
@@ -1511,6 +1703,14 @@ void uncrush_solution(const presolve_info_t<i_t, f_t>& presolve_info,
     }
     settings.log.printf("Post-solve: Handling removed lower bounds %d\n", num_lower_bounds);
   }
+
+  if (presolve_info.negated_variables.size() > 0) {
+    for (const i_t j : presolve_info.negated_variables) {
+      input_x[j] *= -1.0;
+      input_z[j] *= -1.0;
+    }
+  }
+
   assert(uncrushed_x.size() == input_x.size());
   assert(uncrushed_y.size() == input_y.size());
   assert(uncrushed_z.size() == input_z.size());

cpp/src/dual_simplex/presolve.hpp

@@ -202,6 +202,9 @@ struct presolve_info_t {
   std::vector<i_t> removed_constraints;
 
   folding_info_t<i_t, f_t> folding_info;
+
+  // Variables that were negated to handle -inf < x_j <= u_j
+  std::vector<i_t> negated_variables;
 };
 
 template <typename i_t, typename f_t>