From f622513235888c31e14cd884c6ea7743f223a8b5 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 22 May 2025 10:32:41 -0700 Subject: [PATCH 1/3] Bug fix: solution incorrect when dual simplex presolve removes empty cols --- cpp/src/dual_simplex/presolve.cpp | 61 +++++++++++++++++++++++++++++-- cpp/src/dual_simplex/presolve.hpp | 22 ++++++++++- cpp/src/dual_simplex/scaling.cpp | 22 +++++++++++ cpp/src/dual_simplex/scaling.hpp | 7 ++++ cpp/src/dual_simplex/solve.cpp | 13 ++++--- 5 files changed, 114 insertions(+), 11 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 11b86348a..55e84dea3 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -24,7 +24,9 @@ namespace cuopt::linear_programming::dual_simplex { template -i_t remove_empty_cols(lp_problem_t& problem, i_t& num_empty_cols) +i_t remove_empty_cols(lp_problem_t& problem, + i_t& num_empty_cols, + presolve_info_t& presolve_info) { constexpr bool verbose = false; if (verbose) { printf("Removing %d empty columns\n", num_empty_cols); } @@ -35,15 +37,22 @@ i_t remove_empty_cols(lp_problem_t& problem, i_t& num_empty_cols) // sum_{k != j} c_k * x_k + c_j * l_j if c_j > 0 // or // sum_{k != j} c_k * x_k + c_j * u_j if c_j < 0 + presolve_info.removed_variables.reserve(num_empty_cols); + presolve_info.removed_values.reserve(num_empty_cols); + presolve_info.removed_reduced_costs.reserve(num_empty_cols); std::vector col_marker(problem.num_cols); i_t new_cols = 0; for (i_t j = 0; j < problem.num_cols; ++j) { if ((problem.A.col_start[j + 1] - problem.A.col_start[j]) == 0) { col_marker[j] = 1; + presolve_info.removed_variables.push_back(j); + presolve_info.removed_reduced_costs.push_back(problem.objective[j]); if (problem.objective[j] >= 0) { + presolve_info.removed_values.push_back(problem.lower[j]); problem.obj_constant += problem.objective[j] * problem.lower[j]; assert(problem.lower[j] > -inf); } else { + presolve_info.removed_values.push_back(problem.upper[j]); problem.obj_constant += problem.objective[j] * problem.upper[j]; assert(problem.upper[j] < inf); } @@ -52,6 +61,7 @@ i_t remove_empty_cols(lp_problem_t& problem, i_t& num_empty_cols) new_cols++; } } + presolve_info.remaining_variables.reserve(new_cols); problem.A.remove_columns(col_marker); // Clean up objective, lower, upper, and col_names @@ -66,6 +76,7 @@ i_t remove_empty_cols(lp_problem_t& problem, i_t& num_empty_cols) objective[new_j] = problem.objective[j]; lower[new_j] = problem.lower[j]; upper[new_j] = problem.upper[j]; + presolve_info.remaining_variables.push_back(j); new_j++; } else { num_empty_cols--; @@ -574,7 +585,8 @@ void convert_user_problem(const user_problem_t& user_problem, template i_t presolve(const lp_problem_t& original, const simplex_solver_settings_t& settings, - lp_problem_t& problem) + lp_problem_t& problem, + presolve_info_t& presolve_info) { problem = original; std::vector row_sense(problem.num_rows, '='); @@ -595,6 +607,7 @@ i_t presolve(const lp_problem_t& original, } } if (num_empty_rows > 0) { + settings.log.printf("Presolve removing %d empty rows\n", num_empty_rows); i_t i = remove_empty_rows(problem, row_sense, num_empty_rows); if (i != 0) { return -1; } } @@ -606,7 +619,10 @@ i_t presolve(const lp_problem_t& original, if ((problem.A.col_start[j + 1] - problem.A.col_start[j]) == 0) { num_empty_cols++; } } } - if (num_empty_cols > 0) { remove_empty_cols(problem, num_empty_cols); } + if (num_empty_cols > 0) { + settings.log.printf("Presolve removing %d empty cols\n", num_empty_cols); + remove_empty_cols(problem, num_empty_cols, presolve_info); + } // Check for dependent rows constexpr bool check_dependent_rows = false; @@ -826,6 +842,36 @@ void uncrush_primal_solution(const user_problem_t& user_problem, std::copy(solution.begin(), solution.begin() + user_problem.num_cols, user_solution.data()); } +template +void uncrush_solution(const presolve_info_t& presolve_info, + const std::vector& crushed_x, + const std::vector& crushed_z, + std::vector& uncrushed_x, + std::vector& uncrushed_z) +{ + if (presolve_info.removed_variables.size() == 0) { + uncrushed_x = crushed_x; + uncrushed_z = crushed_z; + return; + } + + const i_t n = presolve_info.removed_variables.size() + presolve_info.remaining_variables.size(); + uncrushed_x.resize(n); + uncrushed_z.resize(n); + + i_t k = 0; + for (const i_t j : presolve_info.remaining_variables) { + uncrushed_x[j] = crushed_x[k]; + uncrushed_z[j] = crushed_z[k]; + k++; + } + + for (const i_t j : presolve_info.removed_variables) { + uncrushed_x[j] = presolve_info.removed_values[j]; + uncrushed_z[j] = presolve_info.removed_reduced_costs[j]; + } +} + #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE template void convert_user_problem(const user_problem_t& user_problem, @@ -841,7 +887,9 @@ template void convert_user_lp_with_guess( template int presolve(const lp_problem_t& original, const simplex_solver_settings_t& settings, - lp_problem_t& presolved); + lp_problem_t& presolved, + presolve_info_t& presolve_info); + template void crush_primal_solution(const user_problem_t& user_problem, const lp_problem_t& problem, const std::vector& user_solution, @@ -853,6 +901,11 @@ template void uncrush_primal_solution(const user_problem_t& solution, std::vector& user_solution); +template void uncrush_solution(const presolve_info_t& presolve_info, + const std::vector& crushed_x, + const std::vector& crushed_z, + std::vector& uncrushed_x, + std::vector& uncrushed_z); #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 4b877caef..e6eb542ac 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -49,6 +49,18 @@ struct lp_problem_t { f_t obj_scale; // 1.0 for min, -1.0 for max }; +template +struct presolve_info_t { + // indices of variables in the original problem that remain in the presolved problem + std::vector remaining_variables; + // indicies of variables in the original problem that have been removed in the presolved problem + std::vector removed_variables; + // values of the removed variables + std::vector removed_values; + // values of the removed reduced costs + std::vector removed_reduced_costs; +}; + template void convert_user_problem(const user_problem_t& user_problem, lp_problem_t& problem, @@ -70,7 +82,8 @@ void convert_user_lp_with_guess(const user_problem_t& user_problem, template i_t presolve(const lp_problem_t& original, const simplex_solver_settings_t& settings, - lp_problem_t& presolved); + lp_problem_t& presolved, + presolve_info_t& presolve_info); template void crush_primal_solution(const user_problem_t& user_problem, @@ -102,4 +115,11 @@ void uncrush_primal_solution(const user_problem_t& user_problem, const std::vector& solution, std::vector& user_solution); +template +void uncrush_solution(const presolve_info_t& presolve_info, + const std::vector& crushed_x, + const std::vector& crushed_z, + std::vector& uncrushed_x, + std::vector& uncrushed_z); + } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/scaling.cpp b/cpp/src/dual_simplex/scaling.cpp index bc33a58e1..b158dadef 100644 --- a/cpp/src/dual_simplex/scaling.cpp +++ b/cpp/src/dual_simplex/scaling.cpp @@ -76,6 +76,22 @@ i_t column_scaling(const lp_problem_t& unscaled, return 0; } +template +void unscale_solution(const std::vector& column_scaling, + const std::vector& scaled_x, + const std::vector& scaled_z, + std::vector& unscaled_x, + std::vector& unscaled_z) +{ + const i_t n = scaled_x.size(); + unscaled_x.resize(n); + unscaled_z.resize(n); + for (i_t j = 0; j < n; ++j) { + unscaled_x[j] = scaled_x[j] / column_scaling[j]; + unscaled_z[j] = scaled_z[j] / column_scaling[j]; + } +} + #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE template int column_scaling(const lp_problem_t& unscaled, @@ -83,6 +99,12 @@ template int column_scaling(const lp_problem_t& unscal lp_problem_t& scaled, std::vector& column_scaling); +template void unscale_solution(const std::vector& column_scaling, + const std::vector& scaled_x, + const std::vector& scaled_z, + std::vector& unscaled_x, + std::vector& unscaled_z); + #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/scaling.hpp b/cpp/src/dual_simplex/scaling.hpp index 09f463080..898a001c6 100644 --- a/cpp/src/dual_simplex/scaling.hpp +++ b/cpp/src/dual_simplex/scaling.hpp @@ -32,4 +32,11 @@ i_t column_scaling(const lp_problem_t& unscaled, lp_problem_t& scaled, std::vector& column_scaling); +template +void unscale_solution(const std::vector& column_scaling, + const std::vector& scaled_x, + const std::vector& scaled_z, + std::vector& unscaled_x, + std::vector& unscaled_z); + } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 05c1512fa..11a0f6d42 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -117,7 +117,8 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t& original { lp_status_t lp_status = lp_status_t::UNSET; lp_problem_t presolved_lp(1, 1, 1); - const i_t ok = presolve(original_lp, settings, presolved_lp); + presolve_info_t presolve_info; + const i_t ok = presolve(original_lp, settings, presolved_lp, presolve_info); if (ok == -1) { return lp_status_t::INFEASIBLE; } constexpr bool write_out_matlab = false; @@ -199,11 +200,11 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t& original primal_phase2(2, start_time, lp, settings, vstatus, solution, iter); } if (status == dual::status_t::OPTIMAL) { - // Unscale solution - for (i_t j = 0; j < original_lp.num_cols; j++) { - original_solution.x[j] = solution.x[j] / column_scales[j]; - original_solution.z[j] = solution.z[j] / column_scales[j]; - } + std::vector unscaled_x(lp.num_cols); + std::vector unscaled_z(lp.num_cols); + unscale_solution(column_scales, solution.x, solution.z, unscaled_x, unscaled_z); + uncrush_solution( + presolve_info, unscaled_x, unscaled_z, original_solution.x, original_solution.z); original_solution.y = solution.y; original_solution.objective = solution.objective; original_solution.user_objective = solution.user_objective; From 769771f04ed0ad7843de77bf7f2265bbae9d7eac Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 22 May 2025 14:57:59 -0700 Subject: [PATCH 2/3] Fix loop over removed variables. Add unit test. --- cpp/src/dual_simplex/presolve.cpp | 6 +- cpp/tests/dual_simplex/unit_tests/solve.cpp | 87 +++++++++++++++++++++ 2 files changed, 91 insertions(+), 2 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 55e84dea3..e87c88b40 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -866,9 +866,11 @@ void uncrush_solution(const presolve_info_t& presolve_info, k++; } + k = 0; for (const i_t j : presolve_info.removed_variables) { - uncrushed_x[j] = presolve_info.removed_values[j]; - uncrushed_z[j] = presolve_info.removed_reduced_costs[j]; + uncrushed_x[j] = presolve_info.removed_values[k]; + uncrushed_z[j] = presolve_info.removed_reduced_costs[k]; + k++; } } diff --git a/cpp/tests/dual_simplex/unit_tests/solve.cpp b/cpp/tests/dual_simplex/unit_tests/solve.cpp index 9664e0ac7..3fbfc3766 100644 --- a/cpp/tests/dual_simplex/unit_tests/solve.cpp +++ b/cpp/tests/dual_simplex/unit_tests/solve.cpp @@ -175,4 +175,91 @@ TEST(dual_simplex, burglar) EXPECT_NEAR(solution[5], 1, 1e-6); } +TEST(dual_simplex, empty_columns) +{ + // Same as burglar problem above but with an empty column inserted + constexpr int num_items = 9; + constexpr double max_weight = 102; + + std::vector value({15, 100, 90, 0, 60, 40, 15, 10, 1}); + std::vector weight({2, 20, 20, 0, 30, 40, 30, 60, 10}); + + // maximize sum_i value[i] * take[i] + // sum_i weight[i] * take[i] <= max_weight + // take[i] binary for all i + + cuopt::linear_programming::dual_simplex::user_problem_t user_problem; + constexpr int m = 1; + constexpr int n = num_items; + constexpr int nz = num_items - 1; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective.resize(n); + for (int j = 0; j < num_items; ++j) { + user_problem.objective[j] = -value[j]; + } + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start.resize(n + 1); + int nnz = 0; + for (int j = 0; j < num_items; ++j) { + user_problem.A.col_start[j] = nnz; + if (weight[j] > 0) { + user_problem.A.i[nnz] = 0; + user_problem.A.x[nnz] = weight[j]; + nnz++; + } + } + user_problem.A.col_start[n] = nnz; + user_problem.rhs.resize(m); + user_problem.rhs[0] = max_weight; + user_problem.row_sense.resize(m); + user_problem.row_sense[0] = 'L'; + user_problem.lower.resize(n); + user_problem.upper.resize(n); + for (int j = 0; j < num_items; ++j) { + user_problem.lower[j] = 0.0; + user_problem.upper[j] = 1.0; + } + user_problem.num_range_rows = 0; + user_problem.problem_name = "burglar"; + user_problem.row_names.resize(m); + user_problem.row_names[0] = "weight restriction"; + user_problem.col_names.resize(n); + for (int j = 0; j < num_items; ++j) { + user_problem.col_names[j] = "x"; + } + user_problem.obj_constant = 0.0; + user_problem.var_types.resize(n); + for (int j = 0; j < num_items; ++j) { + user_problem.var_types[j] = + cuopt::linear_programming::dual_simplex::variable_type_t::CONTINUOUS; + } + + cuopt::linear_programming::dual_simplex::simplex_solver_settings_t settings; + + cuopt::linear_programming::dual_simplex::lp_solution_t solution( + user_problem.num_rows, user_problem.num_cols); + EXPECT_EQ((cuopt::linear_programming::dual_simplex::solve_linear_program( + user_problem, settings, solution)), + cuopt::linear_programming::dual_simplex::lp_status_t::OPTIMAL); + double objective = 0.0; + for (int j = 0; j < num_items; ++j) { + objective += value[j] * solution.x[j]; + } + EXPECT_NEAR(objective, 295, 1e-6); + EXPECT_NEAR(solution.x[0], 1, 1e-6); + EXPECT_NEAR(solution.x[1], 1, 1e-6); + EXPECT_NEAR(solution.x[2], 1, 1e-6); + EXPECT_NEAR(solution.x[3], 0, 1e-6); + EXPECT_NEAR(solution.x[4], 1, 1e-6); + EXPECT_NEAR(solution.x[5], 0.75, 1e-6); + EXPECT_NEAR(solution.x[6], 0, 1e-6); + EXPECT_NEAR(solution.x[7], 0, 1e-6); + EXPECT_NEAR(solution.x[8], 0, 1e-6); +} + } // namespace cuopt::linear_programming::dual_simplex::test From de78eb6f790ad54cadf45b476f2654a522ef1a9a Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 22 May 2025 19:07:52 -0700 Subject: [PATCH 3/3] Fix for column scaling on problems where column has norm 0 --- cpp/src/dual_simplex/scaling.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/scaling.cpp b/cpp/src/dual_simplex/scaling.cpp index b158dadef..3a19ec613 100644 --- a/cpp/src/dual_simplex/scaling.cpp +++ b/cpp/src/dual_simplex/scaling.cpp @@ -48,7 +48,7 @@ i_t column_scaling(const lp_problem_t& unscaled, const f_t x = scaled.A.x[p]; sum += x * x; } - f_t col_norm_j = column_scaling[j] = std::sqrt(sum); + f_t col_norm_j = column_scaling[j] = sum > 0 ? std::sqrt(sum) : 1.0; max = std::max(col_norm_j, max); } settings.log.printf("Scaling matrix. Maximum column norm %e\n", max);