Skip to content

Use a single dual simplex pivot to initialize pseudocost #963

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

Open
chris-maes wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Use a single dual simplex pivot to initialize pseudocost #963

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
5 changes: 5 additions & 0 deletions cpp/src/branch_and_bound/branch_and_bound.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2454,10 +2454,15 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_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_);
}

Expand Down
286 changes: 280 additions & 6 deletions cpp/src/branch_and_bound/pseudo_costs.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,256 @@ namespace cuopt::linear_programming::dual_simplex {

namespace {


template <typename f_t>
struct objective_estimate_t {
f_t down_obj;
f_t up_obj;
};

template <typename i_t, typename f_t>
f_t compute_step_length(const simplex_solver_settings_t<i_t, f_t>& settings,
const std::vector<variable_status_t>& vstatus,
const std::vector<f_t>& z,
const std::vector<f_t>& delta_z,
const std::vector<i_t>& 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 <typename i_t, typename f_t>
objective_estimate_t<f_t> single_pivot_objective_estimate(
const lp_problem_t<i_t, f_t>& lp,
const simplex_solver_settings_t<i_t, f_t>& settings,
const csc_matrix_t<i_t, f_t>& A_transpose,
const std::vector<variable_status_t>& vstatus,
i_t variable_j,
i_t basic_j,
f_t root_obj,
const std::vector<f_t>& x,
const std::vector<f_t>& y,
const std::vector<f_t>& z,
const std::vector<i_t>& basic_list,
const std::vector<i_t>& nonbasic_list,
const std::vector<i_t>& nonbasic_mark,
basis_update_mpf_t<i_t, f_t>& basis_factors,
std::vector<i_t>& workspace,
std::vector<f_t>& 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<i_t, f_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<i_t, f_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++; }
Copy link
Contributor

@nguidotti nguidotti Mar 18, 2026

Choose a reason for hiding this comment

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

nitpick: replaced 1E-12 with zero_tol from simplex_settings.

}
work_estimate += nz_delta_y;
const f_t delta_y_nz_percentage = delta_y_nz0 / static_cast<f_t>(lp.num_rows) * 100.0;
const bool use_transpose = delta_y_nz_percentage <= 30.0;
std::vector<i_t> 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<f_t> 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<f_t> 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<i_t, f_t>(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<f_t>{.down_obj = down_obj, .up_obj = up_obj};
}


template <typename i_t, typename f_t>
void initialize_pseudo_costs_with_estimate(const lp_problem_t<i_t, f_t>& lp,
const simplex_solver_settings_t<i_t, f_t>& settings,
const std::vector<variable_status_t>& vstatus,
const std::vector<f_t>& x,
const std::vector<f_t>& y,
const std::vector<f_t>& z,
const std::vector<i_t>& basic_list,
const std::vector<i_t>& nonbasic_list,
const std::vector<i_t>& fractional,
f_t root_obj,
basis_update_mpf_t<i_t, f_t>& basis_factors,
pseudo_costs_t<i_t, f_t>& pc)
{
i_t m = lp.num_rows;
i_t n = lp.num_cols;

csc_matrix_t<i_t, f_t> A_transpose(m, n, 0);
Copy link
Contributor

@nguidotti nguidotti Mar 18, 2026 •
edited
Loading

Choose a reason for hiding this comment

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

We are using $\mathbf{A}^T$ in many places (in dual simplex, bound strengthening and here). Can this move to lp_problem_t so it is visible across the entire code?

lp.A.transpose(A_transpose);

std::vector<f_t> delta_z(n, 0);
std::vector<i_t> workspace(n, 0);

f_t work_estimate = 0;

std::vector<i_t> basic_map(n, -1);
for (i_t i = 0; i < m; i++) {
basic_map[basic_list[i]] = i;
}

std::vector<i_t> 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<f_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;
Copy link
Contributor

@nguidotti nguidotti Mar 18, 2026 •
edited
Loading

Choose a reason for hiding this comment

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

For the pseudocost initilization, should be strong_branch_down/up set to the objective change instead of the objective? This matches the other code paths.

pc.strong_branch_up[k] = estimate.up_obj;
}
}




template <typename i_t, typename f_t>
void strong_branch_helper(i_t start,
i_t end,
Expand Down Expand Up @@ -302,11 +552,16 @@ void strong_branching(const user_problem_t<i_t, f_t>& original_problem,
const simplex_solver_settings_t<i_t, f_t>& settings,
f_t start_time,
const std::vector<variable_type_t>& var_types,
const std::vector<f_t> root_soln,
const std::vector<f_t>& root_x,
const std::vector<f_t>& root_y,
const std::vector<f_t>& root_z,
const std::vector<i_t>& fractional,
f_t root_obj,
const std::vector<variable_status_t>& root_vstatus,
const std::vector<f_t>& edge_norms,
const std::vector<i_t>& basic_list,
const std::vector<i_t>& nonbasic_list,
basis_update_mpf_t<i_t, f_t>& basis_factors,
pseudo_costs_t<i_t, f_t>& pc)
{
pc.resize(original_lp.num_cols);
Expand All @@ -317,7 +572,7 @@ void strong_branching(const user_problem_t<i_t, f_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();
Expand All @@ -328,7 +583,7 @@ void strong_branching(const user_problem_t<i_t, f_t>& original_problem,

// Convert the root_soln to the original problem space
std::vector<f_t> 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<f_t> fraction_values;

Expand Down Expand Up @@ -394,6 +649,20 @@ void strong_branching(const user_problem_t<i_t, f_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,
Expand Down Expand Up @@ -429,7 +698,7 @@ void strong_branching(const user_problem_t<i_t, f_t>& original_problem,
var_types,
fractional,
root_obj,
root_soln,
root_x,
root_vstatus,
edge_norms,
pc);
Expand All @@ -438,7 +707,7 @@ void strong_branching(const user_problem_t<i_t, f_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 <typename i_t, typename f_t>
Expand Down Expand Up @@ -781,11 +1050,16 @@ template void strong_branching<int, double>(const user_problem_t<int, double>& o
const simplex_solver_settings_t<int, double>& settings,
double start_time,
const std::vector<variable_type_t>& var_types,
const std::vector<double> root_soln,
const std::vector<double>& root_x,
const std::vector<double>& root_y,
const std::vector<double>& root_z,
const std::vector<int>& fractional,
double root_obj,
const std::vector<variable_status_t>& root_vstatus,
const std::vector<double>& edge_norms,
const std::vector<int>& basic_list,
const std::vector<int>& nonbasic_list,
basis_update_mpf_t<int, double>& basis_factors,
pseudo_costs_t<int, double>& pc);

#endif
Expand Down
7 changes: 6 additions & 1 deletion cpp/src/branch_and_bound/pseudo_costs.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -522,11 +522,16 @@ void strong_branching(const user_problem_t<i_t, f_t>& original_problem,
const simplex_solver_settings_t<i_t, f_t>& settings,
f_t start_time,
const std::vector<variable_type_t>& var_types,
const std::vector<f_t> root_soln,
const std::vector<f_t>& root_x,
const std::vector<f_t>& root_y,
const std::vector<f_t>& root_z,
const std::vector<i_t>& fractional,
f_t root_obj,
const std::vector<variable_status_t>& root_vstatus,
const std::vector<f_t>& edge_norms,
const std::vector<i_t>& basic_list,
const std::vector<i_t>& nonbasic_list,
basis_update_mpf_t<i_t, f_t>& basis_factors,
pseudo_costs_t<i_t, f_t>& pc);

} // namespace cuopt::linear_programming::dual_simplex
Loading