diff --git a/cpp/include/cuopt/linear_programming/constants.h b/cpp/include/cuopt/linear_programming/constants.h index 09a138eb54..d5457d54b4 100644 --- a/cpp/include/cuopt/linear_programming/constants.h +++ b/cpp/include/cuopt/linear_programming/constants.h @@ -20,64 +20,65 @@ #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_IMPLIED_BOUND_CUTS "mip_implied_bound_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_IMPLIED_BOUND_CUTS "mip_implied_bound_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 c2cd33df98..617a8942f0 100644 --- a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp @@ -95,13 +95,14 @@ class mip_solver_settings_t { i_t knapsack_cuts = -1; i_t clique_cuts = -1; i_t implied_bound_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/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 7106d8dce3..c81273b0f0 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -840,11 +840,10 @@ 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_, upper_bound_, worker_pool_.num_idle_workers(), log); @@ -1176,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; @@ -2499,6 +2497,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_, @@ -2506,11 +2505,15 @@ 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_, fractional, root_objective_, + upper_bound_, 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..52efe96b98 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -22,6 +22,257 @@ namespace cuopt::linear_programming::dual_simplex { namespace { +template +struct objective_change_estimate_t { + f_t down_obj_change; + f_t up_obj_change; +}; + +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_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, + const lp_solution_t& lp_solution, + 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 + assert(variable_j >= 0); + assert(basic_j >= 0); + + // 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]) > 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; + 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, + direction, + 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, + direction, + workspace, + delta_z_indices, + delta_z, + work_estimate); + } + + // Verify dual feasibility +#ifdef CHECK_DUAL_FEASIBILITY + { + std::vector dual_residual = lp_solution.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, 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); + } +#endif + + // Compute the step-length + 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_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]; + } + 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(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", + delta_obj_check, + delta_obj, + step_length); + } +#endif + + settings.log.debug( + "Down branch %d. Step length: %e. Delta obj: %e. \n", variable_j, step_length, delta_obj_down); + + // 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, lp_solution.z, delta_z, delta_z_indices); + + 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. Delta obj: %e.\n", variable_j, step_length, delta_obj_up); + + 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 {.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 lp_solution_t& lp_solution, + const std::vector& basic_list, + const std::vector& nonbasic_list, + const std::vector& fractional, + basis_update_mpf_t& basis_factors, + pseudo_costs_t& pc) +{ + i_t m = lp.num_rows; + i_t n = lp.num_cols; + + 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]; + assert(j >= 0); + + 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; + } +} + +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, @@ -30,10 +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, 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"); @@ -61,7 +314,10 @@ 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); + lp_solution_t solution(original_lp.num_rows, original_lp.num_cols); i_t iter = 0; std::vector vstatus = root_vstatus; @@ -80,7 +336,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", @@ -152,10 +409,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; @@ -165,12 +420,11 @@ 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.cut_off = upper_bound + settings.dual_tol; + child_settings.iteration_limit = iter_limit; child_settings.inside_mip = 2; child_settings.scale_columns = false; + 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; @@ -181,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, @@ -219,15 +473,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]; @@ -302,11 +556,15 @@ 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 lp_solution_t& root_solution, const std::vector& fractional, f_t root_obj, + f_t upper_bound, 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 +575,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 +586,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_solution.x, original_root_soln_x); std::vector fraction_values; @@ -399,46 +657,62 @@ 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; + + if (simplex_iteration_limit < 1) { + 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); - - 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); + 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); + } + + strong_branch_helper(start, + end, + start_time, + original_lp, + settings, + var_types, + fractional, + root_solution.x, + root_vstatus, + edge_norms, + root_obj, + upper_bound, + simplex_iteration_limit, + pc); } - - strong_branch_helper(start, - end, - start_time, - original_lp, - settings, - var_types, - fractional, - root_obj, - root_soln, - 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_soln); + pc.update_pseudo_costs_from_strong_branching(fractional, root_solution.x); } template @@ -529,38 +803,30 @@ 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 upper_bound, 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; - - 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); + 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; 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; + 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) { @@ -580,17 +846,32 @@ i_t pseudo_costs_t::reliable_variable_selection( reliable_threshold = strong_branching_lp_iter < max_reliability_iter ? reliable_threshold : 0; } - std::vector unreliable_list; + // 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(j); + unreliable_list.push_back(std::make_pair(-1, j)); 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; @@ -599,13 +880,15 @@ 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; } - 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); @@ -623,18 +906,77 @@ 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 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 (auto& [score, j] : unreliable_list) { + 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, + 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 { + // 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 { + 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(low, high); } + } + } + + // 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) { log.printf("Time limit reached"); return branch_var; } -#pragma omp taskloop if (num_tasks > 1) 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) { - const i_t j = unreliable_list[i]; + auto [score, j] = unreliable_list[i]; if (toc(start_time) > settings.time_limit) { continue; } @@ -651,17 +993,15 @@ 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, - 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)) { 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]++; } @@ -681,18 +1021,16 @@ 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, 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)) { 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]++; } @@ -701,7 +1039,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); + score = + calculate_pseudocost_score(j, leaf_solution.x, pseudo_cost_up_avg, pseudo_cost_down_avg); score_mutex.lock(); if (score > max_score) { @@ -712,7 +1051,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; } @@ -781,11 +1120,15 @@ 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 lp_solution_t& root_solution, const std::vector& fractional, double root_obj, + double upper_bound, 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..868380a620 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,12 @@ 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 upper_bound, int max_num_tasks, logger_t& log); @@ -504,6 +509,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; @@ -522,11 +528,15 @@ 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 lp_solution_t& root_solution, const std::vector& fractional, f_t root_obj, + f_t upper_bound, 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 a6a1b80b6f..34bf0376b7 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -82,6 +82,129 @@ 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,129 +445,6 @@ 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, const csc_matrix_t& A, @@ -2939,30 +2939,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(); @@ -3641,6 +3641,27 @@ 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/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 1ba65aad06..b87fc09ff9 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -106,6 +106,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), @@ -189,8 +191,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 + 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 183c964a90..84d3ca300a 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -101,6 +101,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {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_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}, diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu index 32ffa778d1..63d561cacb 100644 --- a/cpp/src/mip_heuristics/diversity/lns/rins.cu +++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu @@ -265,12 +265,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 f6e4f172cf..8175293b98 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 5d2f043ee4..c40de4c2e4 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -347,12 +347,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 diff --git a/cpp/tests/mip/determinism_test.cu b/cpp/tests/mip/determinism_test.cu index 1e59fba649..dcd6f7749d 100644 --- a/cpp/tests/mip/determinism_test.cu +++ b/cpp/tests/mip/determinism_test.cu @@ -233,10 +233,10 @@ 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)), + 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);