Skip to content

Adjust basis repair for super basic variables. Repeat root solve if cuts make it infeasible #831

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Merged
rapids-bot merged 6 commits into from
Feb 6, 2026
Merged

Adjust basis repair for super basic variables. Repeat root solve if cuts make it infeasible #831

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
12382b4
Fix issue in crossover where basis repair needs a superbasic slack
chris-maes Feb 5, 2026
a3740e0
Fix assert
chris-maes Feb 5, 2026
65af246
Pop superbasic before basis repair. So that lists are accurate
chris-maes Feb 5, 2026
d19b30f
Fix incorrect classification of cbs-cta as infeasible
chris-maes Feb 5, 2026
48b9e88
Update iterations and objective
chris-maes Feb 5, 2026
4e1d35c
Fix iteration tracking
chris-maes Feb 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
39 changes: 27 additions & 12 deletions cpp/src/dual_simplex/basis_solves.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -619,12 +619,13 @@ i_t basis_repair(const csc_matrix_t<i_t, f_t>& A,
const std::vector<i_t>& slacks_needed,
std::vector<i_t>& basis_list,
std::vector<i_t>& nonbasic_list,
std::vector<i_t>& superbasic_list,
std::vector<variable_status_t>& vstatus)
{
const i_t m = A.m;
const i_t n = A.n;
assert(basis_list.size() == m);
assert(nonbasic_list.size() == n - m);
assert(nonbasic_list.size() + superbasic_list.size() == n - m);

// Create slack_map
std::vector<i_t> slack_map(m); // slack_map[i] = j if column j is e_i
Expand All @@ -649,26 +650,39 @@ i_t basis_repair(const csc_matrix_t<i_t, f_t>& A,
for (i_t k = 0; k < num_nonbasic; ++k) {
nonbasic_map[nonbasic_list[k]] = k;
}
// Create a superbasic_map
std::vector<i_t> superbasic_map(
n, -1); // superbasic_map[j] = p if superbasic[p] = j, -1 if j is basic/nonbasic
const i_t num_superbasic = superbasic_list.size();
for (i_t k = 0; k < num_superbasic; ++k) {
superbasic_map[superbasic_list[k]] = k;
}

const i_t columns_to_replace = deficient.size();
for (i_t k = 0; k < columns_to_replace; ++k) {
const i_t bad_j = basis_list[deficient[k]];
const i_t replace_i = slacks_needed[k];
const i_t replace_j = slack_map[replace_i];
basis_list[deficient[k]] = replace_j;
assert(nonbasic_map[replace_j] != -1);
nonbasic_list[nonbasic_map[replace_j]] = bad_j;
vstatus[replace_j] = variable_status_t::BASIC;
// This is the main issue. What value should bad_j take on.
if (lower[bad_j] == -inf && upper[bad_j] == inf) {
vstatus[bad_j] = variable_status_t::NONBASIC_FREE;
} else if (lower[bad_j] > -inf) {
vstatus[bad_j] = variable_status_t::NONBASIC_LOWER;
} else if (upper[bad_j] < inf) {
vstatus[bad_j] = variable_status_t::NONBASIC_UPPER;
if (nonbasic_map[replace_j] != -1) {
nonbasic_list[nonbasic_map[replace_j]] = bad_j;
// This is the main issue. What value should bad_j take on.
if (lower[bad_j] == -inf && upper[bad_j] == inf) {
vstatus[bad_j] = variable_status_t::NONBASIC_FREE;
} else if (lower[bad_j] > -inf) {
vstatus[bad_j] = variable_status_t::NONBASIC_LOWER;
} else if (upper[bad_j] < inf) {
vstatus[bad_j] = variable_status_t::NONBASIC_UPPER;
} else {
assert(1 == 0);
}
} else if (superbasic_map[replace_j] != -1) {
superbasic_list[superbasic_map[replace_j]] = bad_j;
vstatus[bad_j] = variable_status_t::SUPERBASIC;
} else {
assert(1 == 0);
assert(nonbasic_map[replace_j] != -1 || superbasic_map[replace_j] != -1);
}
vstatus[replace_j] = variable_status_t::BASIC;
}

return 0;
Expand Down Expand Up @@ -865,6 +879,7 @@ template int basis_repair<int, double>(const csc_matrix_t<int, double>& A,
const std::vector<int>& slacks_needed,
std::vector<int>& basis_list,
std::vector<int>& nonbasic_list,
std::vector<int>& superbasic_list,
std::vector<variable_status_t>& vstatus);

template int form_b<int, double>(const csc_matrix_t<int, double>& A,
Expand Down
1 change: 1 addition & 0 deletions cpp/src/dual_simplex/basis_solves.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@ i_t basis_repair(const csc_matrix_t<i_t, f_t>& A,
const std::vector<i_t>& slacks_needed,
std::vector<i_t>& basis_list,
std::vector<i_t>& nonbasic_list,
std::vector<i_t>& superbasic_list,
std::vector<variable_status_t>& vstatus);

// Form the basis matrix B = A(:, basic_list)
Expand Down
13 changes: 11 additions & 2 deletions cpp/src/dual_simplex/basis_updates.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2265,6 +2265,7 @@ int basis_update_mpf_t<i_t, f_t>::refactor_basis(
{
std::vector<i_t> deficient;
std::vector<i_t> slacks_needed;
std::vector<i_t> superbasic_list; // Empty superbasic list

if (L0_.m != A.m) { resize(A.m); }
std::vector<i_t> q;
Expand All @@ -2281,8 +2282,16 @@ int basis_update_mpf_t<i_t, f_t>::refactor_basis(
if (status == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; }
if (status == -1) {
settings.log.debug("Initial factorization failed\n");
basis_repair(
A, settings, lower, upper, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
basis_repair(A,
settings,
lower,
upper,
deficient,
slacks_needed,
basic_list,
nonbasic_list,
superbasic_list,
vstatus);

#ifdef CHECK_BASIS_REPAIR
const i_t m = A.m;
Expand Down
28 changes: 23 additions & 5 deletions cpp/src/dual_simplex/branch_and_bound.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1991,7 +1991,9 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
root_relax_soln_,
iter,
edge_norms_);
f_t dual_phase2_time = toc(dual_phase2_start_time);
exploration_stats_.total_lp_iters += iter;
root_objective_ = compute_objective(original_lp_, root_relax_soln_.x);
f_t dual_phase2_time = toc(dual_phase2_start_time);
if (dual_phase2_time > 1.0) {
settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time);
}
Expand All @@ -2002,11 +2004,27 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
}

if (cut_status != dual::status_t::OPTIMAL) {
settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str());
return mip_status_t::NUMERICAL;
settings_.log.printf("Numerical issue at root node. Resolving from scratch\n");
lp_status_t scratch_status =
solve_linear_program_with_advanced_basis(original_lp_,
exploration_stats_.start_time,
lp_settings,
root_relax_soln_,
basis_update,
basic_list,
nonbasic_list,
root_vstatus_,
edge_norms_);
if (scratch_status == lp_status_t::OPTIMAL) {
// We recovered
cut_status = convert_lp_status_to_dual_status(scratch_status);
exploration_stats_.total_lp_iters += root_relax_soln_.iterations;
root_objective_ = compute_objective(original_lp_, root_relax_soln_.x);
} else {
settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str());
return mip_status_t::NUMERICAL;
}
}
Comment thread
coderabbitai[bot] marked this conversation as resolved.
exploration_stats_.total_lp_iters += root_relax_soln_.iterations;
root_objective_ = compute_objective(original_lp_, root_relax_soln_.x);

local_lower_bounds_.assign(settings_.num_bfs_workers, root_objective_);

Expand Down
98 changes: 60 additions & 38 deletions cpp/src/dual_simplex/crossover.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -516,6 +516,7 @@ i_t dual_push(const lp_problem_t<i_t, f_t>& lp,
slacks_needed,
basic_list,
nonbasic_list,
superbasic_list,
vstatus);
rank =
factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed);
Expand Down Expand Up @@ -585,6 +586,55 @@ f_t primal_residual(const lp_problem_t<i_t, f_t>& lp, const lp_solution_t<i_t, f
return vector_norm_inf<i_t, f_t>(primal_residual);
}

template <typename i_t, typename f_t>
void find_primal_superbasic_variables(const lp_problem_t<i_t, f_t>& lp,
const simplex_solver_settings_t<i_t, f_t>& settings,
const lp_solution_t<i_t, f_t>& initial_solution,
lp_solution_t<i_t, f_t>& solution,
std::vector<variable_status_t>& vstatus,
std::vector<i_t>& nonbasic_list,
std::vector<i_t>& superbasic_list)
{
const i_t n = lp.num_cols;
const f_t fixed_tolerance = settings.fixed_tol;
constexpr f_t basis_threshold = 1e-6;
nonbasic_list.clear();
superbasic_list.clear();

for (i_t j = 0; j < n; ++j) {
if (vstatus[j] != variable_status_t::BASIC) {
const f_t lower_infeas = lp.lower[j] - initial_solution.x[j];
const f_t lower_bound_slack = initial_solution.x[j] - lp.lower[j];
const f_t upper_infeas = initial_solution.x[j] - lp.upper[j];
const f_t upper_bound_slack = lp.upper[j] - initial_solution.x[j];
if (std::abs(lp.lower[j] - lp.upper[j]) < fixed_tolerance) {
vstatus[j] = variable_status_t::NONBASIC_FIXED;
nonbasic_list.push_back(j);
} else if (lower_infeas > 0 && lp.lower[j] > -inf) {
vstatus[j] = variable_status_t::NONBASIC_LOWER;
solution.x[j] = lp.lower[j];
nonbasic_list.push_back(j);
} else if (upper_infeas > 0 && lp.upper[j] < inf) {
vstatus[j] = variable_status_t::NONBASIC_UPPER;
solution.x[j] = lp.upper[j];
nonbasic_list.push_back(j);
} else if (lower_bound_slack < basis_threshold && lp.lower[j] > -inf) {
vstatus[j] = variable_status_t::NONBASIC_LOWER;
nonbasic_list.push_back(j);
} else if (upper_bound_slack < basis_threshold && lp.upper[j] < inf) {
vstatus[j] = variable_status_t::NONBASIC_UPPER;
nonbasic_list.push_back(j);
} else if (lp.lower[j] == -inf && lp.upper[j] == inf) {
vstatus[j] = variable_status_t::NONBASIC_FREE;
nonbasic_list.push_back(j);
} else {
vstatus[j] = variable_status_t::SUPERBASIC;
superbasic_list.push_back(j);
}
}
}
}

template <typename i_t, typename f_t>
f_t primal_ratio_test(const lp_problem_t<i_t, f_t>& lp,
const simplex_solver_settings_t<i_t, f_t>& settings,
Expand Down Expand Up @@ -784,6 +834,7 @@ i_t primal_push(const lp_problem_t<i_t, f_t>& lp,
}
basic_list[basic_leaving_index] = s;
nonbasic_list.push_back(leaving_index);
superbasic_list.pop_back(); // Remove superbasic variable

// Refactor or Update
bool should_refactor = ft.num_updates() > 100;
Expand Down Expand Up @@ -820,7 +871,11 @@ i_t primal_push(const lp_problem_t<i_t, f_t>& lp,
slacks_needed,
basic_list,
nonbasic_list,
superbasic_list,
vstatus);
// We need to be careful. As basis_repair may have changed the superbasic list
find_primal_superbasic_variables(
lp, settings, solution, solution, vstatus, nonbasic_list, superbasic_list);
rank =
factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed);
if (rank == CONCURRENT_HALT_RETURN) {
Expand All @@ -844,11 +899,9 @@ i_t primal_push(const lp_problem_t<i_t, f_t>& lp,
vstatus[s] = variable_status_t::NONBASIC_UPPER;
nonbasic_list.push_back(s);
}
superbasic_list.pop_back(); // Remove superbasic variable
}

// Remove superbasic variable
superbasic_list.pop_back();

num_pushes++;
if (num_pushes % settings.iteration_log_frequency == 0 || toc(last_print_time) > 10.0 ||
superbasic_list.size() == 0) {
Expand Down Expand Up @@ -1178,6 +1231,7 @@ crossover_status_t crossover(const lp_problem_t<i_t, f_t>& lp,
slacks_needed,
basic_list,
nonbasic_list,
superbasic_list,
vstatus);
rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed);
if (rank == CONCURRENT_HALT_RETURN) {
Expand Down Expand Up @@ -1210,41 +1264,8 @@ crossover_status_t crossover(const lp_problem_t<i_t, f_t>& lp,
settings.log.debug("nonbasic list size %ld n - m %d\n", nonbasic_list.size(), n - m);
print_crossover_info(lp, settings, vstatus, solution, "Dual push complete");

nonbasic_list.clear();
superbasic_list.clear();

for (i_t j = 0; j < n; ++j) {
if (vstatus[j] != variable_status_t::BASIC) {
const f_t lower_infeas = lp.lower[j] - initial_solution.x[j];
const f_t lower_bound_slack = initial_solution.x[j] - lp.lower[j];
const f_t upper_infeas = initial_solution.x[j] - lp.upper[j];
const f_t upper_bound_slack = lp.upper[j] - initial_solution.x[j];
if (std::abs(lp.lower[j] - lp.upper[j]) < fixed_tolerance) {
vstatus[j] = variable_status_t::NONBASIC_FIXED;
nonbasic_list.push_back(j);
} else if (lower_infeas > 0 && lp.lower[j] > -inf) {
vstatus[j] = variable_status_t::NONBASIC_LOWER;
solution.x[j] = lp.lower[j];
nonbasic_list.push_back(j);
} else if (upper_infeas > 0 && lp.upper[j] < inf) {
vstatus[j] = variable_status_t::NONBASIC_UPPER;
solution.x[j] = lp.upper[j];
nonbasic_list.push_back(j);
} else if (lower_bound_slack < basis_threshold && lp.lower[j] > -inf) {
vstatus[j] = variable_status_t::NONBASIC_LOWER;
nonbasic_list.push_back(j);
} else if (upper_bound_slack < basis_threshold && lp.upper[j] < inf) {
vstatus[j] = variable_status_t::NONBASIC_UPPER;
nonbasic_list.push_back(j);
} else if (lp.lower[j] == -inf && lp.upper[j] == inf) {
vstatus[j] = variable_status_t::NONBASIC_FREE;
nonbasic_list.push_back(j);
} else {
vstatus[j] = variable_status_t::SUPERBASIC;
superbasic_list.push_back(j);
}
}
}
find_primal_superbasic_variables(
lp, settings, initial_solution, solution, vstatus, nonbasic_list, superbasic_list);

if (superbasic_list.size() > 0) {
std::vector<f_t> save_x = solution.x;
Expand Down Expand Up @@ -1381,6 +1402,7 @@ crossover_status_t crossover(const lp_problem_t<i_t, f_t>& lp,
slacks_needed,
basic_list,
nonbasic_list,
superbasic_list,
vstatus);
rank =
factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed);
Expand Down
1 change: 1 addition & 0 deletions cpp/src/dual_simplex/primal.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -308,6 +308,7 @@ primal::status_t primal_phase2(i_t phase,
slacks_needed,
basic_list,
nonbasic_list,
superbasic_list,
vstatus);
rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed);
if (rank == CONCURRENT_HALT_RETURN) {
Expand Down