diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index c5ef847106..c2ff65a59d 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -851,6 +851,12 @@ i_t presolve(const lp_problem_t& original, i_t removed_free_variables = 0; + // Track which constraint provided each implied bound for dual correction + std::vector lower_bound_constraint(problem.num_cols, -1); + std::vector lower_bound_coefficient(problem.num_cols, 0.0); + std::vector upper_bound_constraint(problem.num_cols, -1); + std::vector upper_bound_coefficient(problem.num_cols, 0.0); + if (constraints_to_check.size() > 0) { // Check if the constraints are feasible csr_matrix_t Arow(0, 0, 0); @@ -928,30 +934,38 @@ i_t presolve(const lp_problem_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 +987,24 @@ i_t presolve(const lp_problem_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}); + } + } + 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++; } @@ -1562,6 +1594,7 @@ void uncrush_dual_solution(const user_problem_t& user_problem, template void uncrush_solution(const presolve_info_t& presolve_info, const simplex_solver_settings_t& settings, + const lp_problem_t& original_problem, const std::vector& crushed_x, const std::vector& crushed_y, const std::vector& crushed_z, @@ -1711,6 +1744,32 @@ void uncrush_solution(const presolve_info_t& presolve_info, } } + // Dual correction for originally-free variables that received implied bounds. + // Barrier produced (u, w) with w_j != 0 satisfying A^T u + w = c + Qx. + // We need corrected (y, z) with z_j = 0: set du = w_j / a_{i*,j}, then + // y_{i*} += du and z_k -= a_{i*,k} * du for all k in constraint i*. + if (!presolve_info.bounded_free_variables.empty()) { + settings.log.printf("Post-solve: Correcting duals for %d bounded free variables\n", + static_cast(presolve_info.bounded_free_variables.size())); + const csc_matrix_t& A = original_problem.A; + for (const auto& bfv : presolve_info.bounded_free_variables) { + 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()); @@ -1769,6 +1828,7 @@ template void uncrush_dual_solution(const user_problem_t(const presolve_info_t& presolve_info, const simplex_solver_settings_t& settings, + const lp_problem_t& original_problem, const std::vector& crushed_x, const std::vector& crushed_y, const std::vector& crushed_z, diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index d0e2d52812..65ee09427d 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/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 +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 struct presolve_info_t { // indices of variables in the original problem that remain in the presolved problem @@ -153,6 +162,9 @@ struct presolve_info_t { // Variables that were negated to handle -inf < x_j <= u_j std::vector negated_variables; + + // Originally-free variables that received implied bounds, with the constraint used + std::vector> bounded_free_variables; }; template @@ -241,6 +253,7 @@ void uncrush_dual_solution(const user_problem_t& user_problem, template void uncrush_solution(const presolve_info_t& presolve_info, const simplex_solver_settings_t& settings, + const lp_problem_t& original_problem, const std::vector& crushed_x, const std::vector& crushed_y, const std::vector& crushed_z, diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 82d922eec3..5ae6fd7f9b 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -297,6 +297,7 @@ lp_status_t solve_linear_program_with_advanced_basis( unscale_solution(column_scales, solution.x, solution.z, unscaled_x, unscaled_z); uncrush_solution(presolve_info, settings, + original_lp, unscaled_x, solution.y, unscaled_z, @@ -439,6 +440,7 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us // Undo presolve uncrush_solution(presolve_info, barrier_settings, + original_lp, unscaled_x, barrier_solution.y, unscaled_z, diff --git a/cpp/tests/dual_simplex/unit_tests/solve.cpp b/cpp/tests/dual_simplex/unit_tests/solve.cpp index 7aed72fe0f..f37aa25bff 100644 --- a/cpp/tests/dual_simplex/unit_tests/solve.cpp +++ b/cpp/tests/dual_simplex/unit_tests/solve.cpp @@ -17,11 +17,13 @@ #include #include +#include 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 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 diff --git a/cpp/tests/dual_simplex/unit_tests/solve_barrier.cu b/cpp/tests/dual_simplex/unit_tests/solve_barrier.cu index abfe37c9fd..5b671f7d4d 100644 --- a/cpp/tests/dual_simplex/unit_tests/solve_barrier.cu +++ b/cpp/tests/dual_simplex/unit_tests/solve_barrier.cu @@ -8,9 +8,13 @@ #include #include +#include #include +#include +#include +#include #include #include #include @@ -20,6 +24,7 @@ #include #include +#include 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::mps_parser::parse_mps(path); + + auto settings = cuopt::linear_programming::pdlp_solver_settings_t{}; + + 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 diff --git a/datasets/quadratic_programming/min_x_squared.mps b/datasets/quadratic_programming/min_x_squared.mps new file mode 100644 index 0000000000..c37b7ea373 --- /dev/null +++ b/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