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 8 commits into
base: release/26.04
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
8 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
3ab9ff7
crossover: hoist delta_zN and delta_expanded out of dual push loop
rg20 Mar 16, 2026
fbbdd35
Added review comments
rg20 Mar 16, 2026
bdea7a2
Remove code duplication
rg20 Mar 16, 2026
cc54886
Merge remote-tracking branch 'upstream/main' into improve_crossover_d…
rg20 Mar 16, 2026
3e2352b
keep the dense vector path alive
rg20 Mar 16, 2026
1ecb41f
Merge remote-tracking branch 'upstream/release/26.04' into improve_cr…
rg20 Mar 20, 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
69 changes: 53 additions & 16 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 @@ -387,6 +388,9 @@ i_t dual_push(const lp_problem_t<i_t, f_t>& lp,
std::vector<f_t>& y = solution.y;
const std::vector<f_t>& x = solution.x;
i_t num_pushes = 0;
std::vector<f_t> delta_zN(n - m);
std::vector<f_t> delta_expanded; // workspace for sparse path (delta_y is sparse enough)
std::vector<f_t> delta_y_dense; // workspace for dense path (delta_y is not sparse enough)
while (superbasic_list.size() > 0) {
const i_t s = superbasic_list.back();
const i_t basic_leaving_index = superbasic_list_index.back();
Expand All @@ -401,11 +405,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 @@ -416,16 +418,38 @@ 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]];
// Choose sparse vs dense method by delta_y sparsity (match dual simplex: sparse if <= 30% nnz)
std::fill(delta_zN.begin(), delta_zN.end(), 0.);
const bool use_sparse = (delta_y_sparse.i.size() * 1.0 / m) <= 0.3;

if (use_sparse) {
delta_expanded.resize(n);
std::fill(delta_expanded.begin(), delta_expanded.end(), 0.);
for (i_t nnz_idx = 0; nnz_idx < static_cast<i_t>(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];
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;
}
}
for (i_t k = 0; k < n - m; ++k) {
delta_zN[k] = -delta_expanded[nonbasic_list[k]];
}
} else {
delta_y_sparse.to_dense(delta_y_dense);
for (i_t k = 0; k < n - m; ++k) {
const i_t j = nonbasic_list[k];
f_t dot = 0.0;
const i_t c_start = lp.A.col_start[j];
const i_t c_end = lp.A.col_start[j + 1];
for (i_t p = c_start; p < c_end; ++p) {
dot += lp.A.x[p] * delta_y_dense[lp.A.i[p]];
}
delta_zN[k] = -dot;
}
delta_zN[k] = -dot;
}

i_t entering_index = -1;
Expand All @@ -435,8 +459,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];
}
Comment on lines +462 to 466
Copy link

@coderabbitai coderabbitai bot Mar 16, 2026

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟡 Minor

Add explicit cast for consistency and to avoid signed/unsigned comparison warning.

Line 428 uses static_cast<i_t>(delta_y_sparse.i.size()) for the same pattern, but line 463 omits it. This inconsistency may trigger compiler warnings under strict settings.

Proposed fix 
-    for (i_t nnz_idx = 0; nnz_idx < delta_y_sparse.i.size(); ++nnz_idx) {
+    for (i_t nnz_idx = 0; nnz_idx < static_cast<i_t>(delta_y_sparse.i.size()); ++nnz_idx) {
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
// 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];
}
// Optimized: Only update non-zero elements from sparse representation
for (i_t nnz_idx = 0; nnz_idx < static_cast<i_t>(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];
}
🤖 Prompt for AI Agents 
Verify each finding against the current code and only fix it if needed.

In `@cpp/src/dual_simplex/crossover.cpp` around lines 462 - 466, The loop over
sparse entries in crossover.cpp uses delta_y_sparse.i.size() as the upper bound
without casting, causing a signed/unsigned comparison warning and inconsistency
with the earlier loop; change the loop condition to use
static_cast<i_t>(delta_y_sparse.i.size()) so the type of nnz_idx (i_t) matches
the cast size, keeping the body unchanged (variables: delta_y_sparse, nnz_idx,
i_t, y, step_length).


// z <- z + step_length * delta z
Expand Down Expand Up @@ -725,7 +751,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 +1027,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 +1216,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 @@ -1331,8 +1360,16 @@ crossover_status_t crossover(const lp_problem_t<i_t, f_t>& lp,
basis_update_mpf_t ft(L, U, p, settings.refactor_frequency);
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);
i_t dual_push_status = dual_push(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