From c65b3c35eb276c82cdb447cd2a64640450db9984 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Mon, 21 Sep 2026 15:29:12 -0700 Subject: [PATCH 1/3] Recompute dual variables and check dual feasibility after basis repair After basis repair, we recompute dual variables using the working objective. We check the new dual variables for feasibility. If infeasible, we try to restore feasibility by flipping boxed nonbasic variables to their opposite bounds. If recovery fails, we return numerical failure. Callers that cannot recover from basis repairs also propagate numerical failure. We report repaired columns separately from remaining factorization deficiencies. On 34 fast to solve MIPLIB problems, basis repair occurred on 13 problems. In a subsequent run of those 13 problems on this branch, 12,367 basis repairs occurred. Of the 12,327 evaluated by dual-feasibility recovery, 1,972 produced dual infeasibility, 1,509 were recovered by flipping bounds, and 463 remained infeasible. --- cpp/src/branch_and_bound/branch_and_bound.cpp | 14 +- cpp/src/cuts/cuts.cpp | 13 +- cpp/src/cuts/cuts.hpp | 1 + cpp/src/dual_simplex/basis_updates.cpp | 5 +- cpp/src/dual_simplex/basis_updates.hpp | 5 +- cpp/src/dual_simplex/phase2.cpp | 124 ++++++++++++++++-- cpp/src/linear_algebra/vector_math.hpp | 9 ++ 7 files changed, 153 insertions(+), 18 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index a320cc0602..de5cabc2c7 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -3259,16 +3259,18 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( assert(nonbasic_list.size() == original_lp_.num_cols - original_lp_.num_rows); } // Populate the basis_update from the crossover vstatus - i_t refactor_status = basis_update.refactor_basis(original_lp_.A, + i_t deficient_repaired = 0; + i_t refactor_status = basis_update.refactor_basis(original_lp_.A, root_crossover_settings, original_lp_.lower, original_lp_.upper, exploration_stats_.start_time, basic_list, nonbasic_list, - crossover_vstatus_); - if (refactor_status != 0) { - settings_.log.printf("Failed to refactor basis. %d deficient columns.\n", refactor_status); + crossover_vstatus_, + deficient_repaired); + if (refactor_status != 0 || deficient_repaired > 0) { + settings_.log.printf("Failed to refactor basis. %d deficient columns.\n", deficient_repaired); assert(refactor_status == 0); root_status = lp_status_t::NUMERICAL_ISSUES; } @@ -3578,6 +3580,10 @@ auto branch_and_bound_t::do_cut_pass( set_final_solution(solution, root_objective_); return cut_pass_action_t::RETURN; } + if (remove_cuts_status != 0) { + solver_status_ = mip_status_t::NUMERICAL; + return cut_pass_action_t::RETURN; + } f_t remove_cuts_time = toc(remove_cuts_start_time); if (remove_cuts_time > 1.0) { diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index e9f51666dc..ce01ed92c1 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -6559,10 +6559,19 @@ i_t remove_cuts(lp_problem_t& lp, lp.A.col_start[lp.A.n]); basis_update.resize(lp.num_rows); - i_t refactor_status = basis_update.refactor_basis( - lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); + i_t deficient_repaired = 0; + i_t refactor_status = basis_update.refactor_basis(lp.A, + settings, + lp.lower, + lp.upper, + start_time, + basic_list, + nonbasic_list, + vstatus, + deficient_repaired); if (refactor_status == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; } if (refactor_status == TIME_LIMIT_RETURN) { return TIME_LIMIT_RETURN; } + if (refactor_status != 0 || deficient_repaired > 0) { return -1; } } return 0; diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index ca87e26c39..12e8fb1ce9 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -1263,6 +1263,7 @@ i_t add_cuts(const simplex::simplex_solver_settings_t& settings, std::vector& vstatus, std::vector& edge_norms); +// Returns -1 on numerical failure, or the halt/time-limit return code. template i_t remove_cuts(simplex::lp_problem_t& lp, const simplex::simplex_solver_settings_t& settings, diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index f81962d054..6efde71413 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -2346,8 +2346,10 @@ int basis_update_mpf_t::refactor_basis( f_t start_time, std::vector& basic_list, std::vector& nonbasic_list, - std::vector& vstatus) + std::vector& vstatus, + i_t& deficient_repaired) { + deficient_repaired = 0; raft::common::nvtx::range scope("LU::refactor_basis"); std::vector deficient; std::vector slacks_needed; @@ -2374,6 +2376,7 @@ int basis_update_mpf_t::refactor_basis( if (status == TIME_LIMIT_RETURN) { return TIME_LIMIT_RETURN; } if (status == -1) { settings.log.debug("Initial factorization failed\n"); + deficient_repaired = static_cast(deficient.size()); basis_repair(A, settings, lower, diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp index d1c623db55..f5e35b978a 100644 --- a/cpp/src/dual_simplex/basis_updates.hpp +++ b/cpp/src/dual_simplex/basis_updates.hpp @@ -377,7 +377,7 @@ class basis_update_mpf_t { void multiply_lu(csc_matrix_t& out) const; - // Compute L*U = A(p, basic_list) + // Compute L*U = A(p, basic_list). Report the number of deficient columns repaired. int refactor_basis(const csc_matrix_t& A, const simplex_solver_settings_t& settings, const std::vector& lower, @@ -385,7 +385,8 @@ class basis_update_mpf_t { f_t start_time, std::vector& basic_list, std::vector& nonbasic_list, - std::vector& vstatus); + std::vector& vstatus, + i_t& deficient_repaired); void set_refactor_frequency(i_t new_frequency) { refactor_frequency_ = new_frequency; } diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 84656953df..468361d4cb 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -1957,6 +1957,48 @@ f_t dual_infeasibility(const lp_problem_t& lp, return sum_infeasible; } +template +bool recover_dual_feasibility_after_repair(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + basis_update_mpf_t& ft, + const std::vector& basic_list, + const std::vector& nonbasic_list, + std::vector& vstatus, + const std::vector& objective, + std::vector& cB, + std::vector& y, + std::vector& z, + f_t& work_estimate) +{ + for (i_t k = 0; k < lp.num_rows; ++k) { + cB[k] = objective[basic_list[k]]; + } + work_estimate += 3 * lp.num_rows; + ft.b_transpose_solve(cB, y); + compute_reduced_costs(objective, lp.A, y, basic_list, nonbasic_list, z, work_estimate); + work_estimate += lp.num_rows + lp.num_cols; + if (!all_finite(y) || !all_finite(z)) { return false; } + const f_t before = + dual_infeasibility(lp, settings, vstatus, z, settings.tight_tol, settings.dual_tol); + work_estimate += 3 * lp.num_cols; + if (before <= settings.dual_tol) { return true; } + for (i_t j : nonbasic_list) { + if (!std::isfinite(lp.lower[j]) || !std::isfinite(lp.upper[j]) || lp.lower[j] == lp.upper[j]) { + continue; + } + if (vstatus[j] == variable_status_t::NONBASIC_LOWER && z[j] < -settings.dual_tol) { + vstatus[j] = variable_status_t::NONBASIC_UPPER; + } else if (vstatus[j] == variable_status_t::NONBASIC_UPPER && z[j] > settings.dual_tol) { + vstatus[j] = variable_status_t::NONBASIC_LOWER; + } + } + const f_t after = + dual_infeasibility(lp, settings, vstatus, z, settings.tight_tol, settings.dual_tol); + work_estimate += 8 * lp.num_cols; + // The caller must rebuild x and its infeasibilities after these bound changes. + return after <= settings.dual_tol; +} + template f_t primal_infeasibility_breakdown(const lp_problem_t& lp, const simplex_solver_settings_t& settings, @@ -2599,9 +2641,17 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, assert(nonbasic_list.size() == n - m); f_t refactor_start_work = ft.work_estimate(); - i_t refactor_status = ft.refactor_basis( - lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); - refactor_work = ft.work_estimate() - refactor_start_work; + i_t deficient_repaired = 0; + i_t refactor_status = ft.refactor_basis(lp.A, + settings, + lp.lower, + lp.upper, + start_time, + basic_list, + nonbasic_list, + vstatus, + deficient_repaired); + refactor_work = ft.work_estimate() - refactor_start_work; if (refactor_status == CONCURRENT_HALT_RETURN) { return dual_status_t::CONCURRENT_LIMIT; } if (refactor_status == TIME_LIMIT_RETURN) { return dual_status_t::TIME_LIMIT; } if (refactor_status > 0) { return dual_status_t::NUMERICAL; } @@ -2934,8 +2984,16 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, iter, ft.num_updates()); f_t refactor_start_work = ft.work_estimate(); - i_t refactor_status = ft.refactor_basis( - lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); + i_t deficient_repaired = 0; + i_t refactor_status = ft.refactor_basis(lp.A, + settings, + lp.lower, + lp.upper, + start_time, + basic_list, + nonbasic_list, + vstatus, + deficient_repaired); if (refactor_status == CONCURRENT_HALT_RETURN) { return dual_status_t::CONCURRENT_LIMIT; } if (refactor_status == TIME_LIMIT_RETURN) { return dual_status_t::TIME_LIMIT; } if (refactor_status > 0) { return dual_status_t::NUMERICAL; } @@ -2945,6 +3003,20 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, basic_list, nonbasic_list, basic_mark, nonbasic_mark, phase2_work_estimate); compute_initial_nonbasic_end(basic_mark, Arow, nonbasic_end); + if (deficient_repaired > 0 && + !phase2::recover_dual_feasibility_after_repair(lp, + settings, + ft, + basic_list, + nonbasic_list, + vstatus, + objective, + c_basic, + y, + z, + phase2_work_estimate)) { + return dual_status_t::NUMERICAL; + } phase2::compute_primal_solution_from_basis( lp, ft, basic_list, nonbasic_list, vstatus, x, xB_workspace, phase2_work_estimate); @@ -2961,6 +3033,8 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, solve_work = 0.0; if (primal_infeasibility > settings.primal_tol) { + obj = phase2::compute_perturbed_objective(objective, x); + phase2_work_estimate += 2 * n; settings.log.printf( "New infeasibilities found after recompute (primal_inf=%.2e). " "Continuing phase 2.\n", @@ -3591,8 +3665,17 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, num_refactors++; bool should_recompute_x = true; // Needed for numerically difficult problems like cbs-cta f_t refactor_start_work = ft.work_estimate(); - i_t refactor_status = ft.refactor_basis( - lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); + i_t deficient_repaired = 0; + i_t refactor_status = ft.refactor_basis(lp.A, + settings, + lp.lower, + lp.upper, + start_time, + basic_list, + nonbasic_list, + vstatus, + deficient_repaired); + const bool did_basis_repair = deficient_repaired > 0; if (refactor_status == CONCURRENT_HALT_RETURN) { return dual_status_t::CONCURRENT_LIMIT; } if (refactor_status == TIME_LIMIT_RETURN) { return dual_status_t::TIME_LIMIT; } if (refactor_status > 0) { @@ -3602,8 +3685,15 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, i_t count = 0; i_t deficient_size = 0; while (true) { - deficient_size = ft.refactor_basis( - lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); + deficient_size = ft.refactor_basis(lp.A, + settings, + lp.lower, + lp.upper, + start_time, + basic_list, + nonbasic_list, + vstatus, + deficient_repaired); if (deficient_size == CONCURRENT_HALT_RETURN) { return dual_status_t::CONCURRENT_LIMIT; } @@ -3628,6 +3718,20 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, phase2::reset_basis_mark( basic_list, nonbasic_list, basic_mark, nonbasic_mark, phase2_work_estimate); compute_initial_nonbasic_end(basic_mark, Arow, nonbasic_end); + if (did_basis_repair && + !phase2::recover_dual_feasibility_after_repair(lp, + settings, + ft, + basic_list, + nonbasic_list, + vstatus, + objective, + c_basic, + y, + z, + phase2_work_estimate)) { + return dual_status_t::NUMERICAL; + } if (should_recompute_x) { std::vector unperturbed_x(n); phase2_work_estimate += n; @@ -3641,6 +3745,8 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, phase2_work_estimate); x = unperturbed_x; phase2_work_estimate += 2 * n; + obj = phase2::compute_perturbed_objective(objective, x); + phase2_work_estimate += 2 * n; } primal_infeasibility_squared = phase2::compute_initial_primal_infeasibilities(lp, diff --git a/cpp/src/linear_algebra/vector_math.hpp b/cpp/src/linear_algebra/vector_math.hpp index 8063bc735b..6c42416435 100644 --- a/cpp/src/linear_algebra/vector_math.hpp +++ b/cpp/src/linear_algebra/vector_math.hpp @@ -13,6 +13,15 @@ namespace cuopt::mathematical_optimization { +template +bool all_finite(const std::vector& in) +{ + for (f_t value : in) { + if (!std::isfinite(value)) { return false; } + } + return true; +} + // Computes || x ||_inf = max_j | x |_j template f_t vector_norm_inf(const std::vector& x) From 08e5df6c1cbc0e4c0248f0ca63776aaec2188562 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Tue, 22 Sep 2026 10:48:15 -0700 Subject: [PATCH 2/3] Format basis repair logging --- cpp/src/branch_and_bound/branch_and_bound.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index de5cabc2c7..ea3e25c694 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -3270,7 +3270,8 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( crossover_vstatus_, deficient_repaired); if (refactor_status != 0 || deficient_repaired > 0) { - settings_.log.printf("Failed to refactor basis. %d deficient columns.\n", deficient_repaired); + settings_.log.printf("Failed to refactor basis. %d deficient columns.\n", + deficient_repaired); assert(refactor_status == 0); root_status = lp_status_t::NUMERICAL_ISSUES; } From 14e7e3182f13c046652972968011fd95c3a7a9d2 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Tue, 22 Sep 2026 10:51:26 -0700 Subject: [PATCH 3/3] Handle TIME_LIMIT in refactorization --- cpp/src/branch_and_bound/branch_and_bound.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index ea3e25c694..fc40888cf6 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -3269,7 +3269,9 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( nonbasic_list, crossover_vstatus_, deficient_repaired); - if (refactor_status != 0 || deficient_repaired > 0) { + if (refactor_status == TIME_LIMIT_RETURN) { + root_status = lp_status_t::TIME_LIMIT; + } else if (refactor_status != 0 || deficient_repaired > 0) { settings_.log.printf("Failed to refactor basis. %d deficient columns.\n", deficient_repaired); assert(refactor_status == 0);