Skip to content

Bug fix: solution incorrect when simplex presolve removes empty columns #29

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
rapids-bot merged 4 commits into from
May 23, 2025
Merged

Bug fix: solution incorrect when simplex presolve removes empty columns #29

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f622513
Bug fix: solution incorrect when dual simplex presolve removes empty …
chris-maes May 22, 2025
769771f
Fix loop over removed variables. Add unit test.
chris-maes May 22, 2025
f7d243e
Merge branch 'branch-25.05-nvidia' into fix_presolve_removed_empty_cols
chris-maes May 22, 2025
de78eb6
Fix for column scaling on problems where column has norm 0
chris-maes May 23, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 59 additions & 4 deletions cpp/src/dual_simplex/presolve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,9 @@
namespace cuopt::linear_programming::dual_simplex {

template <typename i_t, typename f_t>
i_t remove_empty_cols(lp_problem_t<i_t, f_t>& problem, i_t& num_empty_cols)
i_t remove_empty_cols(lp_problem_t<i_t, f_t>& problem,
i_t& num_empty_cols,
presolve_info_t<i_t, f_t>& presolve_info)
{
constexpr bool verbose = false;
if (verbose) { printf("Removing %d empty columns\n", num_empty_cols); }
Expand All @@ -35,15 +37,22 @@ i_t remove_empty_cols(lp_problem_t<i_t, f_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<i_t> 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);
}
Expand All @@ -52,6 +61,7 @@ i_t remove_empty_cols(lp_problem_t<i_t, f_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
Expand All @@ -66,6 +76,7 @@ i_t remove_empty_cols(lp_problem_t<i_t, f_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--;
Expand Down Expand Up @@ -574,7 +585,8 @@ void convert_user_problem(const user_problem_t<i_t, f_t>& user_problem,
template <typename i_t, typename f_t>
i_t presolve(const lp_problem_t<i_t, f_t>& original,
const simplex_solver_settings_t<i_t, f_t>& settings,
lp_problem_t<i_t, f_t>& problem)
lp_problem_t<i_t, f_t>& problem,
presolve_info_t<i_t, f_t>& presolve_info)
{
problem = original;
std::vector<char> row_sense(problem.num_rows, '=');
Expand All @@ -595,6 +607,7 @@ i_t presolve(const lp_problem_t<i_t, f_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; }
}
Expand All @@ -606,7 +619,10 @@ i_t presolve(const lp_problem_t<i_t, f_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;
Expand Down Expand Up @@ -826,6 +842,38 @@ void uncrush_primal_solution(const user_problem_t<i_t, f_t>& user_problem,
std::copy(solution.begin(), solution.begin() + user_problem.num_cols, user_solution.data());
}

template <typename i_t, typename f_t>
void uncrush_solution(const presolve_info_t<i_t, f_t>& presolve_info,
const std::vector<f_t>& crushed_x,
const std::vector<f_t>& crushed_z,
std::vector<f_t>& uncrushed_x,
std::vector<f_t>& 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++;
}

k = 0;
for (const i_t j : presolve_info.removed_variables) {
uncrushed_x[j] = presolve_info.removed_values[k];
uncrushed_z[j] = presolve_info.removed_reduced_costs[k];
k++;
}
}

#ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE

template void convert_user_problem<int, double>(const user_problem_t<int, double>& user_problem,
Expand All @@ -841,7 +889,9 @@ template void convert_user_lp_with_guess<int, double>(

template int presolve<int, double>(const lp_problem_t<int, double>& original,
const simplex_solver_settings_t<int, double>& settings,
lp_problem_t<int, double>& presolved);
lp_problem_t<int, double>& presolved,
presolve_info_t<int, double>& presolve_info);

template void crush_primal_solution<int, double>(const user_problem_t<int, double>& user_problem,
const lp_problem_t<int, double>& problem,
const std::vector<double>& user_solution,
Expand All @@ -853,6 +903,11 @@ template void uncrush_primal_solution<int, double>(const user_problem_t<int, dou
const std::vector<double>& solution,
std::vector<double>& user_solution);

template void uncrush_solution<int, double>(const presolve_info_t<int, double>& presolve_info,
const std::vector<double>& crushed_x,
const std::vector<double>& crushed_z,
std::vector<double>& uncrushed_x,
std::vector<double>& uncrushed_z);
#endif

} // namespace cuopt::linear_programming::dual_simplex
22 changes: 21 additions & 1 deletion cpp/src/dual_simplex/presolve.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,18 @@ struct lp_problem_t {
f_t obj_scale; // 1.0 for min, -1.0 for max
};

template <typename i_t, typename f_t>
struct presolve_info_t {
// indices of variables in the original problem that remain in the presolved problem
std::vector<i_t> remaining_variables;
// indicies of variables in the original problem that have been removed in the presolved problem
std::vector<i_t> removed_variables;
// values of the removed variables
std::vector<f_t> removed_values;
// values of the removed reduced costs
std::vector<f_t> removed_reduced_costs;
};

template <typename i_t, typename f_t>
void convert_user_problem(const user_problem_t<i_t, f_t>& user_problem,
lp_problem_t<i_t, f_t>& problem,
Expand All @@ -70,7 +82,8 @@ void convert_user_lp_with_guess(const user_problem_t<i_t, f_t>& user_problem,
template <typename i_t, typename f_t>
i_t presolve(const lp_problem_t<i_t, f_t>& original,
const simplex_solver_settings_t<i_t, f_t>& settings,
lp_problem_t<i_t, f_t>& presolved);
lp_problem_t<i_t, f_t>& presolved,
presolve_info_t<i_t, f_t>& presolve_info);

template <typename i_t, typename f_t>
void crush_primal_solution(const user_problem_t<i_t, f_t>& user_problem,
Expand Down Expand Up @@ -102,4 +115,11 @@ void uncrush_primal_solution(const user_problem_t<i_t, f_t>& user_problem,
const std::vector<f_t>& solution,
std::vector<f_t>& user_solution);

template <typename i_t, typename f_t>
void uncrush_solution(const presolve_info_t<i_t, f_t>& presolve_info,
const std::vector<f_t>& crushed_x,
const std::vector<f_t>& crushed_z,
std::vector<f_t>& uncrushed_x,
std::vector<f_t>& uncrushed_z);

} // namespace cuopt::linear_programming::dual_simplex
24 changes: 23 additions & 1 deletion cpp/src/dual_simplex/scaling.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ i_t column_scaling(const lp_problem_t<i_t, f_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);
Expand Down Expand Up @@ -76,13 +76,35 @@ i_t column_scaling(const lp_problem_t<i_t, f_t>& unscaled,
return 0;
}

template <typename i_t, typename f_t>
void unscale_solution(const std::vector<f_t>& column_scaling,
const std::vector<f_t>& scaled_x,
const std::vector<f_t>& scaled_z,
std::vector<f_t>& unscaled_x,
std::vector<f_t>& 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<int, double>(const lp_problem_t<int, double>& unscaled,
const simplex_solver_settings_t<int, double>& settings,
lp_problem_t<int, double>& scaled,
std::vector<double>& column_scaling);

template void unscale_solution<int, double>(const std::vector<double>& column_scaling,
const std::vector<double>& scaled_x,
const std::vector<double>& scaled_z,
std::vector<double>& unscaled_x,
std::vector<double>& unscaled_z);

#endif

} // namespace cuopt::linear_programming::dual_simplex
7 changes: 7 additions & 0 deletions cpp/src/dual_simplex/scaling.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,4 +32,11 @@ i_t column_scaling(const lp_problem_t<i_t, f_t>& unscaled,
lp_problem_t<i_t, f_t>& scaled,
std::vector<f_t>& column_scaling);

template <typename i_t, typename f_t>
void unscale_solution(const std::vector<f_t>& column_scaling,
const std::vector<f_t>& scaled_x,
const std::vector<f_t>& scaled_z,
std::vector<f_t>& unscaled_x,
std::vector<f_t>& unscaled_z);

} // namespace cuopt::linear_programming::dual_simplex
13 changes: 7 additions & 6 deletions cpp/src/dual_simplex/solve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,8 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original
{
lp_status_t lp_status = lp_status_t::UNSET;
lp_problem_t<i_t, f_t> presolved_lp(1, 1, 1);
const i_t ok = presolve(original_lp, settings, presolved_lp);
presolve_info_t<i_t, f_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;
Expand Down Expand Up @@ -199,11 +200,11 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_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<f_t> unscaled_x(lp.num_cols);
std::vector<f_t> unscaled_z(lp.num_cols);
unscale_solution<i_t, f_t>(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;
Expand Down
87 changes: 87 additions & 0 deletions cpp/tests/dual_simplex/unit_tests/solve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<double> value({15, 100, 90, 0, 60, 40, 15, 10, 1});
std::vector<double> 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<int, double> 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<int, double> settings;

cuopt::linear_programming::dual_simplex::lp_solution_t<int, double> 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