From 62d7a5f186fef7d0c0de95ae3138006429f7aba6 Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Mon, 16 Mar 2026 14:42:27 +0100 Subject: [PATCH 01/19] pass global upper bound to strong branching. Added an estimation for the cutoff when no incumbent was found. --- cpp/src/branch_and_bound/branch_and_bound.cpp | 2 +- cpp/src/branch_and_bound/pseudo_costs.cpp | 32 ++++++++++++++++++- cpp/src/branch_and_bound/pseudo_costs.hpp | 1 + 3 files changed, 33 insertions(+), 2 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index ba2b63983e..5c7a5092be 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1175,7 +1175,6 @@ std::pair branch_and_bound_t::upd dual::status_t lp_status, Policy& policy) { - constexpr f_t inf = std::numeric_limits::infinity(); const f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; lp_problem_t& leaf_problem = worker->leaf_problem; lp_solution_t& leaf_solution = worker->leaf_solution; @@ -2456,6 +2455,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_soln_.x, fractional, root_objective_, + upper_bound_, root_vstatus_, edge_norms_, pc_); diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index ee7e2f7803..66115cac18 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -31,6 +31,7 @@ void strong_branch_helper(i_t start, const std::vector& var_types, const std::vector& fractional, f_t root_obj, + f_t upper_bound, const std::vector& root_soln, const std::vector& root_vstatus, const std::vector& edge_norms, @@ -62,6 +63,20 @@ void strong_branch_helper(i_t start, if (elapsed_time > settings.time_limit) { break; } child_settings.time_limit = std::max(0.0, settings.time_limit - elapsed_time); child_settings.iteration_limit = 200; + + if (std::isfinite(upper_bound)) { + child_settings.cut_off = upper_bound + settings.dual_tol; + } else { + child_settings.cut_off = 0; + for (i_t i = 0; i < original_lp.num_cols; ++i) { + if (original_lp.objective[i] < 0) { + child_settings.cut_off += original_lp.objective[i] * child_problem.upper[i]; + } else if (original_lp.objective[i] > 0) { + child_settings.cut_off += original_lp.objective[i] * child_problem.lower[i]; + } + } + } + lp_solution_t solution(original_lp.num_rows, original_lp.num_cols); i_t iter = 0; std::vector vstatus = root_vstatus; @@ -168,10 +183,22 @@ f_t trial_branching(const lp_problem_t& original_lp, i_t lp_iter_upper = upper_max_lp_iter; i_t lp_iter_lower = lower_max_lp_iter; child_settings.iteration_limit = std::clamp(bnb_lp_iter_per_node, lp_iter_lower, lp_iter_upper); - child_settings.cut_off = upper_bound + settings.dual_tol; child_settings.inside_mip = 2; child_settings.scale_columns = false; + if (std::isfinite(upper_bound)) { + child_settings.cut_off = upper_bound + settings.dual_tol; + } else { + child_settings.cut_off = 0; + for (i_t i = 0; i < original_lp.num_cols; ++i) { + if (original_lp.objective[i] < 0) { + child_settings.cut_off += original_lp.objective[i] * child_problem.upper[i]; + } else if (original_lp.objective[i] > 0) { + child_settings.cut_off += original_lp.objective[i] * child_problem.lower[i]; + } + } + } + lp_solution_t solution(original_lp.num_rows, original_lp.num_cols); i_t iter = 0; std::vector child_vstatus = vstatus; @@ -305,6 +332,7 @@ void strong_branching(const user_problem_t& original_problem, const std::vector root_soln, const std::vector& fractional, f_t root_obj, + f_t upper_bound, const std::vector& root_vstatus, const std::vector& edge_norms, pseudo_costs_t& pc) @@ -429,6 +457,7 @@ void strong_branching(const user_problem_t& original_problem, var_types, fractional, root_obj, + upper_bound, root_soln, root_vstatus, edge_norms, @@ -784,6 +813,7 @@ template void strong_branching(const user_problem_t& o const std::vector root_soln, const std::vector& fractional, double root_obj, + double upper_bound, const std::vector& root_vstatus, const std::vector& edge_norms, pseudo_costs_t& pc); diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index 6b6c6917b6..ebcbd4abad 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.hpp +++ b/cpp/src/branch_and_bound/pseudo_costs.hpp @@ -525,6 +525,7 @@ void strong_branching(const user_problem_t& original_problem, const std::vector root_soln, const std::vector& fractional, f_t root_obj, + f_t upper_bound, const std::vector& root_vstatus, const std::vector& edge_norms, pseudo_costs_t& pc); From aa4d55c5c9ef7b8ef78d88fea36e381c4efb6f58 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Mon, 16 Mar 2026 16:59:00 -0700 Subject: [PATCH 02/19] Use a single dual simplex pivot to initialize pseudocost --- cpp/src/branch_and_bound/branch_and_bound.cpp | 5 + cpp/src/branch_and_bound/pseudo_costs.cpp | 286 +++++++++++++++- cpp/src/branch_and_bound/pseudo_costs.hpp | 7 +- cpp/src/dual_simplex/phase2.cpp | 306 ++++++++++-------- cpp/src/dual_simplex/phase2.hpp | 23 ++ cpp/src/math_optimization/solver_settings.cu | 2 +- 6 files changed, 481 insertions(+), 148 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index ba2b63983e..dac6c82a24 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2454,10 +2454,15 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut exploration_stats_.start_time, var_types_, root_relax_soln_.x, + root_relax_soln_.y, + root_relax_soln_.z, fractional, root_objective_, root_vstatus_, edge_norms_, + basic_list, + nonbasic_list, + basis_update, pc_); } diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index ee7e2f7803..d86c06740e 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -22,6 +22,256 @@ namespace cuopt::linear_programming::dual_simplex { namespace { + +template +struct objective_estimate_t { + f_t down_obj; + f_t up_obj; +}; + +template +f_t compute_step_length(const simplex_solver_settings_t& settings, + const std::vector& vstatus, + const std::vector& z, + const std::vector& delta_z, + const std::vector& delta_z_indices) +{ + f_t step_length = inf; + f_t pivot_tol = settings.pivot_tol; + const i_t nz = delta_z_indices.size(); + for (i_t h = 0; h < nz; h++) { + const i_t j = delta_z_indices[h]; + if (vstatus[j] == variable_status_t::NONBASIC_LOWER && delta_z[j] < -pivot_tol) { + const f_t ratio = -z[j] / delta_z[j]; + if (ratio < step_length) { + step_length = ratio; + } + } else if (vstatus[j] == variable_status_t::NONBASIC_UPPER && delta_z[j] > pivot_tol) { + const f_t ratio = -z[j] / delta_z[j]; + if (ratio < step_length) { + step_length = ratio; + } + } + } + return step_length; +} + + +template +objective_estimate_t single_pivot_objective_estimate( + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const csc_matrix_t& A_transpose, + const std::vector& vstatus, + i_t variable_j, + i_t basic_j, + f_t root_obj, + const std::vector& x, + const std::vector& y, + const std::vector& z, + const std::vector& basic_list, + const std::vector& nonbasic_list, + const std::vector& nonbasic_mark, + basis_update_mpf_t& basis_factors, + std::vector& workspace, + std::vector& delta_z, + f_t& work_estimate) +{ + // Compute the objective estimate for the down and up branches of variable j + + // Down branch + i_t direction = -1; + sparse_vector_t e_k(lp.num_rows, 0); + e_k.i.push_back(basic_j); + e_k.x.push_back(-f_t(direction)); + + sparse_vector_t delta_y(lp.num_rows, 0); + basis_factors.b_transpose_solve(e_k, delta_y); + + // Compute delta_z_N = -N^T * delta_y + i_t delta_y_nz0 = 0; + const i_t nz_delta_y = delta_y.i.size(); + for (i_t k = 0; k < nz_delta_y; k++) { + if (std::abs(delta_y.x[k]) > 1e-12) { delta_y_nz0++; } + } + work_estimate += nz_delta_y; + const f_t delta_y_nz_percentage = delta_y_nz0 / static_cast(lp.num_rows) * 100.0; + const bool use_transpose = delta_y_nz_percentage <= 30.0; + std::vector delta_z_indices; + // delta_z starts out all zero + if (use_transpose) { + compute_delta_z(A_transpose, + delta_y, + variable_j, + -1, + nonbasic_mark, + workspace, + delta_z_indices, + delta_z, + work_estimate); + } else { + std::vector delta_y_dense(lp.num_rows, 0); + delta_y.to_dense(delta_y_dense); + compute_reduced_cost_update(lp, + basic_list, + nonbasic_list, + delta_y_dense, + variable_j, + -1, + workspace, + delta_z_indices, + delta_z, + work_estimate); + } + + + + // Verify dual feasibility +#ifdef CHECK_DUAL_FEASIBILITY + { + std::vector dual_residual = z; + for (i_t j = 0; j < lp.num_cols; j++) { + dual_residual[j] -= lp.objective[j]; + } + matrix_transpose_vector_multiply(lp.A, 1.0, y, 1.0, dual_residual); + f_t dual_residual_norm = vector_norm_inf(dual_residual); + settings.log.printf("Dual residual norm: %e\n", dual_residual_norm); + } +#endif + + // Compute the step-length + f_t step_length = compute_step_length(settings, vstatus, z, delta_z, delta_z_indices); + + // Handle the leaving variable case + + f_t delta_obj = step_length * (x[variable_j] - std::floor(x[variable_j])); +#ifdef CHECK_DELTA_OBJ + f_t delta_obj_check = 0.0; + for (i_t k = 0; k < delta_y.i.size(); k++) { + delta_obj_check += lp.rhs[delta_y.i[k]] * delta_y.x[k]; + } + for (i_t h = 0; h < delta_z_indices.size(); h++) { + const i_t j = delta_z_indices[h]; + if (vstatus[j] == variable_status_t::NONBASIC_LOWER) { + delta_obj_check += lp.lower[j] * delta_z[j]; + } else if (vstatus[j] == variable_status_t::NONBASIC_UPPER) { + delta_obj_check += lp.upper[j] * delta_z[j]; + } + } + delta_obj_check += std::floor(x[variable_j]) * delta_z[variable_j]; + delta_obj_check *= step_length; + if (std::abs(delta_obj_check - delta_obj) > 1e-6) { + settings.log.printf("Delta obj check %e. Delta obj %e. Step length %e.\n", delta_obj_check, delta_obj, step_length); + } +#endif + + f_t down_obj = root_obj + delta_obj; + settings.log.printf("Down branch %d. Step length: %e. Objective estimate: %e. Delta obj: %e. Root obj: %e.\n", variable_j, step_length, down_obj, delta_obj, root_obj); + + // Up branch + direction = 1; + // Negate delta_z + for (i_t j : delta_z_indices) + { + delta_z[j] *= -1.0; + } + + // Compute the step-length + step_length = compute_step_length(settings, vstatus, z, delta_z, delta_z_indices); + + delta_obj = step_length * (std::ceil(x[variable_j]) - x[variable_j]); + f_t up_obj = root_obj + delta_obj; + settings.log.printf("Up branch %d. Step length: %e. Objective estimate: %e. Delta obj: %e. Root obj: %e.\n", variable_j, step_length, up_obj, delta_obj, root_obj); + + delta_z_indices.push_back(variable_j); + + // Clear delta_z + for (i_t j : delta_z_indices) + { + delta_z[j] = 0.0; + workspace[j] = 0; + } + +#ifdef CHECK_DELTA_Z + for (i_t j = 0; j < lp.num_cols; j++) { + if (delta_z[j] != 0.0) { + settings.log.printf("Delta z %d: %e\n", j, delta_z[j]); + } + } + for (i_t j = 0; j < lp.num_cols; j++) { + if (workspace[j] != 0) { + settings.log.printf("Workspace %d: %d\n", j, workspace[j]); + } + } +#endif + + + return objective_estimate_t{.down_obj = down_obj, .up_obj = up_obj}; +} + + +template +void initialize_pseudo_costs_with_estimate(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& vstatus, + const std::vector& x, + const std::vector& y, + const std::vector& z, + const std::vector& basic_list, + const std::vector& nonbasic_list, + const std::vector& fractional, + f_t root_obj, + basis_update_mpf_t& basis_factors, + pseudo_costs_t& pc) +{ + i_t m = lp.num_rows; + i_t n = lp.num_cols; + + csc_matrix_t A_transpose(m, n, 0); + lp.A.transpose(A_transpose); + + std::vector delta_z(n, 0); + std::vector workspace(n, 0); + + f_t work_estimate = 0; + + std::vector basic_map(n, -1); + for (i_t i = 0; i < m; i++) { + basic_map[basic_list[i]] = i; + } + + std::vector nonbasic_mark(n, -1); + for (i_t i = 0; i < n - m; i++) { + nonbasic_mark[nonbasic_list[i]] = i; + } + + for (i_t k = 0; k < fractional.size(); k++) { + const i_t j = fractional[k]; + objective_estimate_t estimate = single_pivot_objective_estimate(lp, + settings, + A_transpose, + vstatus, + j, + basic_map[j], + root_obj, + x, + y, + z, + basic_list, + nonbasic_list, + nonbasic_mark, + basis_factors, + workspace, + delta_z, + work_estimate); + pc.strong_branch_down[k] = estimate.down_obj; + pc.strong_branch_up[k] = estimate.up_obj; + } +} + + + + template void strong_branch_helper(i_t start, i_t end, @@ -302,11 +552,16 @@ void strong_branching(const user_problem_t& original_problem, const simplex_solver_settings_t& settings, f_t start_time, const std::vector& var_types, - const std::vector root_soln, + const std::vector& root_x, + const std::vector& root_y, + const std::vector& root_z, const std::vector& fractional, f_t root_obj, const std::vector& root_vstatus, const std::vector& edge_norms, + const std::vector& basic_list, + const std::vector& nonbasic_list, + basis_update_mpf_t& basis_factors, pseudo_costs_t& pc) { pc.resize(original_lp.num_cols); @@ -317,7 +572,7 @@ void strong_branching(const user_problem_t& original_problem, const f_t elapsed_time = toc(start_time); if (elapsed_time > settings.time_limit) { return; } - if (settings.mip_batch_pdlp_strong_branching) { + if (settings.mip_batch_pdlp_strong_branching == 1) { settings.log.printf("Batch PDLP strong branching enabled\n"); f_t start_batch = tic(); @@ -328,7 +583,7 @@ void strong_branching(const user_problem_t& original_problem, // Convert the root_soln to the original problem space std::vector original_root_soln_x; - uncrush_primal_solution(original_problem, original_lp, root_soln, original_root_soln_x); + uncrush_primal_solution(original_problem, original_lp, root_x, original_root_soln_x); std::vector fraction_values; @@ -394,6 +649,20 @@ void strong_branching(const user_problem_t& original_problem, pc.strong_branch_down[k] = obj_down - root_obj; pc.strong_branch_up[k] = obj_up - root_obj; } + } else if (settings.mip_batch_pdlp_strong_branching == 2) { + initialize_pseudo_costs_with_estimate(original_lp, + settings, + root_vstatus, + root_x, + root_y, + root_z, + basic_list, + nonbasic_list, + fractional, + root_obj, + basis_factors, + pc); + } else { settings.log.printf("Strong branching using %d threads and %ld fractional variables\n", settings.num_threads, @@ -429,7 +698,7 @@ void strong_branching(const user_problem_t& original_problem, var_types, fractional, root_obj, - root_soln, + root_x, root_vstatus, edge_norms, pc); @@ -438,7 +707,7 @@ void strong_branching(const user_problem_t& original_problem, settings.log.printf("Strong branching completed in %.2fs\n", toc(strong_branching_start_time)); } - pc.update_pseudo_costs_from_strong_branching(fractional, root_soln); + pc.update_pseudo_costs_from_strong_branching(fractional, root_x); } template @@ -781,11 +1050,16 @@ template void strong_branching(const user_problem_t& o const simplex_solver_settings_t& settings, double start_time, const std::vector& var_types, - const std::vector root_soln, + const std::vector& root_x, + const std::vector& root_y, + const std::vector& root_z, const std::vector& fractional, double root_obj, const std::vector& root_vstatus, const std::vector& edge_norms, + const std::vector& basic_list, + const std::vector& nonbasic_list, + basis_update_mpf_t& basis_factors, pseudo_costs_t& pc); #endif diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index 6b6c6917b6..e34dc5b6a4 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.hpp +++ b/cpp/src/branch_and_bound/pseudo_costs.hpp @@ -522,11 +522,16 @@ void strong_branching(const user_problem_t& original_problem, const simplex_solver_settings_t& settings, f_t start_time, const std::vector& var_types, - const std::vector root_soln, + const std::vector& root_x, + const std::vector& root_y, + const std::vector& root_z, const std::vector& fractional, f_t root_obj, const std::vector& root_vstatus, const std::vector& edge_norms, + const std::vector& basic_list, + const std::vector& nonbasic_list, + basis_update_mpf_t& basis_factors, pseudo_costs_t& pc); } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 43429ba2de..eca9dbd251 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -82,6 +82,130 @@ class nvtx_range_guard { bool active_; }; + +template +void compute_reduced_cost_update(const lp_problem_t& lp, + const std::vector& basic_list, + const std::vector& nonbasic_list, + const std::vector& delta_y, + i_t leaving_index, + i_t direction, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + f_t& work_estimate) +{ + const i_t m = lp.num_rows; + const i_t n = lp.num_cols; + + size_t nnzs_processed = 0; + // delta_zB = sigma*ei + for (i_t k = 0; k < m; k++) { + const i_t j = basic_list[k]; + delta_z[j] = 0; + } + work_estimate += 2 * m; + delta_z[leaving_index] = direction; + // delta_zN = -N'*delta_y + const i_t num_nonbasic = n - m; + for (i_t k = 0; k < num_nonbasic; k++) { + const i_t j = nonbasic_list[k]; + // z_j <- -A(:, j)'*delta_y + const i_t col_start = lp.A.col_start[j]; + const i_t col_end = lp.A.col_start[j + 1]; + f_t dot = 0.0; + for (i_t p = col_start; p < col_end; ++p) { + dot += lp.A.x[p] * delta_y[lp.A.i[p]]; + } + nnzs_processed += col_end - col_start; + + delta_z[j] = -dot; + if (dot != 0.0) { + delta_z_indices.push_back(j); // Note delta_z_indices has n elements reserved + delta_z_mark[j] = 1; + } + } + work_estimate += 3 * num_nonbasic; + work_estimate += 3 * nnzs_processed; + work_estimate += 2 * delta_z_indices.size(); +} + +template +void compute_delta_z(const csc_matrix_t& A_transpose, + const sparse_vector_t& delta_y, + i_t leaving_index, + i_t direction, + const std::vector& nonbasic_mark, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + f_t& work_estimate) +{ + // delta_zN = - N'*delta_y + const i_t nz_delta_y = delta_y.i.size(); + size_t nnz_processed = 0; + size_t nonbasic_marked = 0; + for (i_t k = 0; k < nz_delta_y; k++) { + const i_t i = delta_y.i[k]; + const f_t delta_y_i = delta_y.x[k]; + if (std::abs(delta_y_i) < 1e-12) { continue; } + const i_t row_start = A_transpose.col_start[i]; + const i_t row_end = A_transpose.col_start[i + 1]; + nnz_processed += row_end - row_start; + for (i_t p = row_start; p < row_end; ++p) { + const i_t j = A_transpose.i[p]; + if (nonbasic_mark[j] >= 0) { + delta_z[j] -= delta_y_i * A_transpose.x[p]; + nonbasic_marked++; + if (!delta_z_mark[j]) { + delta_z_mark[j] = 1; + delta_z_indices.push_back(j); + } + } + } + } + work_estimate += 4 * nz_delta_y; + work_estimate += 2 * nnz_processed; + work_estimate += 3 * nonbasic_marked; + work_estimate += 2 * delta_z_indices.size(); + + // delta_zB = sigma*ei + delta_z[leaving_index] = direction; + +#ifdef CHECK_CHANGE_IN_REDUCED_COST + const i_t m = A_transpose.n; + const i_t n = A_transpose.m; + std::vector delta_y_dense(m); + delta_y.to_dense(delta_y_dense); + std::vector delta_z_check(n); + std::vector delta_z_mark_check(n, 0); + std::vector delta_z_indices_check; + phase2::compute_reduced_cost_update(lp, + basic_list, + nonbasic_list, + delta_y_dense, + leaving_index, + direction, + delta_z_mark_check, + delta_z_indices_check, + delta_z_check, + work_estimate); + f_t error_check = 0.0; + for (i_t k = 0; k < n; ++k) { + const f_t diff = std::abs(delta_z[k] - delta_z_check[k]); + if (diff > 1e-6) { + printf("delta_z error %d transpose %e no transpose %e diff %e\n", + k, + delta_z[k], + delta_z_check[k], + diff); + } + error_check = std::max(error_check, diff); + } + if (error_check > 1e-6) { printf("delta_z error %e\n", error_check); } +#endif +} + namespace phase2 { // Computes vectors farkas_y, farkas_zl, farkas_zu that satisfy @@ -322,128 +446,7 @@ void initial_perturbation(const lp_problem_t& lp, n); } -template -void compute_reduced_cost_update(const lp_problem_t& lp, - const std::vector& basic_list, - const std::vector& nonbasic_list, - const std::vector& delta_y, - i_t leaving_index, - i_t direction, - std::vector& delta_z_mark, - std::vector& delta_z_indices, - std::vector& delta_z, - f_t& work_estimate) -{ - const i_t m = lp.num_rows; - const i_t n = lp.num_cols; - - size_t nnzs_processed = 0; - // delta_zB = sigma*ei - for (i_t k = 0; k < m; k++) { - const i_t j = basic_list[k]; - delta_z[j] = 0; - } - work_estimate += 2 * m; - delta_z[leaving_index] = direction; - // delta_zN = -N'*delta_y - const i_t num_nonbasic = n - m; - for (i_t k = 0; k < num_nonbasic; k++) { - const i_t j = nonbasic_list[k]; - // z_j <- -A(:, j)'*delta_y - const i_t col_start = lp.A.col_start[j]; - const i_t col_end = lp.A.col_start[j + 1]; - f_t dot = 0.0; - for (i_t p = col_start; p < col_end; ++p) { - dot += lp.A.x[p] * delta_y[lp.A.i[p]]; - } - nnzs_processed += col_end - col_start; - - delta_z[j] = -dot; - if (dot != 0.0) { - delta_z_indices.push_back(j); // Note delta_z_indices has n elements reserved - delta_z_mark[j] = 1; - } - } - work_estimate += 3 * num_nonbasic; - work_estimate += 3 * nnzs_processed; - work_estimate += 2 * delta_z_indices.size(); -} - -template -void compute_delta_z(const csc_matrix_t& A_transpose, - const sparse_vector_t& delta_y, - i_t leaving_index, - i_t direction, - std::vector& nonbasic_mark, - std::vector& delta_z_mark, - std::vector& delta_z_indices, - std::vector& delta_z, - f_t& work_estimate) -{ - // delta_zN = - N'*delta_y - const i_t nz_delta_y = delta_y.i.size(); - size_t nnz_processed = 0; - size_t nonbasic_marked = 0; - for (i_t k = 0; k < nz_delta_y; k++) { - const i_t i = delta_y.i[k]; - const f_t delta_y_i = delta_y.x[k]; - if (std::abs(delta_y_i) < 1e-12) { continue; } - const i_t row_start = A_transpose.col_start[i]; - const i_t row_end = A_transpose.col_start[i + 1]; - nnz_processed += row_end - row_start; - for (i_t p = row_start; p < row_end; ++p) { - const i_t j = A_transpose.i[p]; - if (nonbasic_mark[j] >= 0) { - delta_z[j] -= delta_y_i * A_transpose.x[p]; - nonbasic_marked++; - if (!delta_z_mark[j]) { - delta_z_mark[j] = 1; - delta_z_indices.push_back(j); - } - } - } - } - work_estimate += 4 * nz_delta_y; - work_estimate += 2 * nnz_processed; - work_estimate += 3 * nonbasic_marked; - work_estimate += 2 * delta_z_indices.size(); - // delta_zB = sigma*ei - delta_z[leaving_index] = direction; - -#ifdef CHECK_CHANGE_IN_REDUCED_COST - const i_t m = A_transpose.n; - const i_t n = A_transpose.m; - std::vector delta_y_dense(m); - delta_y.to_dense(delta_y_dense); - std::vector delta_z_check(n); - std::vector delta_z_mark_check(n, 0); - std::vector delta_z_indices_check; - phase2::compute_reduced_cost_update(lp, - basic_list, - nonbasic_list, - delta_y_dense, - leaving_index, - direction, - delta_z_mark_check, - delta_z_indices_check, - delta_z_check, - work_estimate); - f_t error_check = 0.0; - for (i_t k = 0; k < n; ++k) { - const f_t diff = std::abs(delta_z[k] - delta_z_check[k]); - if (diff > 1e-6) { - printf("delta_z error %d transpose %e no transpose %e diff %e\n", - k, - delta_z[k], - delta_z_check[k], - diff); - } - error_check = std::max(error_check, diff); - } - if (error_check > 1e-6) { printf("delta_z error %e\n", error_check); } -#endif -} template void compute_reduced_costs(const std::vector& objective, @@ -2932,30 +2935,30 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, PHASE2_NVTX_RANGE("DualSimplex::delta_z"); if (use_transpose) { sparse_delta_z++; - phase2::compute_delta_z(A_transpose, - delta_y_sparse, - leaving_index, - direction, - nonbasic_mark, - delta_z_mark, - delta_z_indices, - delta_z, - phase2_work_estimate); + compute_delta_z(A_transpose, + delta_y_sparse, + leaving_index, + direction, + nonbasic_mark, + delta_z_mark, + delta_z_indices, + delta_z, + phase2_work_estimate); } else { dense_delta_z++; // delta_zB = sigma*ei delta_y_sparse.to_dense(delta_y); phase2_work_estimate += delta_y.size(); - phase2::compute_reduced_cost_update(lp, - basic_list, - nonbasic_list, - delta_y, - leaving_index, - direction, - delta_z_mark, - delta_z_indices, - delta_z, - phase2_work_estimate); + compute_reduced_cost_update(lp, + basic_list, + nonbasic_list, + delta_y, + leaving_index, + direction, + delta_z_mark, + delta_z_indices, + delta_z, + phase2_work_estimate); } } timers.delta_z_time += timers.stop_timer(); @@ -3634,6 +3637,29 @@ template dual::status_t dual_phase2_with_advanced_basis( int& iter, std::vector& steepest_edge_norms, work_limit_context_t* work_unit_context); + +template +void compute_reduced_cost_update(const lp_problem_t& lp, + const std::vector& basic_list, + const std::vector& nonbasic_list, + const std::vector& delta_y, + int leaving_index, + int direction, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + double& work_estimate); + +template +void compute_delta_z(const csc_matrix_t& A_transpose, + const sparse_vector_t& delta_y, + int leaving_index, + int direction, + const std::vector& nonbasic_mark, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + double& work_estimate); #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/phase2.hpp b/cpp/src/dual_simplex/phase2.hpp index 7f01eb3cf7..5db797449c 100644 --- a/cpp/src/dual_simplex/phase2.hpp +++ b/cpp/src/dual_simplex/phase2.hpp @@ -81,4 +81,27 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, std::vector& delta_y_steepest_edge, work_limit_context_t* work_unit_context = nullptr); +template +void compute_reduced_cost_update(const lp_problem_t& lp, + const std::vector& basic_list, + const std::vector& nonbasic_list, + const std::vector& delta_y, + i_t leaving_index, + i_t direction, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + f_t& work_estimate); + +template +void compute_delta_z(const csc_matrix_t& A_transpose, + const sparse_vector_t& delta_y, + i_t leaving_index, + i_t direction, + const std::vector& nonbasic_mark, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + f_t& work_estimate); + } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index a60d508fac..13a8f74145 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -99,7 +99,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, {CUOPT_NUM_GPUS, &pdlp_settings.num_gpus, 1, 2, 1}, {CUOPT_NUM_GPUS, &mip_settings.num_gpus, 1, 2, 1}, - {CUOPT_MIP_BATCH_PDLP_STRONG_BRANCHING, &mip_settings.mip_batch_pdlp_strong_branching, 0, 1, 0}, + {CUOPT_MIP_BATCH_PDLP_STRONG_BRANCHING, &mip_settings.mip_batch_pdlp_strong_branching, 0, 2, 0}, {CUOPT_PRESOLVE, reinterpret_cast(&pdlp_settings.presolver), CUOPT_PRESOLVE_DEFAULT, CUOPT_PRESOLVE_PSLP, CUOPT_PRESOLVE_DEFAULT}, {CUOPT_PRESOLVE, reinterpret_cast(&mip_settings.presolver), CUOPT_PRESOLVE_DEFAULT, CUOPT_PRESOLVE_PSLP, CUOPT_PRESOLVE_DEFAULT}, {CUOPT_MIP_DETERMINISM_MODE, &mip_settings.determinism_mode, CUOPT_MODE_OPPORTUNISTIC, CUOPT_MODE_DETERMINISTIC, CUOPT_MODE_OPPORTUNISTIC}, From 3ae3f22461c11323c8c786dd3068289e5e19bcfc Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Tue, 17 Mar 2026 18:12:45 +0100 Subject: [PATCH 03/19] reliability branching now rank the candidates based on the objective change estimate via dual simplex single pivot (#963). --- cpp/src/branch_and_bound/branch_and_bound.cpp | 7 +- cpp/src/branch_and_bound/pseudo_costs.cpp | 214 +++++++++++------- cpp/src/branch_and_bound/pseudo_costs.hpp | 17 +- 3 files changed, 153 insertions(+), 85 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 1f24ea225a..b549932ead 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -839,11 +839,11 @@ branch_variable_t branch_and_bound_t::variable_selection( if (settings_.reliability_branching != 0) { branch_var = pc_.reliable_variable_selection(node_ptr, fractional, - solution, - settings_, - var_types_, worker, + var_types_, exploration_stats_, + settings_, + node_ptr->lower_bound, upper_bound_, worker_pool_.num_idle_workers(), log); @@ -2445,6 +2445,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); pc_.resize(original_lp_.num_cols); + original_lp_.A.transpose(pc_.AT); { raft::common::nvtx::range scope_sb("BB::strong_branching"); strong_branching(original_problem_, diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index edeb903da8..499d85437b 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -22,7 +22,6 @@ namespace cuopt::linear_programming::dual_simplex { namespace { - template struct objective_estimate_t { f_t down_obj; @@ -37,26 +36,21 @@ f_t compute_step_length(const simplex_solver_settings_t& settings, const std::vector& delta_z_indices) { f_t step_length = inf; - f_t pivot_tol = settings.pivot_tol; - const i_t nz = delta_z_indices.size(); + f_t pivot_tol = settings.pivot_tol; + const i_t nz = delta_z_indices.size(); for (i_t h = 0; h < nz; h++) { const i_t j = delta_z_indices[h]; if (vstatus[j] == variable_status_t::NONBASIC_LOWER && delta_z[j] < -pivot_tol) { const f_t ratio = -z[j] / delta_z[j]; - if (ratio < step_length) { - step_length = ratio; - } + if (ratio < step_length) { step_length = ratio; } } else if (vstatus[j] == variable_status_t::NONBASIC_UPPER && delta_z[j] > pivot_tol) { const f_t ratio = -z[j] / delta_z[j]; - if (ratio < step_length) { - step_length = ratio; - } + if (ratio < step_length) { step_length = ratio; } } } return step_length; } - template objective_estimate_t single_pivot_objective_estimate( const lp_problem_t& lp, @@ -65,7 +59,7 @@ objective_estimate_t single_pivot_objective_estimate( const std::vector& vstatus, i_t variable_j, i_t basic_j, - f_t root_obj, + f_t leaf_obj, const std::vector& x, const std::vector& y, const std::vector& z, @@ -96,7 +90,7 @@ objective_estimate_t single_pivot_objective_estimate( } work_estimate += nz_delta_y; const f_t delta_y_nz_percentage = delta_y_nz0 / static_cast(lp.num_rows) * 100.0; - const bool use_transpose = delta_y_nz_percentage <= 30.0; + const bool use_transpose = delta_y_nz_percentage <= 30.0; std::vector delta_z_indices; // delta_z starts out all zero if (use_transpose) { @@ -124,8 +118,6 @@ objective_estimate_t single_pivot_objective_estimate( work_estimate); } - - // Verify dual feasibility #ifdef CHECK_DUAL_FEASIBILITY { @@ -161,55 +153,62 @@ objective_estimate_t single_pivot_objective_estimate( delta_obj_check += std::floor(x[variable_j]) * delta_z[variable_j]; delta_obj_check *= step_length; if (std::abs(delta_obj_check - delta_obj) > 1e-6) { - settings.log.printf("Delta obj check %e. Delta obj %e. Step length %e.\n", delta_obj_check, delta_obj, step_length); + settings.log.printf("Delta obj check %e. Delta obj %e. Step length %e.\n", + delta_obj_check, + delta_obj, + step_length); } #endif - f_t down_obj = root_obj + delta_obj; - settings.log.printf("Down branch %d. Step length: %e. Objective estimate: %e. Delta obj: %e. Root obj: %e.\n", variable_j, step_length, down_obj, delta_obj, root_obj); + f_t down_obj = leaf_obj + std::max(delta_obj, 0); + settings.log.debug( + "Down branch %d. Step length: %e. Objective estimate: %e. Delta obj: %e. Root obj: %e.\n", + variable_j, + step_length, + down_obj, + delta_obj, + leaf_obj); // Up branch direction = 1; // Negate delta_z - for (i_t j : delta_z_indices) - { + for (i_t j : delta_z_indices) { delta_z[j] *= -1.0; } // Compute the step-length step_length = compute_step_length(settings, vstatus, z, delta_z, delta_z_indices); - delta_obj = step_length * (std::ceil(x[variable_j]) - x[variable_j]); - f_t up_obj = root_obj + delta_obj; - settings.log.printf("Up branch %d. Step length: %e. Objective estimate: %e. Delta obj: %e. Root obj: %e.\n", variable_j, step_length, up_obj, delta_obj, root_obj); + delta_obj = step_length * (std::ceil(x[variable_j]) - x[variable_j]); + f_t up_obj = leaf_obj + std::max(delta_obj, 0); + settings.log.debug( + "Up branch %d. Step length: %e. Objective estimate: %e. Delta obj: %e. Root obj: %e.\n", + variable_j, + step_length, + up_obj, + delta_obj, + leaf_obj); delta_z_indices.push_back(variable_j); // Clear delta_z - for (i_t j : delta_z_indices) - { - delta_z[j] = 0.0; + for (i_t j : delta_z_indices) { + delta_z[j] = 0.0; workspace[j] = 0; } #ifdef CHECK_DELTA_Z for (i_t j = 0; j < lp.num_cols; j++) { - if (delta_z[j] != 0.0) { - settings.log.printf("Delta z %d: %e\n", j, delta_z[j]); - } + if (delta_z[j] != 0.0) { settings.log.printf("Delta z %d: %e\n", j, delta_z[j]); } } for (i_t j = 0; j < lp.num_cols; j++) { - if (workspace[j] != 0) { - settings.log.printf("Workspace %d: %d\n", j, workspace[j]); - } + if (workspace[j] != 0) { settings.log.printf("Workspace %d: %d\n", j, workspace[j]); } } #endif - return objective_estimate_t{.down_obj = down_obj, .up_obj = up_obj}; } - template void initialize_pseudo_costs_with_estimate(const lp_problem_t& lp, const simplex_solver_settings_t& settings, @@ -227,9 +226,6 @@ void initialize_pseudo_costs_with_estimate(const lp_problem_t& lp, i_t m = lp.num_rows; i_t n = lp.num_cols; - csc_matrix_t A_transpose(m, n, 0); - lp.A.transpose(A_transpose); - std::vector delta_z(n, 0); std::vector workspace(n, 0); @@ -246,32 +242,29 @@ void initialize_pseudo_costs_with_estimate(const lp_problem_t& lp, } for (i_t k = 0; k < fractional.size(); k++) { - const i_t j = fractional[k]; + const i_t j = fractional[k]; objective_estimate_t estimate = single_pivot_objective_estimate(lp, - settings, - A_transpose, - vstatus, - j, - basic_map[j], - root_obj, - x, - y, - z, - basic_list, - nonbasic_list, - nonbasic_mark, - basis_factors, - workspace, - delta_z, - work_estimate); - pc.strong_branch_down[k] = estimate.down_obj; - pc.strong_branch_up[k] = estimate.up_obj; + settings, + pc.AT, + vstatus, + j, + basic_map[j], + root_obj, + x, + y, + z, + basic_list, + nonbasic_list, + nonbasic_mark, + basis_factors, + workspace, + delta_z, + work_estimate); + pc.strong_branch_down[k] = estimate.down_obj; + pc.strong_branch_up[k] = estimate.up_obj; } } - - - template void strong_branch_helper(i_t start, i_t end, @@ -496,15 +489,15 @@ f_t trial_branching(const lp_problem_t& original_lp, } // namespace template -static cuopt::mps_parser::mps_data_model_t simplex_problem_to_mps_data_model( - const dual_simplex::user_problem_t& user_problem) +static mps_parser::mps_data_model_t simplex_problem_to_mps_data_model( + const user_problem_t& user_problem) { - cuopt::mps_parser::mps_data_model_t mps_model; + mps_parser::mps_data_model_t mps_model; int m = user_problem.num_rows; int n = user_problem.num_cols; // Convert CSC to CSR using built-in method - dual_simplex::csr_matrix_t csr_A(m, n, 0); + csr_matrix_t csr_A(m, n, 0); user_problem.A.to_compressed_row(csr_A); int nz = csr_A.row_start[m]; @@ -827,13 +820,13 @@ i_t pseudo_costs_t::variable_selection(const std::vector& fractio template i_t pseudo_costs_t::reliable_variable_selection( - mip_node_t* node_ptr, + const mip_node_t* node_ptr, const std::vector& fractional, - const std::vector& solution, - const simplex_solver_settings_t& settings, - const std::vector& var_types, branch_and_bound_worker_t* worker, + const std::vector& var_types, const branch_and_bound_stats_t& bnb_stats, + const simplex_solver_settings_t& settings, + f_t leaf_obj, f_t upper_bound, int max_num_tasks, logger_t& log) @@ -847,6 +840,8 @@ i_t pseudo_costs_t::reliable_variable_selection( f_t pseudo_cost_down_avg; f_t pseudo_cost_up_avg; + lp_solution_t& leaf_solution = worker->leaf_solution; + initialized(num_initialized_down, num_initialized_up, pseudo_cost_down_avg, pseudo_cost_up_avg); log.printf("PC: num initialized down %d up %d avg down %e up %e\n", @@ -888,7 +883,8 @@ i_t pseudo_costs_t::reliable_variable_selection( continue; } - f_t score = calculate_pseudocost_score(j, solution, pseudo_cost_up_avg, pseudo_cost_down_avg); + f_t score = + calculate_pseudocost_score(j, leaf_solution.x, pseudo_cost_up_avg, pseudo_cost_down_avg); if (score > max_score) { max_score = score; @@ -897,8 +893,10 @@ i_t pseudo_costs_t::reliable_variable_selection( } if (unreliable_list.empty()) { - log.printf( - "pc branching on %d. Value %e. Score %e\n", branch_var, solution[branch_var], max_score); + log.printf("pc branching on %d. Value %e. Score %e\n", + branch_var, + leaf_solution.x[branch_var], + max_score); return branch_var; } @@ -921,8 +919,69 @@ i_t pseudo_costs_t::reliable_variable_selection( num_tasks, reliable_threshold); - // Shuffle the unreliable list so every variable has the same chance to be selected. - if (unreliable_list.size() > max_num_candidates) { worker->rng.shuffle(unreliable_list); } + if (unreliable_list.size() > max_num_candidates) { + if (reliability_branching_settings.rank_candidates_with_dual_pivot) { + i_t m = worker->leaf_problem.num_rows; + i_t n = worker->leaf_problem.num_cols; + f_t work_estimate = 0; + + std::vector delta_z(n, 0); + std::vector workspace(n, 0); + std::vector> unreliable_list_with_obj; // {index, objective change} + unreliable_list_with_obj.reserve(unreliable_list.size()); + + std::vector basic_map(n, -1); + for (i_t i = 0; i < m; i++) { + basic_map[worker->basic_list[i]] = i; + } + + std::vector nonbasic_mark(n, -1); + for (i_t i = 0; i < n - m; i++) { + nonbasic_mark[worker->nonbasic_list[i]] = i; + } + + for (i_t i = 0; i < unreliable_list.size(); ++i) { + i_t j = unreliable_list[i]; + + auto [up_obj, down_obj] = single_pivot_objective_estimate(worker->leaf_problem, + settings, + AT, + node_ptr->vstatus, + j, + basic_map[j], + leaf_obj, + leaf_solution.x, + leaf_solution.y, + leaf_solution.z, + worker->basic_list, + worker->nonbasic_list, + nonbasic_mark, + worker->basis_factors, + workspace, + delta_z, + work_estimate); + + unreliable_list_with_obj.push_back({j, std::min(up_obj, down_obj)}); + } + + // We only need to get the top-k elements in the list, where + // k = num_candidates + std::partial_sort(unreliable_list_with_obj.begin(), + unreliable_list_with_obj.begin() + num_candidates, + unreliable_list_with_obj.end(), + [](auto el1, auto el2) { return el1.second < el2.second; }); + + unreliable_list.clear(); + for (i_t i = 0; i < num_candidates; ++i) { + auto [j, obj] = unreliable_list_with_obj[i]; + unreliable_list.push_back(j); + } + + } else { + // Shuffle the unreliable list so every variable has the same chance to be selected. + worker->rng.shuffle(unreliable_list); + } + } if (toc(start_time) > settings.time_limit) { log.printf("Time limit reached"); @@ -949,7 +1008,7 @@ i_t pseudo_costs_t::reliable_variable_selection( worker->nonbasic_list, j, worker->leaf_problem.lower[j], - std::floor(solution[j]), + std::floor(leaf_solution.x[j]), upper_bound, branch_and_bound_lp_iter_per_node, start_time, @@ -959,7 +1018,7 @@ i_t pseudo_costs_t::reliable_variable_selection( if (!std::isnan(obj)) { f_t change_in_obj = std::max(obj - node_ptr->lower_bound, eps); - f_t change_in_x = solution[j] - std::floor(solution[j]); + f_t change_in_x = leaf_solution.x[j] - std::floor(leaf_solution.x[j]); pseudo_cost_sum_down[j] += change_in_obj / change_in_x; pseudo_cost_num_down[j]++; } @@ -979,7 +1038,7 @@ i_t pseudo_costs_t::reliable_variable_selection( worker->basic_list, worker->nonbasic_list, j, - std::ceil(solution[j]), + std::ceil(leaf_solution.x[j]), worker->leaf_problem.upper[j], upper_bound, branch_and_bound_lp_iter_per_node, @@ -990,7 +1049,7 @@ i_t pseudo_costs_t::reliable_variable_selection( if (!std::isnan(obj)) { f_t change_in_obj = std::max(obj - node_ptr->lower_bound, eps); - f_t change_in_x = std::ceil(solution[j]) - solution[j]; + f_t change_in_x = std::ceil(leaf_solution.x[j]) - leaf_solution.x[j]; pseudo_cost_sum_up[j] += change_in_obj / change_in_x; pseudo_cost_num_up[j]++; } @@ -999,7 +1058,8 @@ i_t pseudo_costs_t::reliable_variable_selection( if (toc(start_time) > settings.time_limit) { continue; } - f_t score = calculate_pseudocost_score(j, solution, pseudo_cost_up_avg, pseudo_cost_down_avg); + f_t score = + calculate_pseudocost_score(j, leaf_solution.x, pseudo_cost_up_avg, pseudo_cost_down_avg); score_mutex.lock(); if (score > max_score) { @@ -1010,7 +1070,7 @@ i_t pseudo_costs_t::reliable_variable_selection( } log.printf( - "pc branching on %d. Value %e. Score %e\n", branch_var, solution[branch_var], max_score); + "pc branching on %d. Value %e. Score %e\n", branch_var, leaf_solution.x[branch_var], max_score); return branch_var; } diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index 37b7a98060..c72f060a76 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.hpp +++ b/cpp/src/branch_and_bound/pseudo_costs.hpp @@ -403,6 +403,11 @@ struct reliability_branching_settings_t { // Only used when `reliable_threshold` is negative i_t max_reliable_threshold = 5; i_t min_reliable_threshold = 1; + + // Estimate the objective change of each fractional variable + // using a single pivot of dual simplex. Then rank the candidates + // based on this estimation. + bool rank_candidates_with_dual_pivot = true; }; template @@ -414,7 +419,8 @@ class pseudo_costs_t { pseudo_cost_num_down(num_variables), pseudo_cost_num_up(num_variables), pseudo_cost_mutex_up(num_variables), - pseudo_cost_mutex_down(num_variables) + pseudo_cost_mutex_down(num_variables), + AT(1, 1, 1) { } @@ -472,13 +478,13 @@ class pseudo_costs_t { const std::vector& solution, logger_t& log); - i_t reliable_variable_selection(mip_node_t* node_ptr, + i_t reliable_variable_selection(const mip_node_t* node_ptr, const std::vector& fractional, - const std::vector& solution, - const simplex_solver_settings_t& settings, - const std::vector& var_types, branch_and_bound_worker_t* worker, + const std::vector& var_types, const branch_and_bound_stats_t& bnb_stats, + const simplex_solver_settings_t& settings, + f_t leaf_obj, f_t upper_bound, int max_num_tasks, logger_t& log); @@ -504,6 +510,7 @@ class pseudo_costs_t { reliability_branching_settings_t reliability_branching_settings; + csc_matrix_t AT; // Transpose of the constraint matrix A std::vector> pseudo_cost_sum_up; std::vector> pseudo_cost_sum_down; std::vector> pseudo_cost_num_up; From 2b2f94dec6fbe328eb388decddd6d3540b29e369 Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Wed, 18 Mar 2026 18:13:58 +0100 Subject: [PATCH 04/19] `single_pivot_objective_change_estimate` now returns the objective change instead of the objective. fixed candidate ranking in reliability branching. --- cpp/src/branch_and_bound/branch_and_bound.cpp | 12 +- cpp/src/branch_and_bound/pseudo_costs.cpp | 218 ++++++++---------- cpp/src/branch_and_bound/pseudo_costs.hpp | 6 +- cpp/src/math_optimization/solver_settings.cu | 2 +- 4 files changed, 108 insertions(+), 130 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index b549932ead..2b5450cb4f 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1971,6 +1971,14 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.log.printf("Reduced cost strengthening enabled: %d\n", settings_.reduced_cost_strengthening); + if (settings_.reliability_branching == -2) { + settings_.log.printf("Enabled objective estimation via dual simplex\n"); + pc_.reliability_branching_settings.rank_candidates_with_dual_pivot = true; + } else { + settings_.log.printf("Disabled objective estimation via dual simplex\n"); + pc_.reliability_branching_settings.rank_candidates_with_dual_pivot = false; + } + variable_bounds_t variable_bounds( original_lp_, settings_, var_types_, Arow_, new_slacks_); @@ -2453,9 +2461,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_, exploration_stats_.start_time, var_types_, - root_relax_soln_.x, - root_relax_soln_.y, - root_relax_soln_.z, + root_relax_soln_, fractional, root_objective_, upper_bound_, diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index 499d85437b..ea9d4009f9 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -23,9 +23,9 @@ namespace cuopt::linear_programming::dual_simplex { namespace { template -struct objective_estimate_t { - f_t down_obj; - f_t up_obj; +struct objective_change_estimate_t { + f_t down_obj_change; + f_t up_obj_change; }; template @@ -52,17 +52,14 @@ f_t compute_step_length(const simplex_solver_settings_t& settings, } template -objective_estimate_t single_pivot_objective_estimate( +objective_change_estimate_t single_pivot_objective_change_estimate( const lp_problem_t& lp, const simplex_solver_settings_t& settings, const csc_matrix_t& A_transpose, const std::vector& vstatus, i_t variable_j, i_t basic_j, - f_t leaf_obj, - const std::vector& x, - const std::vector& y, - const std::vector& z, + const lp_solution_t& lp_solution, const std::vector& basic_list, const std::vector& nonbasic_list, const std::vector& nonbasic_mark, @@ -86,7 +83,7 @@ objective_estimate_t single_pivot_objective_estimate( i_t delta_y_nz0 = 0; const i_t nz_delta_y = delta_y.i.size(); for (i_t k = 0; k < nz_delta_y; k++) { - if (std::abs(delta_y.x[k]) > 1e-12) { delta_y_nz0++; } + if (std::abs(delta_y.x[k]) > settings.zero_tol) { delta_y_nz0++; } } work_estimate += nz_delta_y; const f_t delta_y_nz_percentage = delta_y_nz0 / static_cast(lp.num_rows) * 100.0; @@ -121,7 +118,7 @@ objective_estimate_t single_pivot_objective_estimate( // Verify dual feasibility #ifdef CHECK_DUAL_FEASIBILITY { - std::vector dual_residual = z; + std::vector dual_residual = lp_solution.z; for (i_t j = 0; j < lp.num_cols; j++) { dual_residual[j] -= lp.objective[j]; } @@ -132,15 +129,16 @@ objective_estimate_t single_pivot_objective_estimate( #endif // Compute the step-length - f_t step_length = compute_step_length(settings, vstatus, z, delta_z, delta_z_indices); + f_t step_length = compute_step_length(settings, vstatus, lp_solution.z, delta_z, delta_z_indices); // Handle the leaving variable case - f_t delta_obj = step_length * (x[variable_j] - std::floor(x[variable_j])); + f_t delta_obj_down = + step_length * (lp_solution.x[variable_j] - std::floor(lp_solution.x[variable_j])); #ifdef CHECK_DELTA_OBJ f_t delta_obj_check = 0.0; for (i_t k = 0; k < delta_y.i.size(); k++) { - delta_obj_check += lp.rhs[delta_y.i[k]] * delta_y.x[k]; + delta_obj_check += lp.rhs[delta_y.i[k]] * delta_y.lp_solution.x[k]; } for (i_t h = 0; h < delta_z_indices.size(); h++) { const i_t j = delta_z_indices[h]; @@ -150,7 +148,7 @@ objective_estimate_t single_pivot_objective_estimate( delta_obj_check += lp.upper[j] * delta_z[j]; } } - delta_obj_check += std::floor(x[variable_j]) * delta_z[variable_j]; + delta_obj_check += std::floor(lp_solution.x[variable_j]) * delta_z[variable_j]; delta_obj_check *= step_length; if (std::abs(delta_obj_check - delta_obj) > 1e-6) { settings.log.printf("Delta obj check %e. Delta obj %e. Step length %e.\n", @@ -160,14 +158,8 @@ objective_estimate_t single_pivot_objective_estimate( } #endif - f_t down_obj = leaf_obj + std::max(delta_obj, 0); settings.log.debug( - "Down branch %d. Step length: %e. Objective estimate: %e. Delta obj: %e. Root obj: %e.\n", - variable_j, - step_length, - down_obj, - delta_obj, - leaf_obj); + "Down branch %d. Step length: %e. Delta obj: %e. \n", variable_j, step_length, delta_obj_down); // Up branch direction = 1; @@ -177,17 +169,12 @@ objective_estimate_t single_pivot_objective_estimate( } // Compute the step-length - step_length = compute_step_length(settings, vstatus, z, delta_z, delta_z_indices); + step_length = compute_step_length(settings, vstatus, lp_solution.z, delta_z, delta_z_indices); - delta_obj = step_length * (std::ceil(x[variable_j]) - x[variable_j]); - f_t up_obj = leaf_obj + std::max(delta_obj, 0); + f_t delta_obj_up = + step_length * (std::ceil(lp_solution.x[variable_j]) - lp_solution.x[variable_j]); settings.log.debug( - "Up branch %d. Step length: %e. Objective estimate: %e. Delta obj: %e. Root obj: %e.\n", - variable_j, - step_length, - up_obj, - delta_obj, - leaf_obj); + "Up branch %d. Step length: %e. Delta obj: %e.\n", variable_j, step_length, delta_obj_up); delta_z_indices.push_back(variable_j); @@ -206,20 +193,18 @@ objective_estimate_t single_pivot_objective_estimate( } #endif - return objective_estimate_t{.down_obj = down_obj, .up_obj = up_obj}; + return {.down_obj_change = std::max(delta_obj_down, 0), + .up_obj_change = std::max(delta_obj_up, 0)}; } template void initialize_pseudo_costs_with_estimate(const lp_problem_t& lp, const simplex_solver_settings_t& settings, const std::vector& vstatus, - const std::vector& x, - const std::vector& y, - const std::vector& z, + const lp_solution_t& lp_solution, const std::vector& basic_list, const std::vector& nonbasic_list, const std::vector& fractional, - f_t root_obj, basis_update_mpf_t& basis_factors, pseudo_costs_t& pc) { @@ -242,26 +227,24 @@ void initialize_pseudo_costs_with_estimate(const lp_problem_t& lp, } for (i_t k = 0; k < fractional.size(); k++) { - const i_t j = fractional[k]; - objective_estimate_t estimate = single_pivot_objective_estimate(lp, - settings, - pc.AT, - vstatus, - j, - basic_map[j], - root_obj, - x, - y, - z, - basic_list, - nonbasic_list, - nonbasic_mark, - basis_factors, - workspace, - delta_z, - work_estimate); - pc.strong_branch_down[k] = estimate.down_obj; - pc.strong_branch_up[k] = estimate.up_obj; + const i_t j = fractional[k]; + objective_change_estimate_t estimate = + single_pivot_objective_change_estimate(lp, + settings, + pc.AT, + vstatus, + j, + basic_map[j], + lp_solution, + basic_list, + nonbasic_list, + nonbasic_mark, + basis_factors, + workspace, + delta_z, + work_estimate); + pc.strong_branch_down[k] = estimate.down_obj_change; + pc.strong_branch_up[k] = estimate.up_obj_change; } } @@ -572,9 +555,7 @@ void strong_branching(const user_problem_t& original_problem, const simplex_solver_settings_t& settings, f_t start_time, const std::vector& var_types, - const std::vector& root_x, - const std::vector& root_y, - const std::vector& root_z, + const lp_solution_t& root_solution, const std::vector& fractional, f_t root_obj, f_t upper_bound, @@ -604,7 +585,7 @@ void strong_branching(const user_problem_t& original_problem, // Convert the root_soln to the original problem space std::vector original_root_soln_x; - uncrush_primal_solution(original_problem, original_lp, root_x, original_root_soln_x); + uncrush_primal_solution(original_problem, original_lp, root_solution.x, original_root_soln_x); std::vector fraction_values; @@ -674,13 +655,10 @@ void strong_branching(const user_problem_t& original_problem, initialize_pseudo_costs_with_estimate(original_lp, settings, root_vstatus, - root_x, - root_y, - root_z, + root_solution, basic_list, nonbasic_list, fractional, - root_obj, basis_factors, pc); @@ -720,7 +698,7 @@ void strong_branching(const user_problem_t& original_problem, fractional, root_obj, upper_bound, - root_x, + root_solution.x, root_vstatus, edge_norms, pc); @@ -729,7 +707,7 @@ void strong_branching(const user_problem_t& original_problem, settings.log.printf("Strong branching completed in %.2fs\n", toc(strong_branching_start_time)); } - pc.update_pseudo_costs_from_strong_branching(fractional, root_x); + pc.update_pseudo_costs_from_strong_branching(fractional, root_solution.x); } template @@ -851,7 +829,6 @@ i_t pseudo_costs_t::reliable_variable_selection( pseudo_cost_up_avg); const int64_t branch_and_bound_lp_iters = bnb_stats.total_lp_iters; - const int64_t branch_and_bound_explored = bnb_stats.nodes_explored; const i_t branch_and_bound_lp_iter_per_node = branch_and_bound_lp_iters / bnb_stats.nodes_explored; @@ -873,22 +850,22 @@ i_t pseudo_costs_t::reliable_variable_selection( reliable_threshold = strong_branching_lp_iter < max_reliability_iter ? reliable_threshold : 0; } - std::vector unreliable_list; + std::vector> unreliable_list; omp_mutex_t score_mutex; for (i_t j : fractional) { - if (pseudo_cost_num_down[j] < reliable_threshold || - pseudo_cost_num_up[j] < reliable_threshold) { - unreliable_list.push_back(j); - continue; - } - f_t score = calculate_pseudocost_score(j, leaf_solution.x, pseudo_cost_up_avg, pseudo_cost_down_avg); - if (score > max_score) { - max_score = score; - branch_var = j; + if (pseudo_cost_num_down[j] < reliable_threshold || + pseudo_cost_num_up[j] < reliable_threshold) { + if (pseudo_cost_num_down[j] == 0 || pseudo_cost_num_up[j] == 0) { score = -1; } + unreliable_list.push_back(std::make_pair(score, j)); + } else { + if (score > max_score) { + max_score = score; + branch_var = j; + } } } @@ -901,7 +878,7 @@ i_t pseudo_costs_t::reliable_variable_selection( return branch_var; } - const int num_tasks = std::max(max_num_tasks, 1); + const int num_tasks = std::max(max_num_tasks, 10); const int task_priority = reliability_branching_settings.task_priority; const i_t max_num_candidates = reliability_branching_settings.max_num_candidates; const i_t num_candidates = std::min(unreliable_list.size(), max_num_candidates); @@ -920,6 +897,15 @@ i_t pseudo_costs_t::reliable_variable_selection( reliable_threshold); if (unreliable_list.size() > max_num_candidates) { + std::cout << std::format( + "RB iters = {}, B&B iters = {}, unreliable = {}, num_tasks = {}, " + "reliable_threshold = {}\n", + strong_branching_lp_iter.load(), + branch_and_bound_lp_iters, + unreliable_list.size(), + num_tasks, + reliable_threshold); + if (reliability_branching_settings.rank_candidates_with_dual_pivot) { i_t m = worker->leaf_problem.num_rows; i_t n = worker->leaf_problem.num_cols; @@ -927,8 +913,6 @@ i_t pseudo_costs_t::reliable_variable_selection( std::vector delta_z(n, 0); std::vector workspace(n, 0); - std::vector> unreliable_list_with_obj; // {index, objective change} - unreliable_list_with_obj.reserve(unreliable_list.size()); std::vector basic_map(n, -1); for (i_t i = 0; i < m; i++) { @@ -940,47 +924,39 @@ i_t pseudo_costs_t::reliable_variable_selection( nonbasic_mark[worker->nonbasic_list[i]] = i; } - for (i_t i = 0; i < unreliable_list.size(); ++i) { - i_t j = unreliable_list[i]; - - auto [up_obj, down_obj] = single_pivot_objective_estimate(worker->leaf_problem, - settings, - AT, - node_ptr->vstatus, - j, - basic_map[j], - leaf_obj, - leaf_solution.x, - leaf_solution.y, - leaf_solution.z, - worker->basic_list, - worker->nonbasic_list, - nonbasic_mark, - worker->basis_factors, - workspace, - delta_z, - work_estimate); - - unreliable_list_with_obj.push_back({j, std::min(up_obj, down_obj)}); - } - - // We only need to get the top-k elements in the list, where - // k = num_candidates - std::partial_sort(unreliable_list_with_obj.begin(), - unreliable_list_with_obj.begin() + num_candidates, - unreliable_list_with_obj.end(), - [](auto el1, auto el2) { return el1.second < el2.second; }); - - unreliable_list.clear(); - for (i_t i = 0; i < num_candidates; ++i) { - auto [j, obj] = unreliable_list_with_obj[i]; - unreliable_list.push_back(j); + for (auto& [score, j] : unreliable_list) { + if (score == -1) { + objective_change_estimate_t estimate = + single_pivot_objective_change_estimate(worker->leaf_problem, + settings, + AT, + node_ptr->vstatus, + j, + basic_map[j], + leaf_solution, + worker->basic_list, + worker->nonbasic_list, + nonbasic_mark, + worker->basis_factors, + workspace, + delta_z, + work_estimate); + + score = std::max(estimate.up_obj_change, eps) * std::max(estimate.down_obj_change, eps); + } } - } else { - // Shuffle the unreliable list so every variable has the same chance to be selected. - worker->rng.shuffle(unreliable_list); + for (auto& [score, j] : unreliable_list) { + if (score == -1) { score = worker->rng.uniform(eps, max_score); } + } } + + // We only need to get the top-k elements in the list, where + // k = num_candidates + std::partial_sort(unreliable_list.begin(), + unreliable_list.begin() + num_candidates, + unreliable_list.end(), + [](auto el1, auto el2) { return el1.first > el2.first; }); } if (toc(start_time) > settings.time_limit) { @@ -991,7 +967,7 @@ i_t pseudo_costs_t::reliable_variable_selection( #pragma omp taskloop if (num_tasks > 1) priority(task_priority) num_tasks(num_tasks) \ shared(score_mutex) for (i_t i = 0; i < num_candidates; ++i) { - const i_t j = unreliable_list[i]; + auto [score, j] = unreliable_list[i]; if (toc(start_time) > settings.time_limit) { continue; } @@ -1058,7 +1034,7 @@ i_t pseudo_costs_t::reliable_variable_selection( if (toc(start_time) > settings.time_limit) { continue; } - f_t score = + score = calculate_pseudocost_score(j, leaf_solution.x, pseudo_cost_up_avg, pseudo_cost_down_avg); score_mutex.lock(); @@ -1139,9 +1115,7 @@ template void strong_branching(const user_problem_t& o const simplex_solver_settings_t& settings, double start_time, const std::vector& var_types, - const std::vector& root_x, - const std::vector& root_y, - const std::vector& root_z, + const lp_solution_t& root_solution, const std::vector& fractional, double root_obj, double upper_bound, diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index c72f060a76..f7970b30d1 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.hpp +++ b/cpp/src/branch_and_bound/pseudo_costs.hpp @@ -407,7 +407,7 @@ struct reliability_branching_settings_t { // Estimate the objective change of each fractional variable // using a single pivot of dual simplex. Then rank the candidates // based on this estimation. - bool rank_candidates_with_dual_pivot = true; + bool rank_candidates_with_dual_pivot = false; }; template @@ -529,9 +529,7 @@ void strong_branching(const user_problem_t& original_problem, const simplex_solver_settings_t& settings, f_t start_time, const std::vector& var_types, - const std::vector& root_x, - const std::vector& root_y, - const std::vector& root_z, + const lp_solution_t& root_solution, const std::vector& fractional, f_t root_obj, f_t upper_bound, diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 13a8f74145..3bb10096d5 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -104,7 +104,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_PRESOLVE, reinterpret_cast(&mip_settings.presolver), CUOPT_PRESOLVE_DEFAULT, CUOPT_PRESOLVE_PSLP, CUOPT_PRESOLVE_DEFAULT}, {CUOPT_MIP_DETERMINISM_MODE, &mip_settings.determinism_mode, CUOPT_MODE_OPPORTUNISTIC, CUOPT_MODE_DETERMINISTIC, CUOPT_MODE_OPPORTUNISTIC}, {CUOPT_RANDOM_SEED, &mip_settings.seed, -1, std::numeric_limits::max(), -1}, - {CUOPT_MIP_RELIABILITY_BRANCHING, &mip_settings.reliability_branching, -1, std::numeric_limits::max(), -1}, + {CUOPT_MIP_RELIABILITY_BRANCHING, &mip_settings.reliability_branching, -2, std::numeric_limits::max(), -1}, {CUOPT_PDLP_PRECISION, reinterpret_cast(&pdlp_settings.pdlp_precision), CUOPT_PDLP_DEFAULT_PRECISION, CUOPT_PDLP_MIXED_PRECISION, CUOPT_PDLP_DEFAULT_PRECISION} }; From 925d36e157fc760f050cf9a2a5e89c489fcba5f5 Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Thu, 19 Mar 2026 11:34:29 +0100 Subject: [PATCH 05/19] removed new cutoff in trial branching. adjusted reliability branching task settings. --- cpp/src/branch_and_bound/branch_and_bound.cpp | 1 - cpp/src/branch_and_bound/pseudo_costs.cpp | 51 +++++++++---------- cpp/src/branch_and_bound/pseudo_costs.hpp | 1 - 3 files changed, 25 insertions(+), 28 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 2b5450cb4f..1ed6de776d 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -843,7 +843,6 @@ branch_variable_t branch_and_bound_t::variable_selection( var_types_, exploration_stats_, settings_, - node_ptr->lower_bound, upper_bound_, worker_pool_.num_idle_workers(), log); diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index ea9d4009f9..03758f9680 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -290,18 +290,18 @@ void strong_branch_helper(i_t start, child_settings.time_limit = std::max(0.0, settings.time_limit - elapsed_time); child_settings.iteration_limit = 200; - if (std::isfinite(upper_bound)) { - child_settings.cut_off = upper_bound + settings.dual_tol; - } else { - child_settings.cut_off = 0; - for (i_t i = 0; i < original_lp.num_cols; ++i) { - if (original_lp.objective[i] < 0) { - child_settings.cut_off += original_lp.objective[i] * child_problem.upper[i]; - } else if (original_lp.objective[i] > 0) { - child_settings.cut_off += original_lp.objective[i] * child_problem.lower[i]; - } - } - } + // if (std::isfinite(upper_bound)) { + child_settings.cut_off = upper_bound + settings.dual_tol; + // } else { + // child_settings.cut_off = 0; + // for (i_t i = 0; i < original_lp.num_cols; ++i) { + // if (original_lp.objective[i] < 0) { + // child_settings.cut_off += original_lp.objective[i] * child_problem.upper[i]; + // } else if (original_lp.objective[i] > 0) { + // child_settings.cut_off += original_lp.objective[i] * child_problem.lower[i]; + // } + // } + // } lp_solution_t solution(original_lp.num_rows, original_lp.num_cols); i_t iter = 0; @@ -412,18 +412,18 @@ f_t trial_branching(const lp_problem_t& original_lp, child_settings.inside_mip = 2; child_settings.scale_columns = false; - if (std::isfinite(upper_bound)) { - child_settings.cut_off = upper_bound + settings.dual_tol; - } else { - child_settings.cut_off = 0; - for (i_t i = 0; i < original_lp.num_cols; ++i) { - if (original_lp.objective[i] < 0) { - child_settings.cut_off += original_lp.objective[i] * child_problem.upper[i]; - } else if (original_lp.objective[i] > 0) { - child_settings.cut_off += original_lp.objective[i] * child_problem.lower[i]; - } - } - } + // if (std::isfinite(upper_bound)) { + child_settings.cut_off = upper_bound + settings.dual_tol; + // } else { + // child_settings.cut_off = 0; + // for (i_t i = 0; i < original_lp.num_cols; ++i) { + // if (original_lp.objective[i] < 0) { + // child_settings.cut_off += original_lp.objective[i] * child_problem.upper[i]; + // } else if (original_lp.objective[i] > 0) { + // child_settings.cut_off += original_lp.objective[i] * child_problem.lower[i]; + // } + // } + // } lp_solution_t solution(original_lp.num_rows, original_lp.num_cols); i_t iter = 0; @@ -804,7 +804,6 @@ i_t pseudo_costs_t::reliable_variable_selection( const std::vector& var_types, const branch_and_bound_stats_t& bnb_stats, const simplex_solver_settings_t& settings, - f_t leaf_obj, f_t upper_bound, int max_num_tasks, logger_t& log) @@ -964,7 +963,7 @@ i_t pseudo_costs_t::reliable_variable_selection( return branch_var; } -#pragma omp taskloop if (num_tasks > 1) priority(task_priority) num_tasks(num_tasks) \ +#pragma omp taskloop if (num_candidates > 5) priority(task_priority) num_tasks(num_tasks) \ shared(score_mutex) for (i_t i = 0; i < num_candidates; ++i) { auto [score, j] = unreliable_list[i]; diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index f7970b30d1..e21034f08c 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.hpp +++ b/cpp/src/branch_and_bound/pseudo_costs.hpp @@ -484,7 +484,6 @@ class pseudo_costs_t { const std::vector& var_types, const branch_and_bound_stats_t& bnb_stats, const simplex_solver_settings_t& settings, - f_t leaf_obj, f_t upper_bound, int max_num_tasks, logger_t& log); From c6b0e3c0556e1977844c8a71eb2bfb52fd2773ab Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Thu, 19 Mar 2026 15:35:55 +0100 Subject: [PATCH 06/19] small fixes --- cpp/src/branch_and_bound/pseudo_costs.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index 03758f9680..7ef276630d 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -94,7 +94,7 @@ objective_change_estimate_t single_pivot_objective_change_estimate( compute_delta_z(A_transpose, delta_y, variable_j, - -1, + direction, nonbasic_mark, workspace, delta_z_indices, @@ -108,7 +108,7 @@ objective_change_estimate_t single_pivot_objective_change_estimate( nonbasic_list, delta_y_dense, variable_j, - -1, + direction, workspace, delta_z_indices, delta_z, @@ -829,7 +829,7 @@ i_t pseudo_costs_t::reliable_variable_selection( const int64_t branch_and_bound_lp_iters = bnb_stats.total_lp_iters; const i_t branch_and_bound_lp_iter_per_node = - branch_and_bound_lp_iters / bnb_stats.nodes_explored; + bnb_stats.nodes_explored > 0 ? branch_and_bound_lp_iters / bnb_stats.nodes_explored : 0; i_t reliable_threshold = settings.reliability_branching; if (reliable_threshold < 0) { @@ -945,8 +945,11 @@ i_t pseudo_costs_t::reliable_variable_selection( } } } else { + f_t high = max_score > 0 ? max_score : 1; + f_t low = 0; + for (auto& [score, j] : unreliable_list) { - if (score == -1) { score = worker->rng.uniform(eps, max_score); } + if (score == -1) { score = worker->rng.uniform(low, high); } } } From 703f30c588d9be59ad259d088459ccb9df5eaad6 Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Sat, 21 Mar 2026 16:16:00 +0100 Subject: [PATCH 07/19] re-enable cutoff for trial branching --- cpp/src/branch_and_bound/pseudo_costs.cpp | 55 ++++++++++++----------- 1 file changed, 30 insertions(+), 25 deletions(-) diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index 7ef276630d..fc9cf79343 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -69,6 +69,8 @@ objective_change_estimate_t single_pivot_objective_change_estimate( f_t& work_estimate) { // Compute the objective estimate for the down and up branches of variable j + assert(variable_j >= 0); + assert(basic_j >= 0); // Down branch i_t direction = -1; @@ -228,6 +230,8 @@ void initialize_pseudo_costs_with_estimate(const lp_problem_t& lp, for (i_t k = 0; k < fractional.size(); k++) { const i_t j = fractional[k]; + assert(j >= 0); + objective_change_estimate_t estimate = single_pivot_objective_change_estimate(lp, settings, @@ -290,18 +294,18 @@ void strong_branch_helper(i_t start, child_settings.time_limit = std::max(0.0, settings.time_limit - elapsed_time); child_settings.iteration_limit = 200; - // if (std::isfinite(upper_bound)) { - child_settings.cut_off = upper_bound + settings.dual_tol; - // } else { - // child_settings.cut_off = 0; - // for (i_t i = 0; i < original_lp.num_cols; ++i) { - // if (original_lp.objective[i] < 0) { - // child_settings.cut_off += original_lp.objective[i] * child_problem.upper[i]; - // } else if (original_lp.objective[i] > 0) { - // child_settings.cut_off += original_lp.objective[i] * child_problem.lower[i]; - // } - // } - // } + if (std::isfinite(upper_bound)) { + child_settings.cut_off = upper_bound + settings.dual_tol; + } else { + child_settings.cut_off = 0; + for (i_t i = 0; i < original_lp.num_cols; ++i) { + if (original_lp.objective[i] < 0) { + child_settings.cut_off += original_lp.objective[i] * child_problem.upper[i]; + } else if (original_lp.objective[i] > 0) { + child_settings.cut_off += original_lp.objective[i] * child_problem.lower[i]; + } + } + } lp_solution_t solution(original_lp.num_rows, original_lp.num_cols); i_t iter = 0; @@ -321,7 +325,8 @@ void strong_branch_helper(i_t start, if (status == dual::status_t::DUAL_UNBOUNDED) { // LP was infeasible obj = std::numeric_limits::infinity(); - } else if (status == dual::status_t::OPTIMAL || status == dual::status_t::ITERATION_LIMIT) { + } else if (status == dual::status_t::OPTIMAL || status == dual::status_t::ITERATION_LIMIT || + status == dual::status_t::CUTOFF) { obj = compute_objective(child_problem, solution.x); } else { settings.log.debug("Thread id %2d remaining %d variable %d branch %d status %d\n", @@ -412,18 +417,18 @@ f_t trial_branching(const lp_problem_t& original_lp, child_settings.inside_mip = 2; child_settings.scale_columns = false; - // if (std::isfinite(upper_bound)) { - child_settings.cut_off = upper_bound + settings.dual_tol; - // } else { - // child_settings.cut_off = 0; - // for (i_t i = 0; i < original_lp.num_cols; ++i) { - // if (original_lp.objective[i] < 0) { - // child_settings.cut_off += original_lp.objective[i] * child_problem.upper[i]; - // } else if (original_lp.objective[i] > 0) { - // child_settings.cut_off += original_lp.objective[i] * child_problem.lower[i]; - // } - // } - // } + if (std::isfinite(upper_bound)) { + child_settings.cut_off = upper_bound + settings.dual_tol; + } else { + child_settings.cut_off = 0; + for (i_t i = 0; i < original_lp.num_cols; ++i) { + if (original_lp.objective[i] < 0) { + child_settings.cut_off += original_lp.objective[i] * child_problem.upper[i]; + } else if (original_lp.objective[i] > 0) { + child_settings.cut_off += original_lp.objective[i] * child_problem.lower[i]; + } + } + } lp_solution_t solution(original_lp.num_rows, original_lp.num_cols); i_t iter = 0; From bb71cde93ad34c34e9f4c6b4001adabd59e722df Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Mon, 23 Mar 2026 14:30:13 +0100 Subject: [PATCH 08/19] fixed style --- cpp/src/branch_and_bound/branch_and_bound.cpp | 19 ++++++++ cpp/src/dual_simplex/phase2.cpp | 45 +++++++++---------- 2 files changed, 39 insertions(+), 25 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 1ed6de776d..d6792320d0 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -789,6 +789,25 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, upper_bound_ = leaf_objective; report(feasible_solution_symbol(thread_type), leaf_objective, get_lower_bound(), leaf_depth, 0); send_solution = true; + + f_t upper = upper_bound_; + f_t lower = get_lower_bound(); + f_t user_upper = compute_user_objective(original_lp_, upper); + f_t user_lower = compute_user_objective(original_lp_, lower); + + std::cout << std::format( + "{}: user_obj={:.3g}, solver_obj={:.3g}, user_lower={:.3g}, " + "solver_lower={:.3g}, user_gap={:.3g}, " + "solver_gap={:.3g}, tol={:.3g}", + feasible_solution_symbol(thread_type), + user_upper, + upper, + user_lower, + lower, + std::abs(user_upper - user_lower), + std::abs(upper - lower), + settings_.absolute_mip_gap_tol) + << std::endl; } if (send_solution && settings_.solution_callback != nullptr) { diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index eca9dbd251..244060506f 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -82,7 +82,6 @@ class nvtx_range_guard { bool active_; }; - template void compute_reduced_cost_update(const lp_problem_t& lp, const std::vector& basic_list, @@ -446,8 +445,6 @@ void initial_perturbation(const lp_problem_t& lp, n); } - - template void compute_reduced_costs(const std::vector& objective, const csc_matrix_t& A, @@ -3638,28 +3635,26 @@ template dual::status_t dual_phase2_with_advanced_basis( std::vector& steepest_edge_norms, work_limit_context_t* work_unit_context); -template -void compute_reduced_cost_update(const lp_problem_t& lp, - const std::vector& basic_list, - const std::vector& nonbasic_list, - const std::vector& delta_y, - int leaving_index, - int direction, - std::vector& delta_z_mark, - std::vector& delta_z_indices, - std::vector& delta_z, - double& work_estimate); - -template -void compute_delta_z(const csc_matrix_t& A_transpose, - const sparse_vector_t& delta_y, - int leaving_index, - int direction, - const std::vector& nonbasic_mark, - std::vector& delta_z_mark, - std::vector& delta_z_indices, - std::vector& delta_z, - double& work_estimate); +template void compute_reduced_cost_update(const lp_problem_t& lp, + const std::vector& basic_list, + const std::vector& nonbasic_list, + const std::vector& delta_y, + int leaving_index, + int direction, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + double& work_estimate); + +template void compute_delta_z(const csc_matrix_t& A_transpose, + const sparse_vector_t& delta_y, + int leaving_index, + int direction, + const std::vector& nonbasic_mark, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + double& work_estimate); #endif } // namespace cuopt::linear_programming::dual_simplex From a3193c94a07dae6bb45715c18d76cb72a101c05a Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Mon, 23 Mar 2026 14:33:11 +0100 Subject: [PATCH 09/19] fixed invalid symbols --- cpp/src/branch_and_bound/pseudo_costs.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index fc9cf79343..04a51033e1 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -124,7 +124,7 @@ objective_change_estimate_t single_pivot_objective_change_estimate( for (i_t j = 0; j < lp.num_cols; j++) { dual_residual[j] -= lp.objective[j]; } - matrix_transpose_vector_multiply(lp.A, 1.0, y, 1.0, dual_residual); + matrix_transpose_vector_multiply(lp.A, 1.0, lp_solution.y, 1.0, dual_residual); f_t dual_residual_norm = vector_norm_inf(dual_residual); settings.log.printf("Dual residual norm: %e\n", dual_residual_norm); } @@ -140,7 +140,7 @@ objective_change_estimate_t single_pivot_objective_change_estimate( #ifdef CHECK_DELTA_OBJ f_t delta_obj_check = 0.0; for (i_t k = 0; k < delta_y.i.size(); k++) { - delta_obj_check += lp.rhs[delta_y.i[k]] * delta_y.lp_solution.x[k]; + delta_obj_check += lp.rhs[delta_y.i[k]] * delta_y.x[k]; } for (i_t h = 0; h < delta_z_indices.size(); h++) { const i_t j = delta_z_indices[h]; From f535b497c91f4e2fba5e994f13076fb749252324 Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Tue, 24 Mar 2026 15:56:46 +0100 Subject: [PATCH 10/19] removed debug code Signed-off-by: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> --- cpp/src/branch_and_bound/branch_and_bound.cpp | 19 ------------------- cpp/src/branch_and_bound/pseudo_costs.cpp | 11 +---------- 2 files changed, 1 insertion(+), 29 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index d6792320d0..1ed6de776d 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -789,25 +789,6 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, upper_bound_ = leaf_objective; report(feasible_solution_symbol(thread_type), leaf_objective, get_lower_bound(), leaf_depth, 0); send_solution = true; - - f_t upper = upper_bound_; - f_t lower = get_lower_bound(); - f_t user_upper = compute_user_objective(original_lp_, upper); - f_t user_lower = compute_user_objective(original_lp_, lower); - - std::cout << std::format( - "{}: user_obj={:.3g}, solver_obj={:.3g}, user_lower={:.3g}, " - "solver_lower={:.3g}, user_gap={:.3g}, " - "solver_gap={:.3g}, tol={:.3g}", - feasible_solution_symbol(thread_type), - user_upper, - upper, - user_lower, - lower, - std::abs(user_upper - user_lower), - std::abs(upper - lower), - settings_.absolute_mip_gap_tol) - << std::endl; } if (send_solution && settings_.solution_callback != nullptr) { diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index 04a51033e1..bf2e78d572 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -901,15 +901,6 @@ i_t pseudo_costs_t::reliable_variable_selection( reliable_threshold); if (unreliable_list.size() > max_num_candidates) { - std::cout << std::format( - "RB iters = {}, B&B iters = {}, unreliable = {}, num_tasks = {}, " - "reliable_threshold = {}\n", - strong_branching_lp_iter.load(), - branch_and_bound_lp_iters, - unreliable_list.size(), - num_tasks, - reliable_threshold); - if (reliability_branching_settings.rank_candidates_with_dual_pivot) { i_t m = worker->leaf_problem.num_rows; i_t n = worker->leaf_problem.num_cols; @@ -971,7 +962,7 @@ i_t pseudo_costs_t::reliable_variable_selection( return branch_var; } -#pragma omp taskloop if (num_candidates > 5) priority(task_priority) num_tasks(num_tasks) \ +#pragma omp taskloop if (num_candidates > 1) priority(task_priority) num_tasks(num_tasks) \ shared(score_mutex) for (i_t i = 0; i < num_candidates; ++i) { auto [score, j] = unreliable_list[i]; From e8dea1e525e969aa53f42827d519dead41d3f467 Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Tue, 24 Mar 2026 16:08:41 +0100 Subject: [PATCH 11/19] enable by default the candidate ranking via simplex pivot Signed-off-by: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> --- cpp/src/branch_and_bound/pseudo_costs.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index e21034f08c..868380a620 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.hpp +++ b/cpp/src/branch_and_bound/pseudo_costs.hpp @@ -407,7 +407,7 @@ struct reliability_branching_settings_t { // Estimate the objective change of each fractional variable // using a single pivot of dual simplex. Then rank the candidates // based on this estimation. - bool rank_candidates_with_dual_pivot = false; + bool rank_candidates_with_dual_pivot = true; }; template From 839d56203585dd901e38856914fd540ce0ad0b78 Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Tue, 24 Mar 2026 16:19:01 +0100 Subject: [PATCH 12/19] remove additional debug code Signed-off-by: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> --- cpp/src/branch_and_bound/branch_and_bound.cpp | 8 -------- cpp/src/math_optimization/solver_settings.cu | 2 +- 2 files changed, 1 insertion(+), 9 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 1ed6de776d..8f952db8c3 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1970,14 +1970,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.log.printf("Reduced cost strengthening enabled: %d\n", settings_.reduced_cost_strengthening); - if (settings_.reliability_branching == -2) { - settings_.log.printf("Enabled objective estimation via dual simplex\n"); - pc_.reliability_branching_settings.rank_candidates_with_dual_pivot = true; - } else { - settings_.log.printf("Disabled objective estimation via dual simplex\n"); - pc_.reliability_branching_settings.rank_candidates_with_dual_pivot = false; - } - variable_bounds_t variable_bounds( original_lp_, settings_, var_types_, Arow_, new_slacks_); diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 3bb10096d5..13a8f74145 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -104,7 +104,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_PRESOLVE, reinterpret_cast(&mip_settings.presolver), CUOPT_PRESOLVE_DEFAULT, CUOPT_PRESOLVE_PSLP, CUOPT_PRESOLVE_DEFAULT}, {CUOPT_MIP_DETERMINISM_MODE, &mip_settings.determinism_mode, CUOPT_MODE_OPPORTUNISTIC, CUOPT_MODE_DETERMINISTIC, CUOPT_MODE_OPPORTUNISTIC}, {CUOPT_RANDOM_SEED, &mip_settings.seed, -1, std::numeric_limits::max(), -1}, - {CUOPT_MIP_RELIABILITY_BRANCHING, &mip_settings.reliability_branching, -2, std::numeric_limits::max(), -1}, + {CUOPT_MIP_RELIABILITY_BRANCHING, &mip_settings.reliability_branching, -1, std::numeric_limits::max(), -1}, {CUOPT_PDLP_PRECISION, reinterpret_cast(&pdlp_settings.pdlp_precision), CUOPT_PDLP_DEFAULT_PRECISION, CUOPT_PDLP_MIXED_PRECISION, CUOPT_PDLP_DEFAULT_PRECISION} }; From 6c68a2b6942ffdf0aa08d69ee8bf938cfdda9b24 Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Wed, 25 Mar 2026 15:45:58 +0100 Subject: [PATCH 13/19] small refactor Signed-off-by: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> --- cpp/src/branch_and_bound/pseudo_costs.cpp | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index bf2e78d572..439dad4aae 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -398,10 +398,8 @@ f_t trial_branching(const lp_problem_t& original_lp, f_t branch_var_lower, f_t branch_var_upper, f_t upper_bound, - i_t bnb_lp_iter_per_node, f_t start_time, - i_t upper_max_lp_iter, - i_t lower_max_lp_iter, + i_t iter_limit, omp_atomic_t& total_lp_iter) { lp_problem_t child_problem = original_lp; @@ -411,9 +409,7 @@ f_t trial_branching(const lp_problem_t& original_lp, const bool initialize_basis = false; simplex_solver_settings_t child_settings = settings; child_settings.set_log(false); - i_t lp_iter_upper = upper_max_lp_iter; - i_t lp_iter_lower = lower_max_lp_iter; - child_settings.iteration_limit = std::clamp(bnb_lp_iter_per_node, lp_iter_lower, lp_iter_upper); + child_settings.iteration_limit = iter_limit; child_settings.inside_mip = 2; child_settings.scale_columns = false; @@ -439,7 +435,7 @@ f_t trial_branching(const lp_problem_t& original_lp, basis_update_mpf_t child_basis_factors = basis_factors; // Only refactor the basis if we encounter numerical issues. - child_basis_factors.set_refactor_frequency(upper_max_lp_iter); + child_basis_factors.set_refactor_frequency(iter_limit); dual::status_t status = dual_phase2_with_advanced_basis(2, 0, @@ -821,11 +817,9 @@ i_t pseudo_costs_t::reliable_variable_selection( i_t num_initialized_up; f_t pseudo_cost_down_avg; f_t pseudo_cost_up_avg; - lp_solution_t& leaf_solution = worker->leaf_solution; initialized(num_initialized_down, num_initialized_up, pseudo_cost_down_avg, pseudo_cost_up_avg); - log.printf("PC: num initialized down %d up %d avg down %e up %e\n", num_initialized_down, num_initialized_up, @@ -835,6 +829,9 @@ i_t pseudo_costs_t::reliable_variable_selection( const int64_t branch_and_bound_lp_iters = bnb_stats.total_lp_iters; const i_t branch_and_bound_lp_iter_per_node = bnb_stats.nodes_explored > 0 ? branch_and_bound_lp_iters / bnb_stats.nodes_explored : 0; + const i_t iter_limit_per_trial = std::clamp(2 * branch_and_bound_lp_iter_per_node, + reliability_branching_settings.lower_max_lp_iter, + reliability_branching_settings.upper_max_lp_iter); i_t reliable_threshold = settings.reliability_branching; if (reliable_threshold < 0) { @@ -984,10 +981,8 @@ i_t pseudo_costs_t::reliable_variable_selection( worker->leaf_problem.lower[j], std::floor(leaf_solution.x[j]), upper_bound, - branch_and_bound_lp_iter_per_node, start_time, - reliability_branching_settings.upper_max_lp_iter, - reliability_branching_settings.lower_max_lp_iter, + iter_limit_per_trial, strong_branching_lp_iter); if (!std::isnan(obj)) { @@ -1015,10 +1010,8 @@ i_t pseudo_costs_t::reliable_variable_selection( std::ceil(leaf_solution.x[j]), worker->leaf_problem.upper[j], upper_bound, - branch_and_bound_lp_iter_per_node, start_time, - reliability_branching_settings.upper_max_lp_iter, - reliability_branching_settings.lower_max_lp_iter, + iter_limit_per_trial, strong_branching_lp_iter); if (!std::isnan(obj)) { From 15604c7e40451fd294b9145a1c6096faf40f98a9 Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Thu, 26 Mar 2026 10:36:14 +0100 Subject: [PATCH 14/19] small refactor Signed-off-by: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> --- cpp/src/branch_and_bound/pseudo_costs.cpp | 60 +++++++++++++---------- 1 file changed, 35 insertions(+), 25 deletions(-) diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index 439dad4aae..570384dfdf 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -809,23 +809,14 @@ i_t pseudo_costs_t::reliable_variable_selection( int max_num_tasks, logger_t& log) { - constexpr f_t eps = 1e-6; - f_t start_time = bnb_stats.start_time; - i_t branch_var = fractional[0]; - f_t max_score = -1; - i_t num_initialized_down; - i_t num_initialized_up; - f_t pseudo_cost_down_avg; - f_t pseudo_cost_up_avg; + constexpr f_t eps = 1e-6; + f_t start_time = bnb_stats.start_time; + i_t branch_var = fractional[0]; + f_t max_score = -1; + f_t pseudo_cost_down_avg = -1; + f_t pseudo_cost_up_avg = -1; lp_solution_t& leaf_solution = worker->leaf_solution; - initialized(num_initialized_down, num_initialized_up, pseudo_cost_down_avg, pseudo_cost_up_avg); - log.printf("PC: num initialized down %d up %d avg down %e up %e\n", - num_initialized_down, - num_initialized_up, - pseudo_cost_down_avg, - pseudo_cost_up_avg); - const int64_t branch_and_bound_lp_iters = bnb_stats.total_lp_iters; const i_t branch_and_bound_lp_iter_per_node = bnb_stats.nodes_explored > 0 ? branch_and_bound_lp_iters / bnb_stats.nodes_explored : 0; @@ -851,22 +842,36 @@ i_t pseudo_costs_t::reliable_variable_selection( reliable_threshold = strong_branching_lp_iter < max_reliability_iter ? reliable_threshold : 0; } + // If `reliable_threshold == 0`, then we set the uninitialized pseudocosts to the average. + // Otherwise, the best ones are initialized via strong branching, while the other are ignored. // + // In the latter, we are not using the average pseudocost (which calculated in the `initialized` + // method). + if (reliable_threshold == 0) { + i_t num_initialized_up; + i_t num_initialized_down; + initialized(num_initialized_down, num_initialized_up, pseudo_cost_down_avg, pseudo_cost_up_avg); + log.printf("PC: num initialized down %d up %d avg down %e up %e\n", + num_initialized_down, + num_initialized_up, + pseudo_cost_down_avg, + pseudo_cost_up_avg); + } + std::vector> unreliable_list; omp_mutex_t score_mutex; for (i_t j : fractional) { + if (pseudo_cost_num_down[j] < reliable_threshold || + pseudo_cost_num_up[j] < reliable_threshold) { + unreliable_list.push_back(std::make_pair(-1, j)); + continue; + } f_t score = calculate_pseudocost_score(j, leaf_solution.x, pseudo_cost_up_avg, pseudo_cost_down_avg); - if (pseudo_cost_num_down[j] < reliable_threshold || - pseudo_cost_num_up[j] < reliable_threshold) { - if (pseudo_cost_num_down[j] == 0 || pseudo_cost_num_up[j] == 0) { score = -1; } - unreliable_list.push_back(std::make_pair(score, j)); - } else { - if (score > max_score) { - max_score = score; - branch_var = j; - } + if (score > max_score) { + max_score = score; + branch_var = j; } } @@ -917,7 +922,8 @@ i_t pseudo_costs_t::reliable_variable_selection( } for (auto& [score, j] : unreliable_list) { - if (score == -1) { + if (pseudo_cost_num_down[j] == 0 || pseudo_cost_num_up[j] == 0) { + // Estimate the objective change by performing a single pivot of dual simplex. objective_change_estimate_t estimate = single_pivot_objective_change_estimate(worker->leaf_problem, settings, @@ -935,6 +941,10 @@ i_t pseudo_costs_t::reliable_variable_selection( work_estimate); score = std::max(estimate.up_obj_change, eps) * std::max(estimate.down_obj_change, eps); + } else { + // Use the previous score, even if it is unreliable + score = calculate_pseudocost_score( + j, leaf_solution.x, pseudo_cost_up_avg, pseudo_cost_down_avg); } } } else { From 9c9d019c4d2ea55a2660d6ef98fb12608926150e Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Mon, 30 Mar 2026 16:16:36 +0200 Subject: [PATCH 15/19] increasing time limit for `seymour1` instance for the determinism test Signed-off-by: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> --- cpp/tests/mip/determinism_test.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/tests/mip/determinism_test.cu b/cpp/tests/mip/determinism_test.cu index 1e59fba649..7519911192 100644 --- a/cpp/tests/mip/determinism_test.cu +++ b/cpp/tests/mip/determinism_test.cu @@ -233,7 +233,7 @@ INSTANTIATE_TEST_SUITE_P( std::make_tuple("/mip/gen-ip054.mps", 128, 120.0, 1), std::make_tuple("/mip/bb_optimality.mps", 4, 60.0, 4), std::make_tuple("/mip/neos5.mps", 16, 60.0, 1), - std::make_tuple("/mip/seymour1.mps", 16, 60.0, 1), + std::make_tuple("/mip/seymour1.mps", 16, 120.0, 1), // too heavy for CI // std::make_tuple("/mip/n2seq36q.mps", 16, 60.0, 4), std::make_tuple("/mip/gmu-35-50.mps", 32, 60.0, 3)), From bb20dd3d031a0fb2ede7d72da54b1df2a8a04770 Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> Date: Wed, 1 Apr 2026 11:17:14 +0200 Subject: [PATCH 16/19] fixed incorrect cut off calculation. added pivot estimation for strong branching as a setting Signed-off-by: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> --- cpp/src/branch_and_bound/pseudo_costs.cpp | 135 +++++++++--------- .../dual_simplex/simplex_solver_settings.hpp | 5 + 2 files changed, 73 insertions(+), 67 deletions(-) diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index 570384dfdf..9eedff0a0d 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -252,6 +252,27 @@ void initialize_pseudo_costs_with_estimate(const lp_problem_t& lp, } } +template +f_t objective_upper_bound(const lp_problem_t& lp, f_t upper_bound, f_t dual_tol) +{ + f_t cut_off = 0; + + if (std::isfinite(upper_bound)) { + cut_off = upper_bound + dual_tol; + } else { + cut_off = 0; + for (i_t j = 0; j < lp.num_cols; ++j) { + if (lp.objective[j] > 0) { + cut_off += lp.objective[j] * lp.upper[j]; + } else if (lp.objective[j] < 0) { + cut_off += lp.objective[j] * lp.lower[j]; + } + } + } + + return cut_off; +} + template void strong_branch_helper(i_t start, i_t end, @@ -293,19 +314,8 @@ void strong_branch_helper(i_t start, if (elapsed_time > settings.time_limit) { break; } child_settings.time_limit = std::max(0.0, settings.time_limit - elapsed_time); child_settings.iteration_limit = 200; - - if (std::isfinite(upper_bound)) { - child_settings.cut_off = upper_bound + settings.dual_tol; - } else { - child_settings.cut_off = 0; - for (i_t i = 0; i < original_lp.num_cols; ++i) { - if (original_lp.objective[i] < 0) { - child_settings.cut_off += original_lp.objective[i] * child_problem.upper[i]; - } else if (original_lp.objective[i] > 0) { - child_settings.cut_off += original_lp.objective[i] * child_problem.lower[i]; - } - } - } + child_settings.cut_off = + objective_upper_bound(child_problem, upper_bound, child_settings.dual_tol); lp_solution_t solution(original_lp.num_rows, original_lp.num_cols); i_t iter = 0; @@ -412,19 +422,8 @@ f_t trial_branching(const lp_problem_t& original_lp, child_settings.iteration_limit = iter_limit; child_settings.inside_mip = 2; child_settings.scale_columns = false; - - if (std::isfinite(upper_bound)) { - child_settings.cut_off = upper_bound + settings.dual_tol; - } else { - child_settings.cut_off = 0; - for (i_t i = 0; i < original_lp.num_cols; ++i) { - if (original_lp.objective[i] < 0) { - child_settings.cut_off += original_lp.objective[i] * child_problem.upper[i]; - } else if (original_lp.objective[i] > 0) { - child_settings.cut_off += original_lp.objective[i] * child_problem.lower[i]; - } - } - } + child_settings.cut_off = + objective_upper_bound(child_problem, upper_bound, child_settings.dual_tol); lp_solution_t solution(original_lp.num_rows, original_lp.num_cols); i_t iter = 0; @@ -652,60 +651,62 @@ void strong_branching(const user_problem_t& original_problem, pc.strong_branch_down[k] = obj_down - root_obj; pc.strong_branch_up[k] = obj_up - root_obj; } - } else if (settings.mip_batch_pdlp_strong_branching == 2) { - initialize_pseudo_costs_with_estimate(original_lp, - settings, - root_vstatus, - root_solution, - basic_list, - nonbasic_list, - fractional, - basis_factors, - pc); - } else { settings.log.printf("Strong branching using %d threads and %ld fractional variables\n", settings.num_threads, fractional.size()); f_t strong_branching_start_time = tic(); + if (settings.mip_strong_branching_use_pivot_estimation) { + initialize_pseudo_costs_with_estimate(original_lp, + settings, + root_vstatus, + root_solution, + basic_list, + nonbasic_list, + fractional, + basis_factors, + pc); + } else { #pragma omp parallel num_threads(settings.num_threads) - { - i_t n = std::min(4 * settings.num_threads, fractional.size()); + { + i_t n = std::min(4 * settings.num_threads, fractional.size()); - // Here we are creating more tasks than the number of threads - // such that they can be scheduled dynamically to the threads. + // Here we are creating more tasks than the number of threads + // such that they can be scheduled dynamically to the threads. #pragma omp for schedule(dynamic, 1) - for (i_t k = 0; k < n; k++) { - i_t start = std::floor(k * fractional.size() / n); - i_t end = std::floor((k + 1) * fractional.size() / n); + for (i_t k = 0; k < n; k++) { + i_t start = std::floor(k * fractional.size() / n); + i_t end = std::floor((k + 1) * fractional.size() / n); + + constexpr bool verbose = false; + if (verbose) { + settings.log.printf("Thread id %d task id %d start %d end %d. size %d\n", + omp_get_thread_num(), + k, + start, + end, + end - start); + } - constexpr bool verbose = false; - if (verbose) { - settings.log.printf("Thread id %d task id %d start %d end %d. size %d\n", - omp_get_thread_num(), - k, - start, - end, - end - start); + strong_branch_helper(start, + end, + start_time, + original_lp, + settings, + var_types, + fractional, + root_obj, + upper_bound, + root_solution.x, + root_vstatus, + edge_norms, + pc); } - - strong_branch_helper(start, - end, - start_time, - original_lp, - settings, - var_types, - fractional, - root_obj, - upper_bound, - root_solution.x, - root_vstatus, - edge_norms, - pc); } + settings.log.printf("Strong branching completed in %.2fs\n", + toc(strong_branching_start_time)); } - settings.log.printf("Strong branching completed in %.2fs\n", toc(strong_branching_start_time)); } pc.update_pseudo_costs_from_strong_branching(fractional, root_solution.x); diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index eadd93040c..089fc62f7c 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -190,6 +190,11 @@ struct simplex_solver_settings_t { i_t mip_batch_pdlp_strong_branching{0}; // 0 if not using batch PDLP for strong branching, 1 if // using batch PDLP for strong branching + // Instead of using the full dual simplex for solving each trial branch in strong branching, + // use just a single pivot from dual simplex to estimate the objective change. + // TODO: Need to pass a proper benchmark before it can be enabled. + bool mip_strong_branching_use_pivot_estimation = false; + diving_heuristics_settings_t diving_settings; // Settings for the diving heuristics // Settings for the reliability branching. From 042eba24f59b038910f44887f73df600b5734cfc Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti Date: Thu, 2 Apr 2026 14:29:02 +0200 Subject: [PATCH 17/19] added CLI parameter to control the maximum number of simplex iteration in strong branching Signed-off-by: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> --- .../cuopt/linear_programming/constants.h | 115 +++++++++--------- .../mip/solver_settings.hpp | 15 +-- cpp/src/branch_and_bound/pseudo_costs.cpp | 17 ++- .../dual_simplex/simplex_solver_settings.hpp | 18 +-- cpp/src/math_optimization/solver_settings.cu | 3 +- 5 files changed, 90 insertions(+), 78 deletions(-) diff --git a/cpp/include/cuopt/linear_programming/constants.h b/cpp/include/cuopt/linear_programming/constants.h index e1097455fa..324b573c82 100644 --- a/cpp/include/cuopt/linear_programming/constants.h +++ b/cpp/include/cuopt/linear_programming/constants.h @@ -20,63 +20,64 @@ #define CUOPT_INSTANTIATE_INT64 0 /* @brief LP/MIP parameter string constants */ -#define CUOPT_ABSOLUTE_DUAL_TOLERANCE "absolute_dual_tolerance" -#define CUOPT_RELATIVE_DUAL_TOLERANCE "relative_dual_tolerance" -#define CUOPT_ABSOLUTE_PRIMAL_TOLERANCE "absolute_primal_tolerance" -#define CUOPT_RELATIVE_PRIMAL_TOLERANCE "relative_primal_tolerance" -#define CUOPT_ABSOLUTE_GAP_TOLERANCE "absolute_gap_tolerance" -#define CUOPT_RELATIVE_GAP_TOLERANCE "relative_gap_tolerance" -#define CUOPT_INFEASIBILITY_DETECTION "infeasibility_detection" -#define CUOPT_STRICT_INFEASIBILITY "strict_infeasibility" -#define CUOPT_PRIMAL_INFEASIBLE_TOLERANCE "primal_infeasible_tolerance" -#define CUOPT_DUAL_INFEASIBLE_TOLERANCE "dual_infeasible_tolerance" -#define CUOPT_ITERATION_LIMIT "iteration_limit" -#define CUOPT_TIME_LIMIT "time_limit" -#define CUOPT_WORK_LIMIT "work_limit" -#define CUOPT_PDLP_SOLVER_MODE "pdlp_solver_mode" -#define CUOPT_METHOD "method" -#define CUOPT_PER_CONSTRAINT_RESIDUAL "per_constraint_residual" -#define CUOPT_SAVE_BEST_PRIMAL_SO_FAR "save_best_primal_so_far" -#define CUOPT_FIRST_PRIMAL_FEASIBLE "first_primal_feasible" -#define CUOPT_LOG_FILE "log_file" -#define CUOPT_LOG_TO_CONSOLE "log_to_console" -#define CUOPT_CROSSOVER "crossover" -#define CUOPT_FOLDING "folding" -#define CUOPT_AUGMENTED "augmented" -#define CUOPT_DUALIZE "dualize" -#define CUOPT_ORDERING "ordering" -#define CUOPT_BARRIER_DUAL_INITIAL_POINT "barrier_dual_initial_point" -#define CUOPT_ELIMINATE_DENSE_COLUMNS "eliminate_dense_columns" -#define CUOPT_CUDSS_DETERMINISTIC "cudss_deterministic" -#define CUOPT_PRESOLVE "presolve" -#define CUOPT_DUAL_POSTSOLVE "dual_postsolve" -#define CUOPT_MIP_DETERMINISM_MODE "mip_determinism_mode" -#define CUOPT_MIP_ABSOLUTE_TOLERANCE "mip_absolute_tolerance" -#define CUOPT_MIP_RELATIVE_TOLERANCE "mip_relative_tolerance" -#define CUOPT_MIP_INTEGRALITY_TOLERANCE "mip_integrality_tolerance" -#define CUOPT_MIP_ABSOLUTE_GAP "mip_absolute_gap" -#define CUOPT_MIP_RELATIVE_GAP "mip_relative_gap" -#define CUOPT_MIP_HEURISTICS_ONLY "mip_heuristics_only" -#define CUOPT_MIP_SCALING "mip_scaling" -#define CUOPT_MIP_PRESOLVE "mip_presolve" -#define CUOPT_MIP_RELIABILITY_BRANCHING "mip_reliability_branching" -#define CUOPT_MIP_CUT_PASSES "mip_cut_passes" -#define CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS "mip_mixed_integer_rounding_cuts" -#define CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS "mip_mixed_integer_gomory_cuts" -#define CUOPT_MIP_KNAPSACK_CUTS "mip_knapsack_cuts" -#define CUOPT_MIP_CLIQUE_CUTS "mip_clique_cuts" -#define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS "mip_strong_chvatal_gomory_cuts" -#define CUOPT_MIP_REDUCED_COST_STRENGTHENING "mip_reduced_cost_strengthening" -#define CUOPT_MIP_CUT_CHANGE_THRESHOLD "mip_cut_change_threshold" -#define CUOPT_MIP_CUT_MIN_ORTHOGONALITY "mip_cut_min_orthogonality" -#define CUOPT_MIP_BATCH_PDLP_STRONG_BRANCHING "mip_batch_pdlp_strong_branching" -#define CUOPT_SOLUTION_FILE "solution_file" -#define CUOPT_NUM_CPU_THREADS "num_cpu_threads" -#define CUOPT_NUM_GPUS "num_gpus" -#define CUOPT_USER_PROBLEM_FILE "user_problem_file" -#define CUOPT_PRESOLVE_FILE "presolve_file" -#define CUOPT_RANDOM_SEED "random_seed" -#define CUOPT_PDLP_PRECISION "pdlp_precision" +#define CUOPT_ABSOLUTE_DUAL_TOLERANCE "absolute_dual_tolerance" +#define CUOPT_RELATIVE_DUAL_TOLERANCE "relative_dual_tolerance" +#define CUOPT_ABSOLUTE_PRIMAL_TOLERANCE "absolute_primal_tolerance" +#define CUOPT_RELATIVE_PRIMAL_TOLERANCE "relative_primal_tolerance" +#define CUOPT_ABSOLUTE_GAP_TOLERANCE "absolute_gap_tolerance" +#define CUOPT_RELATIVE_GAP_TOLERANCE "relative_gap_tolerance" +#define CUOPT_INFEASIBILITY_DETECTION "infeasibility_detection" +#define CUOPT_STRICT_INFEASIBILITY "strict_infeasibility" +#define CUOPT_PRIMAL_INFEASIBLE_TOLERANCE "primal_infeasible_tolerance" +#define CUOPT_DUAL_INFEASIBLE_TOLERANCE "dual_infeasible_tolerance" +#define CUOPT_ITERATION_LIMIT "iteration_limit" +#define CUOPT_TIME_LIMIT "time_limit" +#define CUOPT_WORK_LIMIT "work_limit" +#define CUOPT_PDLP_SOLVER_MODE "pdlp_solver_mode" +#define CUOPT_METHOD "method" +#define CUOPT_PER_CONSTRAINT_RESIDUAL "per_constraint_residual" +#define CUOPT_SAVE_BEST_PRIMAL_SO_FAR "save_best_primal_so_far" +#define CUOPT_FIRST_PRIMAL_FEASIBLE "first_primal_feasible" +#define CUOPT_LOG_FILE "log_file" +#define CUOPT_LOG_TO_CONSOLE "log_to_console" +#define CUOPT_CROSSOVER "crossover" +#define CUOPT_FOLDING "folding" +#define CUOPT_AUGMENTED "augmented" +#define CUOPT_DUALIZE "dualize" +#define CUOPT_ORDERING "ordering" +#define CUOPT_BARRIER_DUAL_INITIAL_POINT "barrier_dual_initial_point" +#define CUOPT_ELIMINATE_DENSE_COLUMNS "eliminate_dense_columns" +#define CUOPT_CUDSS_DETERMINISTIC "cudss_deterministic" +#define CUOPT_PRESOLVE "presolve" +#define CUOPT_DUAL_POSTSOLVE "dual_postsolve" +#define CUOPT_MIP_DETERMINISM_MODE "mip_determinism_mode" +#define CUOPT_MIP_ABSOLUTE_TOLERANCE "mip_absolute_tolerance" +#define CUOPT_MIP_RELATIVE_TOLERANCE "mip_relative_tolerance" +#define CUOPT_MIP_INTEGRALITY_TOLERANCE "mip_integrality_tolerance" +#define CUOPT_MIP_ABSOLUTE_GAP "mip_absolute_gap" +#define CUOPT_MIP_RELATIVE_GAP "mip_relative_gap" +#define CUOPT_MIP_HEURISTICS_ONLY "mip_heuristics_only" +#define CUOPT_MIP_SCALING "mip_scaling" +#define CUOPT_MIP_PRESOLVE "mip_presolve" +#define CUOPT_MIP_RELIABILITY_BRANCHING "mip_reliability_branching" +#define CUOPT_MIP_CUT_PASSES "mip_cut_passes" +#define CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS "mip_mixed_integer_rounding_cuts" +#define CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS "mip_mixed_integer_gomory_cuts" +#define CUOPT_MIP_KNAPSACK_CUTS "mip_knapsack_cuts" +#define CUOPT_MIP_CLIQUE_CUTS "mip_clique_cuts" +#define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS "mip_strong_chvatal_gomory_cuts" +#define CUOPT_MIP_REDUCED_COST_STRENGTHENING "mip_reduced_cost_strengthening" +#define CUOPT_MIP_CUT_CHANGE_THRESHOLD "mip_cut_change_threshold" +#define CUOPT_MIP_CUT_MIN_ORTHOGONALITY "mip_cut_min_orthogonality" +#define CUOPT_MIP_BATCH_PDLP_STRONG_BRANCHING "mip_batch_pdlp_strong_branching" +#define CUOPT_MIP_STRONG_BRANCHING_SIMPLEX_ITER_LIMIT "mip_strong_branching_simplex_iter_limit" +#define CUOPT_SOLUTION_FILE "solution_file" +#define CUOPT_NUM_CPU_THREADS "num_cpu_threads" +#define CUOPT_NUM_GPUS "num_gpus" +#define CUOPT_USER_PROBLEM_FILE "user_problem_file" +#define CUOPT_PRESOLVE_FILE "presolve_file" +#define CUOPT_RANDOM_SEED "random_seed" +#define CUOPT_PDLP_PRECISION "pdlp_precision" /* @brief MIP determinism mode constants */ #define CUOPT_MODE_OPPORTUNISTIC 0 diff --git a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp index 8f4a6c1d78..7516f3e2c8 100644 --- a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp @@ -94,13 +94,14 @@ class mip_solver_settings_t { i_t mixed_integer_gomory_cuts = -1; i_t knapsack_cuts = -1; i_t clique_cuts = -1; - i_t strong_chvatal_gomory_cuts = -1; - i_t reduced_cost_strengthening = -1; - f_t cut_change_threshold = -1.0; - f_t cut_min_orthogonality = 0.5; - i_t mip_batch_pdlp_strong_branching = 0; - i_t num_gpus = 1; - bool log_to_console = true; + i_t strong_chvatal_gomory_cuts = -1; + i_t reduced_cost_strengthening = -1; + f_t cut_change_threshold = -1.0; + f_t cut_min_orthogonality = 0.5; + i_t mip_batch_pdlp_strong_branching = 0; + i_t strong_branching_simplex_iteration_limit = -1; + i_t num_gpus = 1; + bool log_to_console = true; std::string log_file; std::string sol_file; diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index 9eedff0a0d..a8f59e3f35 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -281,11 +281,12 @@ void strong_branch_helper(i_t start, const simplex_solver_settings_t& settings, const std::vector& var_types, const std::vector& fractional, - f_t root_obj, - f_t upper_bound, const std::vector& root_soln, const std::vector& root_vstatus, const std::vector& edge_norms, + f_t root_obj, + f_t upper_bound, + i_t iter_limit, pseudo_costs_t& pc) { raft::common::nvtx::range scope("BB::strong_branch_helper"); @@ -313,7 +314,7 @@ void strong_branch_helper(i_t start, f_t elapsed_time = toc(start_time); if (elapsed_time > settings.time_limit) { break; } child_settings.time_limit = std::max(0.0, settings.time_limit - elapsed_time); - child_settings.iteration_limit = 200; + child_settings.iteration_limit = iter_limit; child_settings.cut_off = objective_upper_bound(child_problem, upper_bound, child_settings.dual_tol); @@ -656,8 +657,11 @@ void strong_branching(const user_problem_t& original_problem, settings.num_threads, fractional.size()); f_t strong_branching_start_time = tic(); + i_t simplex_iteration_limit = settings.strong_branching_simplex_iteration_limit < 0 + ? 200 + : settings.strong_branching_simplex_iteration_limit; - if (settings.mip_strong_branching_use_pivot_estimation) { + if (simplex_iteration_limit < 1) { initialize_pseudo_costs_with_estimate(original_lp, settings, root_vstatus, @@ -696,11 +700,12 @@ void strong_branching(const user_problem_t& original_problem, settings, var_types, fractional, - root_obj, - upper_bound, root_solution.x, root_vstatus, edge_norms, + root_obj, + upper_bound, + simplex_iteration_limit, pc); } } diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 089fc62f7c..5acdffc2ea 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -105,6 +105,8 @@ struct simplex_solver_settings_t { reduced_cost_strengthening(-1), cut_change_threshold(1e-3), cut_min_orthogonality(0.5), + mip_batch_pdlp_strong_branching(0), + strong_branching_simplex_iteration_limit(-1), random_seed(0), reliability_branching(-1), inside_mip(0), @@ -187,13 +189,15 @@ struct simplex_solver_settings_t { // strengthening f_t cut_change_threshold; // threshold for cut change f_t cut_min_orthogonality; // minimum orthogonality for cuts - i_t mip_batch_pdlp_strong_branching{0}; // 0 if not using batch PDLP for strong branching, 1 if - // using batch PDLP for strong branching - - // Instead of using the full dual simplex for solving each trial branch in strong branching, - // use just a single pivot from dual simplex to estimate the objective change. - // TODO: Need to pass a proper benchmark before it can be enabled. - bool mip_strong_branching_use_pivot_estimation = false; + i_t mip_batch_pdlp_strong_branching; // 0 if not using batch PDLP for strong branching, 1 if + // using batch PDLP for strong branching + + // Set the maximum number of simplex iterations allowed per trial branch when applying + // strong branching to the root node. + // -1 - Automatic (iteration limit = 200) + // 0, 1 - Estimate the objective change using a single pivot of dual simplex + // >1 - Set as the iteration limit in dual simplex + i_t strong_branching_simplex_iteration_limit; diving_heuristics_settings_t diving_settings; // Settings for the diving heuristics diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index ba9a46cb3e..a1266836ac 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -99,7 +99,8 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, {CUOPT_NUM_GPUS, &pdlp_settings.num_gpus, 1, 2, 1}, {CUOPT_NUM_GPUS, &mip_settings.num_gpus, 1, 2, 1}, - {CUOPT_MIP_BATCH_PDLP_STRONG_BRANCHING, &mip_settings.mip_batch_pdlp_strong_branching, 0, 2, 0}, + {CUOPT_MIP_BATCH_PDLP_STRONG_BRANCHING, &mip_settings.mip_batch_pdlp_strong_branching, 0, 1, 0}, + {CUOPT_MIP_STRONG_BRANCHING_SIMPLEX_ITER_LIMIT, &mip_settings.strong_branching_simplex_iteration_limit, -1,std::numeric_limits::max(), -1}, {CUOPT_PRESOLVE, reinterpret_cast(&pdlp_settings.presolver), CUOPT_PRESOLVE_DEFAULT, CUOPT_PRESOLVE_PSLP, CUOPT_PRESOLVE_DEFAULT}, {CUOPT_PRESOLVE, reinterpret_cast(&mip_settings.presolver), CUOPT_PRESOLVE_DEFAULT, CUOPT_PRESOLVE_PSLP, CUOPT_PRESOLVE_DEFAULT}, {CUOPT_MIP_DETERMINISM_MODE, &mip_settings.determinism_mode, CUOPT_MODE_OPPORTUNISTIC, CUOPT_MODE_DETERMINISTIC, CUOPT_MODE_OPPORTUNISTIC}, From fa9e90231551fc48498a96eb2129d97e04a07fb5 Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti Date: Thu, 2 Apr 2026 14:39:37 +0200 Subject: [PATCH 18/19] set simplex iteration limits to RINS and SubMIP Signed-off-by: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> --- cpp/src/branch_and_bound/pseudo_costs.cpp | 4 +--- cpp/src/mip_heuristics/diversity/lns/rins.cu | 13 +++++++------ .../diversity/recombiners/sub_mip.cuh | 11 ++++++----- cpp/src/mip_heuristics/solver.cu | 8 ++++++-- 4 files changed, 20 insertions(+), 16 deletions(-) diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index a8f59e3f35..52efe96b98 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -657,9 +657,7 @@ void strong_branching(const user_problem_t& original_problem, settings.num_threads, fractional.size()); f_t strong_branching_start_time = tic(); - i_t simplex_iteration_limit = settings.strong_branching_simplex_iteration_limit < 0 - ? 200 - : settings.strong_branching_simplex_iteration_limit; + i_t simplex_iteration_limit = settings.strong_branching_simplex_iteration_limit; if (simplex_iteration_limit < 1) { initialize_pseudo_costs_with_estimate(original_lp, diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu index d7d7601014..07d8a83a4d 100644 --- a/cpp/src/mip_heuristics/diversity/lns/rins.cu +++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu @@ -264,12 +264,13 @@ void rins_t::run_rins() std::min(current_mip_gap, (f_t)settings.target_mip_gap); branch_and_bound_settings.integer_tol = context.settings.tolerances.integrality_tolerance; branch_and_bound_settings.num_threads = 1; - branch_and_bound_settings.reliability_branching = 0; - branch_and_bound_settings.max_cut_passes = 0; - branch_and_bound_settings.clique_cuts = 0; - branch_and_bound_settings.sub_mip = 1; - branch_and_bound_settings.log.log = false; - branch_and_bound_settings.log.log_prefix = "[RINS] "; + branch_and_bound_settings.reliability_branching = 0; + branch_and_bound_settings.max_cut_passes = 0; + branch_and_bound_settings.clique_cuts = 0; + branch_and_bound_settings.sub_mip = 1; + branch_and_bound_settings.strong_branching_simplex_iteration_limit = 200; + branch_and_bound_settings.log.log = false; + branch_and_bound_settings.log.log_prefix = "[RINS] "; branch_and_bound_settings.solution_callback = [&rins_solution_queue](std::vector& solution, f_t objective) { rins_solution_queue.push_back(solution); diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh index a867141d0a..2a7bc3ff99 100644 --- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh @@ -106,11 +106,12 @@ class sub_mip_recombiner_t : public recombiner_t { branch_and_bound_settings.relative_mip_gap_tol = context.settings.tolerances.relative_mip_gap; branch_and_bound_settings.integer_tol = context.settings.tolerances.integrality_tolerance; branch_and_bound_settings.num_threads = 1; - branch_and_bound_settings.reliability_branching = 0; - branch_and_bound_settings.max_cut_passes = 0; - branch_and_bound_settings.clique_cuts = 0; - branch_and_bound_settings.sub_mip = 1; - branch_and_bound_settings.solution_callback = [this](std::vector& solution, + branch_and_bound_settings.reliability_branching = 0; + branch_and_bound_settings.max_cut_passes = 0; + branch_and_bound_settings.clique_cuts = 0; + branch_and_bound_settings.sub_mip = 1; + branch_and_bound_settings.strong_branching_simplex_iteration_limit = 200; + branch_and_bound_settings.solution_callback = [this](std::vector& solution, f_t objective) { this->solution_callback(solution, objective); }; diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index ee27c79052..02b60c1985 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -241,12 +241,16 @@ solution_t mip_solver_t::run_solver() branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; branch_and_bound_settings.strong_chvatal_gomory_cuts = context.settings.strong_chvatal_gomory_cuts; - branch_and_bound_settings.reduced_cost_strengthening = - context.settings.reduced_cost_strengthening; branch_and_bound_settings.cut_change_threshold = context.settings.cut_change_threshold; branch_and_bound_settings.cut_min_orthogonality = context.settings.cut_min_orthogonality; branch_and_bound_settings.mip_batch_pdlp_strong_branching = context.settings.mip_batch_pdlp_strong_branching; + + branch_and_bound_settings.strong_branching_simplex_iteration_limit = + context.settings.strong_branching_simplex_iteration_limit < 0 + ? 200 + : context.settings.strong_branching_simplex_iteration_limit; + branch_and_bound_settings.reduced_cost_strengthening = context.settings.reduced_cost_strengthening == -1 ? 2 From 64d58f29889146ac0edf15707c5ae64f5f9fef7c Mon Sep 17 00:00:00 2001 From: Nicolas Guidotti Date: Thu, 2 Apr 2026 19:03:40 +0200 Subject: [PATCH 19/19] decreasing the work limit from 3 to 2 for `gmu-35-50` in the deterministic_test Signed-off-by: Nicolas Guidotti <224634272+nguidotti@users.noreply.github.com> --- cpp/tests/mip/determinism_test.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/tests/mip/determinism_test.cu b/cpp/tests/mip/determinism_test.cu index 7519911192..dcd6f7749d 100644 --- a/cpp/tests/mip/determinism_test.cu +++ b/cpp/tests/mip/determinism_test.cu @@ -236,7 +236,7 @@ INSTANTIATE_TEST_SUITE_P( std::make_tuple("/mip/seymour1.mps", 16, 120.0, 1), // too heavy for CI // std::make_tuple("/mip/n2seq36q.mps", 16, 60.0, 4), - std::make_tuple("/mip/gmu-35-50.mps", 32, 60.0, 3)), + std::make_tuple("/mip/gmu-35-50.mps", 32, 60.0, 2)), [](const ::testing::TestParamInfo& info) { const auto& path = std::get<0>(info.param); int threads = std::get<1>(info.param);