Skip to content

Improve crossover dual simplex #948

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
rg20 wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Improve crossover dual simplex #948

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
00c58df
Optimize right-looking LU factorization with O(1) degree-bucket ops …
rg20 Mar 10, 2026
5f377a2
Take advantage of hyper sparsity in dual push
rg20 Mar 2, 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
40 changes: 26 additions & 14 deletions cpp/src/dual_simplex/crossover.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -331,6 +331,7 @@ void compute_dual_solution_from_basis(const lp_problem_t<i_t, f_t>& lp,

template <typename i_t, typename f_t>
i_t dual_push(const lp_problem_t<i_t, f_t>& lp,
const csr_matrix_t<i_t, f_t>& Arow,
const simplex_solver_settings_t<i_t, f_t>& settings,
f_t start_time,
lp_solution_t<i_t, f_t>& solution,
Expand Down Expand Up @@ -401,11 +402,9 @@ i_t dual_push(const lp_problem_t<i_t, f_t>& lp,
es_sparse.x[0] = -delta_zs;

// B^T delta_y = -delta_zs*es
std::vector<f_t> delta_y(m);
sparse_vector_t<i_t, f_t> delta_y_sparse(m, 1);
sparse_vector_t<i_t, f_t> UTsol_sparse(m, 1);
ft.b_transpose_solve(es_sparse, delta_y_sparse, UTsol_sparse);
delta_y_sparse.scatter(delta_y);

// We solved B^T delta_y = -delta_zs*es, but for the update we need
// U^T*etilde = es.
Expand All @@ -417,15 +416,23 @@ i_t dual_push(const lp_problem_t<i_t, f_t>& lp,

// delta_zN = -N^T delta_y
std::vector<f_t> delta_zN(n - m);
for (i_t k = 0; k < n - m; ++k) {
const i_t j = nonbasic_list[k];
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]];
std::vector<f_t> delta_expanded(n, 0.);

// Iterate directly over sparse delta_y instead of checking zeros
for (i_t nnz_idx = 0; nnz_idx < delta_y_sparse.i.size(); ++nnz_idx) {
const i_t row = delta_y_sparse.i[nnz_idx];
const f_t val = delta_y_sparse.x[nnz_idx];

// Accumulate contributions from this row to all columns
const i_t row_start = Arow.row_start[row];
const i_t row_end = Arow.row_start[row + 1];
for (i_t p = row_start; p < row_end; ++p) {
const i_t col = Arow.j[p];
delta_expanded[col] += Arow.x[p] * val;
}
delta_zN[k] = -dot;
}
for (i_t k = 0; k < n - m; ++k) {
delta_zN[k] = -delta_expanded[nonbasic_list[k]];
}

i_t entering_index = -1;
Expand All @@ -435,8 +442,10 @@ i_t dual_push(const lp_problem_t<i_t, f_t>& lp,
assert(step_length >= -1e-6);

// y <- y + step_length * delta_y
for (i_t i = 0; i < m; ++i) {
y[i] += step_length * delta_y[i];
// Optimized: Only update non-zero elements from sparse representation
for (i_t nnz_idx = 0; nnz_idx < delta_y_sparse.i.size(); ++nnz_idx) {
const i_t i = delta_y_sparse.i[nnz_idx];
y[i] += step_length * delta_y_sparse.x[nnz_idx];
}

// z <- z + step_length * delta z
Expand Down Expand Up @@ -725,7 +734,6 @@ i_t primal_push(const lp_problem_t<i_t, f_t>& lp,
{
const i_t m = lp.num_rows;
const i_t n = lp.num_cols;

settings.log.debug("Primal push: superbasic %ld\n", superbasic_list.size());

std::vector<f_t>& x = solution.x;
Expand Down Expand Up @@ -1002,6 +1010,7 @@ i_t primal_push(const lp_problem_t<i_t, f_t>& lp,
}
solution.x = x_compare;
solution.iterations += num_pushes;

return 0;
}

Expand Down Expand Up @@ -1190,6 +1199,9 @@ crossover_status_t crossover(const lp_problem_t<i_t, f_t>& lp,
f_t crossover_start = tic();
f_t work_estimate = 0;

csr_matrix_t<i_t, f_t> Arow(m, n, 1);
lp.A.to_compressed_row(Arow);

settings.log.printf("\n");
settings.log.printf("Starting crossover\n");

Expand Down Expand Up @@ -1332,7 +1344,7 @@ crossover_status_t crossover(const lp_problem_t<i_t, f_t>& lp,
verify_basis<i_t, f_t>(m, n, vstatus);
compare_vstatus_with_lists<i_t, f_t>(m, n, basic_list, nonbasic_list, vstatus);
i_t dual_push_status = dual_push(
lp, settings, start_time, solution, ft, basic_list, nonbasic_list, superbasic_list, vstatus);
lp, Arow, settings, start_time, solution, ft, basic_list, nonbasic_list, superbasic_list, vstatus);
if (dual_push_status < 0) { return return_to_status(dual_push_status); }
settings.log.debug("basic list size %ld m %d\n", basic_list.size(), m);
settings.log.debug("nonbasic list size %ld n - m %d\n", nonbasic_list.size(), n - m);
Expand Down
Loading
Loading