Fix a bug in dual variables and reduced costs when we add implied bounds on free variables

cpp/src/dual_simplex/presolve.cpp

@@ -851,6 +851,13 @@ i_t presolve(const lp_problem_t<i_t, f_t>& original,
     }
 
     i_t removed_free_variables = 0;
+
+    // Track which constraint provided each implied bound for dual correction
+    std::vector<i_t> lower_bound_constraint(problem.num_cols, -1);
+    std::vector<f_t> lower_bound_coefficient(problem.num_cols, 0.0);
+    std::vector<i_t> upper_bound_constraint(problem.num_cols, -1);
+    std::vector<f_t> upper_bound_coefficient(problem.num_cols, 0.0);
+
     if (!constraints_to_check.empty()) {
       // Check if the constraints are feasible
       csr_matrix_t<i_t, f_t> Arow(0, 0, 0);
@@ -928,30 +935,38 @@ i_t presolve(const lp_problem_t<i_t, f_t>& original,
           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;
+              problem.upper[j]           = new_upper;
+              upper_bound_constraint[j]  = i;
+              upper_bound_coefficient[j] = a_ij;
+              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;
+              problem.lower[j]           = new_lower;
+              lower_bound_constraint[j]  = i;
+              lower_bound_coefficient[j] = a_ij;
+              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;
+              problem.lower[j]           = new_lower;
+              lower_bound_constraint[j]  = i;
+              lower_bound_coefficient[j] = a_ij;
+              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;
+              problem.upper[j]           = new_upper;
+              upper_bound_constraint[j]  = i;
+              upper_bound_coefficient[j] = a_ij;
+              bounded                    = true;
             }
           }
         }
@@ -973,6 +988,24 @@ i_t presolve(const lp_problem_t<i_t, f_t>& original,
       }
     }
 
+    // Record bounded free variables for dual correction in uncrush.
+    // After the keep-one-bound logic, each bounded variable has exactly one finite bound.
+    for (i_t j : current_free_variables) {
+      i_t bounding_constraint  = -1;
+      f_t bounding_coefficient = 0.0;
+      if (problem.lower[j] > -inf && lower_bound_constraint[j] != -1) {
+        bounding_constraint  = lower_bound_constraint[j];
+        bounding_coefficient = lower_bound_coefficient[j];
+      } else if (problem.upper[j] < inf && upper_bound_constraint[j] != -1) {
+        bounding_constraint  = upper_bound_constraint[j];
+        bounding_coefficient = upper_bound_coefficient[j];
+      }
+      if (bounding_constraint != -1) {
+        presolve_info.bounded_free_variables.push_back(
+          {j, bounding_constraint, bounding_coefficient});
+      }
+    }
+
     free_variables = 0;
     for (i_t j = 0; j < problem.num_cols; j++) {
       if (problem.lower[j] == -inf && problem.upper[j] == inf) { free_variables++; }
@@ -1571,6 +1604,7 @@ void uncrush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
 template <typename i_t, typename f_t>
 void uncrush_solution(const presolve_info_t<i_t, f_t>& presolve_info,
                       const simplex_solver_settings_t<i_t, f_t>& settings,
+                      const lp_problem_t<i_t, f_t>& original_problem,
                       const std::vector<f_t>& crushed_x,
                       const std::vector<f_t>& crushed_y,
                       const std::vector<f_t>& crushed_z,
@@ -1720,6 +1754,46 @@ void uncrush_solution(const presolve_info_t<i_t, f_t>& presolve_info,
     }
   }
 
+  // Dual correction for originally free variables that received implied bounds.
+  // Barrier produced (y, z) with z_j != 0 satisfying A^T y + z = c + Qx.
+  // We need corrected (y_bar, z_bar) with z_bar_j = 0 for all j in F_b where
+  // F_b = { j | x_j was initially free and is now bounded }
+  //
+  // For a given j_f in F_b, let i* be the constraint that implied the bound on x_j_f.
+  // Compute the scalar delta_u = z_j_f / a_{i*,j_f}.
+  // Set y_bar = y + delta_u e_i* and
+  // Let R_i* = { j | a_{i*, j} != 0 }
+  // z_bar_j = z_j - delta_u a_{i*,j} for all j in R_i*.
+  // z_bar_j = z_j                    for all j not in R_i*.
+  //
+  // Then you can show that A^T y_bar + z_bar = c + Qx and
+  // z_bar_{j_f} = 0.
+  if (!presolve_info.bounded_free_variables.empty()) {
+    settings.log.printf("Post-solve: Correcting duals for %d bounded free variables\n",
+                        static_cast<i_t>(presolve_info.bounded_free_variables.size()));
+    const csc_matrix_t<i_t, f_t>& A = original_problem.A;
+    // Traverse in reverse order, to ensure that all z_j = 0 after the correction
+    for (auto it = presolve_info.bounded_free_variables.rbegin();
+         it != presolve_info.bounded_free_variables.rend();
+         ++it) {
+      const auto& bfv = *it;
+      const f_t w_j   = input_z[bfv.variable];
+      if (w_j == 0.0) { continue; }
+      const f_t du = w_j / bfv.coefficient;
+      input_y[bfv.constraint] += du;
+      for (i_t j = 0; j < A.n; j++) {
+        const i_t col_start = A.col_start[j];
+        const i_t col_end   = A.col_start[j + 1];
+        for (i_t p = col_start; p < col_end; p++) {
+          if (A.i[p] == bfv.constraint) {
+            input_z[j] -= A.x[p] * du;
+            break;
+          }
+        }
+      }
+    }
+  }
+
   assert(uncrushed_x.size() == input_x.size());
   assert(uncrushed_y.size() == input_y.size());
   assert(uncrushed_z.size() == input_z.size());
@@ -1778,6 +1852,7 @@ template void uncrush_dual_solution<int, double>(const user_problem_t<int, doubl
 
 template void uncrush_solution<int, double>(const presolve_info_t<int, double>& presolve_info,
                                             const simplex_solver_settings_t<int, double>& settings,
+                                            const lp_problem_t<int, double>& original_problem,
                                             const std::vector<double>& crushed_x,
                                             const std::vector<double>& crushed_y,
                                             const std::vector<double>& crushed_z,

cpp/src/dual_simplex/presolve.hpp

@@ -129,6 +129,15 @@ struct folding_info_t {
   bool is_folded;
 };
 
+// Free variable that received an implied bound during presolve.
+// Stores the bounding constraint and coefficient for dual correction in uncrush.
+template <typename i_t, typename f_t>
+struct bounded_free_var_t {
+  i_t variable;     // j: the originally-free variable
+  i_t constraint;   // i*: the constraint that implied the bound
+  f_t coefficient;  // a_{i*,j}: the coefficient of x_j in constraint i*
+};
+
 template <typename i_t, typename f_t>
 struct presolve_info_t {
   // indices of variables in the original problem that remain in the presolved problem
@@ -153,7 +162,8 @@ struct presolve_info_t {
 
   // Variables that were negated to handle -inf < x_j <= u_j
   std::vector<i_t> negated_variables;
-
+  // Originally-free variables that received implied bounds, with the constraint used
+  std::vector<bounded_free_var_t<i_t, f_t>> bounded_free_variables;
   // Free variable indices for QP augmented system (not split, handled natively)
   std::vector<i_t> free_variable_indices;
 };
@@ -244,6 +254,7 @@ void uncrush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
 template <typename i_t, typename f_t>
 void uncrush_solution(const presolve_info_t<i_t, f_t>& presolve_info,
                       const simplex_solver_settings_t<i_t, f_t>& settings,
+                      const lp_problem_t<i_t, f_t>& original_problem,
                       const std::vector<f_t>& crushed_x,
                       const std::vector<f_t>& crushed_y,
                       const std::vector<f_t>& crushed_z,

cpp/src/dual_simplex/solve.cpp

@@ -297,6 +297,7 @@ lp_status_t solve_linear_program_with_advanced_basis(
       unscale_solution<i_t, f_t>(column_scales, solution.x, solution.z, unscaled_x, unscaled_z);
       uncrush_solution(presolve_info,
                        settings,
+                       original_lp,
                        unscaled_x,
                        solution.y,
                        unscaled_z,
@@ -420,6 +421,7 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t<i_t, f_t>& us
     // Undo presolve
     uncrush_solution(presolve_info,
                      barrier_settings,
+                     original_lp,
                      unscaled_x,
                      barrier_solution.y,
                      unscaled_z,

cpp/tests/dual_simplex/unit_tests/solve.cpp

@@ -17,11 +17,13 @@
 #include <dual_simplex/user_problem.hpp>
 
 #include <cuopt/linear_programming/io/parser.hpp>
+#include <utilities/logger.hpp>
 
 namespace cuopt::linear_programming::dual_simplex::test {
 
 TEST(dual_simplex, chess_set)
 {
+  cuopt::init_logger_t log("", true);
   namespace dual_simplex = cuopt::linear_programming::dual_simplex;
   raft::handle_t handle{};
   dual_simplex::user_problem_t<int, double> user_problem(&handle);
@@ -95,6 +97,7 @@ TEST(dual_simplex, chess_set)
 
 TEST(dual_simplex, burglar)
 {
+  cuopt::init_logger_t log("", true);
   constexpr int num_items     = 8;
   constexpr double max_weight = 102;
 
@@ -169,6 +172,7 @@ TEST(dual_simplex, burglar)
 
 TEST(dual_simplex, empty_columns)
 {
+  cuopt::init_logger_t log("", true);
   // Same as burglar problem above but with an empty column inserted
   constexpr int num_items     = 9;
   constexpr double max_weight = 102;
@@ -257,6 +261,7 @@ TEST(dual_simplex, empty_columns)
 
 TEST(dual_simplex, dual_variable_greater_than)
 {
+  cuopt::init_logger_t log("", true);
   // minimize   3*x0 + 2 * x1
   // subject to  x0 + x1  >= 1
   //             x0 + 2x1 >= 3

cpp/tests/dual_simplex/unit_tests/solve_barrier.cu

@@ -8,9 +8,13 @@
 #include <cstdio>
 
 #include <utilities/common_utils.hpp>
+#include <utilities/copy_helpers.hpp>
 
 #include <gtest/gtest.h>
 
+#include <cuopt/linear_programming/constants.h>
+#include <cuopt/linear_programming/pdlp/solver_settings.hpp>
+#include <cuopt/linear_programming/solve.hpp>
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/solve.hpp>
 #include <dual_simplex/tic_toc.hpp>
@@ -20,6 +24,7 @@
 #include <raft/core/cusparse_macros.hpp>
 
 #include <cuopt/linear_programming/io/parser.hpp>
+#include <utilities/logger.hpp>
 
 namespace cuopt::linear_programming::dual_simplex::test {
 
@@ -35,6 +40,7 @@ static void init_handler(const raft::handle_t* handle_ptr)
 
 TEST(barrier, chess_set)
 {
+  cuopt::init_logger_t log("", true);
   namespace dual_simplex = cuopt::linear_programming::dual_simplex;
   raft::handle_t handle{};
   init_handler(&handle);
@@ -104,6 +110,7 @@ TEST(barrier, chess_set)
 
 TEST(barrier, dual_variable_greater_than)
 {
+  cuopt::init_logger_t log("", true);
   // minimize   3*x0 + 2 * x1
   // subject to  x0 + x1  >= 1
   //             x0 + 2x1 >= 3
@@ -174,4 +181,39 @@ TEST(barrier, dual_variable_greater_than)
   EXPECT_NEAR(solution.z[1], 0.0, 1e-5);
 }
 
+TEST(barrier, min_x_squared_free_variable_dual_correction)
+{
+  // minimize   x^2         (Q = [2.0], so 0.5 * x^T Q x = x^2)
+  // subject to x >= 1
+  // x is free
+  //
+  // Optimal: x = 1, obj = 1, y[0] = 2, z[0] = 0
+  // This tests the dual correction for originally-free variables that
+  // received implied bounds during presolve.
+
+  const raft::handle_t handle{};
+  init_handler(&handle);
+
+  auto path =
+    cuopt::test::get_rapids_dataset_root_dir() + "/quadratic_programming/min_x_squared.mps";
+  auto mps_data = cuopt::linear_programming::io::parse_mps<int, double>(path);
+
+  auto settings = cuopt::linear_programming::pdlp_solver_settings_t<int, double>{};
+
+  auto solution = cuopt::linear_programming::solve_lp(&handle, mps_data, settings);
+
+  EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERMINATION_STATUS_OPTIMAL);
+
+  auto h_x = cuopt::host_copy(solution.get_primal_solution(), handle.get_stream());
+  auto h_y = cuopt::host_copy(solution.get_dual_solution(), handle.get_stream());
+  auto h_z = cuopt::host_copy(solution.get_reduced_cost(), handle.get_stream());
+
+  printf("x %e y %e z %e\n", h_x[0], h_y[0], h_z[0]);
+
+  const double tol = 1e-5;
+  EXPECT_NEAR(h_x[0], 1.0, tol);
+  EXPECT_NEAR(h_y[0], 2.0, tol);
+  EXPECT_NEAR(h_z[0], 0.0, tol);
+}
+
 }  // namespace cuopt::linear_programming::dual_simplex::test

datasets/quadratic_programming/min_x_squared.mps

@@ -0,0 +1,13 @@
+NAME        min_x_squared
+ROWS
+ N  obj
+ G  c1
+COLUMNS
+    x         c1        1
+RHS
+    RHS_V     c1        1
+BOUNDS
+ FR BOUND     x
+QUADOBJ
+    x         x         2
+ENDATA