Skip to content

Fix knapsack rational coefficient scaling #1754

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

Draft
aliceb-nv wants to merge 5 commits into
base: main
Choose a base branch
from
Draft

Fix knapsack rational coefficient scaling #1754

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
e190192
Reject non-integral rows in rational_coefficients
aliceb-nv Aug 20, 2026
5d796ba
Keep full precision weights in the knapsack separation DP
aliceb-nv Aug 20, 2026
1d545fe
Check the cover property before emitting a knapsack cut
aliceb-nv Aug 20, 2026
440e3a0
Add knapsack cover regression tests
aliceb-nv Aug 20, 2026
df84179
Drop the knapsack test-only entry points
aliceb-nv Aug 20, 2026
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
62 changes: 52 additions & 10 deletions cpp/src/cuts/cuts.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2406,6 +2406,7 @@ i_t knapsack_generation_t<i_t, f_t>::generate_knapsack_cut(
f_t objective_constant = 0.0;
std::vector<i_t> fixed_variables;
std::vector<f_t> fixed_values;
std::vector<f_t> fixed_weights;
const f_t x_tol = 1e-5;
for (i_t k = 0; k < knapsack_inequality.size(); k++) {
const i_t j = knapsack_inequality.index(k);
Expand All @@ -2417,6 +2418,7 @@ i_t knapsack_generation_t<i_t, f_t>::generate_knapsack_cut(
// if xstar_j is close to 0, then we can fix z to zero
fixed_variables.push_back(j);
fixed_values.push_back(0.0);
fixed_weights.push_back(knapsack_inequality.vector.x[k]);
seperation_rhs -= knapsack_inequality.vector.x[k];
// No need to adjust the objective constant
continue;
Expand All @@ -2425,6 +2427,7 @@ i_t knapsack_generation_t<i_t, f_t>::generate_knapsack_cut(
// if xstar_j is close to 1, then we can fix z to 1
fixed_variables.push_back(j);
fixed_values.push_back(1.0);
fixed_weights.push_back(knapsack_inequality.vector.x[k]);
// Note seperation rhs is unchanged
objective_constant += vj;
continue;
Expand Down Expand Up @@ -2466,12 +2469,28 @@ i_t knapsack_generation_t<i_t, f_t>::generate_knapsack_cut(
return -1;
}

i_t cover_size = 0;
i_t cover_size = 0;
f_t cover_weight = 0.0;
for (i_t k = 0; k < solution.size(); k++) {
if (solution[k] == 0.0) { cover_size++; }
if (solution[k] == 0.0) {
cover_size++;
cover_weight += weights[k];
}
}
for (i_t k = 0; k < fixed_values.size(); k++) {
if (fixed_values[k] == 1.0) { cover_size++; }
if (fixed_values[k] == 1.0) {
cover_size++;
cover_weight += fixed_weights[k];
}
}

// sum_{j in C} a_j > beta is what makes sum_{j in C} x_j <= |C| - 1 valid. The coefficients are
// integral here, so demand a full unit rather than letting rounding in the sums decide.
const bool is_cover = cover_weight >= knapsack_inequality.rhs + 1.0 - tol;
cuopt_assert(is_cover, "knapsack separation produced a set that is not a cover");
if (!is_cover) {
restore_complemented(complemented_variables);
return -1;
}

cut.reserve(cover_size);
Expand Down Expand Up @@ -2641,6 +2660,9 @@ void knapsack_generation_t<i_t, f_t>::minimal_cover_and_partition(
}
}

cuopt_assert(cover_sum >= beta + 1.0 - 1e-6,
"minimal cover reduction dropped an item the cover needed");

// Go through and correct cover_indicies and cover_coefficients
for (i_t k = 0; k < cover_coefficients.size();) {
if (cover_coefficients[k] == 0.0) {
Expand Down Expand Up @@ -2958,8 +2980,7 @@ f_t knapsack_generation_t<i_t, f_t>::solve_knapsack_problem(const std::vector<f_
}
}

i_t sum_value = std::accumulate(scaled_values.begin(), scaled_values.end(), 0);
const i_t INT_INF = std::numeric_limits<i_t>::max() / 2;
i_t sum_value = std::accumulate(scaled_values.begin(), scaled_values.end(), 0);
if (verbose) { settings_.log.printf("sum value %d\n", sum_value); }
const i_t max_size = 10000;
if (sum_value <= 0.0 || sum_value >= max_size) {
Expand All @@ -2972,10 +2993,12 @@ f_t knapsack_generation_t<i_t, f_t>::solve_knapsack_problem(const std::vector<f_

solution.assign(n, 0.0);

// dp(j, v) = minimum weight using first j items to get value v
dense_matrix_t<i_t, i_t> dp(n + 1, sum_value + 1, INT_INF);
// dp(j, v) = minimum weight using first j items to get value v.
// The weights are carried at full precision: rounding one down would let the DP return a set

@chris-maes chris-maes Aug 20, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure I understand why we need to use floats here. Where is the "rounding down" occurring that is causing the issue? And why does floats fix it?

@chris-maes chris-maes Aug 20, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah ok. Reading further, it sounds like we tried to convert an inequality with floating point coefficients into a knapsack constraint with integer coefficients, but when we do this we can make some small error if we round down. Then when we do dynamic programming we can think we are under the weight limit, when we are actually over.

But I'm not sure the fix should be use floating point in the DP table. That can make the table a lot bigger. Maybe we should reject inequalities with large errors when converting to integer coefficients in our classification of knapsack constraints, so we don't arrive at this case in the first place? Or if this occurs when the error is quite small, maybe we just need to do a check at the end to see if the selected weights violate the capacity (similar to what is done in the assert now)

@chris-maes chris-maes Aug 20, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or just use ceil(weights[j-1]) instead of floor(weights[j-1])

@mlubin mlubin Aug 21, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah I think the typical solution here is careful rounding.

// that violates the capacity, and the caller reads the complement of that set as a cover.
dense_matrix_t<i_t, f_t> dp(n + 1, sum_value + 1, inf);
dense_matrix_t<i_t, uint8_t> take(n + 1, sum_value + 1, 0);
dp(0, 0) = 0;
dp(0, 0) = 0.0;

// 4. Dynamic programming
for (i_t j = 1; j <= n; ++j) {
Expand All @@ -2985,8 +3008,7 @@ f_t knapsack_generation_t<i_t, f_t>::solve_knapsack_problem(const std::vector<f_

// Take item j-1 if possible
if (v >= scaled_values[j - 1]) {
i_t candidate =
dp(j - 1, v - scaled_values[j - 1]) + static_cast<i_t>(std::floor(weights[j - 1]));

@chris-maes chris-maes Aug 20, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems like the floor is the big issue here. Because it causes us to underestimate the weight. Could the issue be solved by just using std::ceil(weights[j-1]) here? If weights[k] is an integer than ceil(weights[k]) == floor(weights[k]). But if it isn't exact, we will round up, overestimate how much weight we can carry, and thus still produce a valid solution.

f_t candidate = dp(j - 1, v - scaled_values[j - 1]) + weights[j - 1];
if (candidate < dp(j, v)) {
dp(j, v) = candidate;
take(j, v) = 1;
Expand All @@ -3012,6 +3034,15 @@ f_t knapsack_generation_t<i_t, f_t>::solve_knapsack_problem(const std::vector<f_
}
}

#ifdef ASSERT_MODE
f_t selected_weight = 0.0;
for (i_t j = 0; j < n; ++j) {
selected_weight += solution[j] * weights[j];
}
cuopt_assert(selected_weight <= rhs + settings_.primal_tol,
"knapsack dynamic program returned a solution over capacity");
#endif

objective = best_value * scale;
return objective;
}
Expand Down Expand Up @@ -4629,6 +4660,17 @@ bool rational_coefficients(const std::vector<variable_type_t>& var_types,

rational_inequality.scale(scalar);

// The scaled product can land an ulp off the integer it represents. Callers rely on the

@chris-maes chris-maes Aug 20, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This code is a bit strange. First, the integral_tol used here, doesn't match our default integer_tol. It would probably be best to pass settings in and use settings.integer_tol. Unless, this should be a separate tolerance.

Second, I think this function is used by more than knapsack cuts. So it's a bit weird to reference knapsack cuts in this code. I think this code is just trying to generate an inequality with rational coefficients, whereas knapsack needs an inequality with integer coefficients. Maybe this code should be moved into a post-processing step, say in a function called integer_coefficients (that first calls rationalize_coefficients)

Finally, rounding introduces some error. That just multiplying by a scalar does not. So, we probably don't want to do this in all cases. Maybe we also want to track the total amount of error introduced (instead of just tracking the error on individual coefficients).

// integer-variable coefficients being exact integers: the knapsack cover test
// sum_C a_j > beta is only equivalent to sum_C a_j >= beta + 1 for integral a_j.
constexpr f_t integral_tol = 1e-6;
for (i_t k : indices) {
const f_t scaled = rational_inequality.vector.x[k];
const f_t rounded = std::round(scaled);
if (std::abs(scaled - rounded) > integral_tol) { return false; }
rational_inequality.vector.x[k] = rounded;
}

return true;
}

Expand Down
69 changes: 69 additions & 0 deletions cpp/tests/mip/cuts_test.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -442,6 +442,19 @@ void disable_all_cuts(mip_solver_settings_t<int, double>& settings)
settings.strong_chvatal_gomory_cuts = 0;
}

void disable_non_knapsack_cuts(mip_solver_settings_t<int, double>& settings)
{
settings.max_cut_passes = 10;
settings.knapsack_cuts = 1;
settings.clique_cuts = 0;
settings.zero_half_cuts = 0;
settings.mixed_integer_gomory_cuts = 0;
settings.mir_cuts = 0;
settings.strong_chvatal_gomory_cuts = 0;
settings.flow_cover_cuts = 0;
settings.implied_bound_cuts = 0;
}

bool cut_is_invalid_for_incumbent(const std::vector<int>& cut_vars,
const std::vector<double>& incumbent,
double tol)
Expand Down Expand Up @@ -934,6 +947,62 @@ TEST(cuts, test_cuts_2)
EXPECT_EQ(solution.get_num_nodes(), 0);
}

io::mps_data_model_t<int, double> create_knapsack_cover_floor_problem()
{
// The odd cycle over z1, z2, z3 makes z = (0.5, 0.5, 0.5), w = 0 the unique LP optimum, which
// puts y at 0.412078. Integrality moves the optimum to w = 1 with y = 0 and objective 3.
//
// The capacity coefficients are load bearing: scaling that row to integers multiplies by 100,

@chris-maes chris-maes Aug 20, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nit: "load bearing" sounds like an AI generated phrase. Could we reword?

// and 135.45 and 135.42 do not land on integers when scaled, which is what the knapsack
// separator needs them to do.
return cuopt::test::parse_inline_lp(R"LP(
Minimize
obj: 2 y + z1 + z2 + z3 + 3 w
Subject To
capacity: -450 y + 135.45 z1 + 135.42 z2 + 100 z3 <= 0
tri12: z1 + z2 + w >= 1
tri13: z1 + z3 + w >= 1
tri23: z2 + z3 + w >= 1
Binaries
y
z1
z2
z3
w
End
)LP");
}

TEST(cuts, knapsack_cover_floor_regression)
{
const raft::handle_t handle_{};
auto problem = create_knapsack_cover_floor_problem();

mip_solver_settings_t<int, double> settings;
settings.time_limit = 10.;
disable_non_knapsack_cuts(settings);
settings.presolver = presolver_t::None;

mip_solution_t<int, double> solution = solve_mip(&handle_, problem, settings);
EXPECT_EQ(solution.get_termination_status(), mip_termination_status_t::Optimal);
EXPECT_NEAR(3.0, solution.get_objective_value(), 1e-6);
}

TEST(cuts, knapsack_cover_floor_regression_reference)
{
const raft::handle_t handle_{};
auto problem = create_knapsack_cover_floor_problem();

mip_solver_settings_t<int, double> settings;
settings.time_limit = 10.;
disable_all_cuts(settings);
settings.presolver = presolver_t::None;

mip_solution_t<int, double> solution = solve_mip(&handle_, problem, settings);
EXPECT_EQ(solution.get_termination_status(), mip_termination_status_t::Optimal);
EXPECT_NEAR(3.0, solution.get_objective_value(), 1e-6);
}

TEST(cuts, test_duplicate_cuts_detection)
{
simplex::simplex_solver_settings_t<int, double> settings;
Expand Down