From e1b24f2a93137d34fe1ac3b7db0a5808a757f611 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Fri, 14 Aug 2026 15:38:37 -0700 Subject: [PATCH 01/23] Add objective gap termination criteria for QP --- cpp/src/barrier/barrier.cu | 51 +++++++++++-------- cpp/src/barrier/barrier.hpp | 1 + .../dual_simplex/simplex_solver_settings.hpp | 4 ++ cpp/src/pdlp/solve.cu | 1 + 4 files changed, 37 insertions(+), 20 deletions(-) diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index c164296a25..ef346eebc1 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -4012,13 +4012,16 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( f_t& relative_primal_residual, f_t& relative_dual_residual, f_t& relative_complementarity_residual, + f_t& relative_objective_gap, lp_solution_t& solution) { raft::common::nvtx::range fun_scope("Barrier: check_for_suboptimal_solution"); + bool small_gap = (!data.has_cones() && data.Q.n == 0) || + relative_objective_gap < settings.barrier_relaxed_objective_gap_tol; if (relative_primal_residual < settings.barrier_relaxed_feasibility_tol && relative_dual_residual < settings.barrier_relaxed_optimality_tol && relative_complementarity_residual < settings.barrier_relaxed_complementarity_tol && - primal_objective == primal_objective) { + small_gap) { raft::copy(data.x.data(), data.d_x_.data(), data.d_x_.size(), stream_view_); raft::copy(data.y.data(), data.d_y_.data(), data.d_y_.size(), stream_view_); raft::copy(data.z.data(), data.d_z_.data(), data.d_z_.size(), stream_view_); @@ -4220,14 +4223,14 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t upper(lp.upper); data.gather_upper_bounds(upper, data.restrict_u_); @@ -4236,11 +4239,12 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t settings.time_limit) { @@ -4354,6 +4362,7 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t& solution); const simplex::lp_problem_t& lp; diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index c7fe06ed4a..1ff9766e41 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -48,9 +48,11 @@ struct simplex_solver_settings_t { barrier_relative_feasibility_tol(1e-8), barrier_relative_optimality_tol(1e-8), barrier_relative_complementarity_tol(1e-8), + barrier_relative_objective_gap_tol(1e-6), barrier_relaxed_feasibility_tol(1e-4), barrier_relaxed_optimality_tol(1e-4), barrier_relaxed_complementarity_tol(1e-4), + barrier_relaxed_objective_gap_tol(1e-4), cut_off(std::numeric_limits::infinity()), steepest_edge_ratio(0.5), steepest_edge_primal_tol(1e-9), @@ -140,9 +142,11 @@ struct simplex_solver_settings_t { f_t barrier_relative_optimality_tol; // Relative optimality tolerance for barrier method f_t barrier_relative_complementarity_tol; // Relative complementarity tolerance for barrier method + f_t barrier_relative_objective_gap_tol; // Relative objective gap tolerance for barrier method f_t barrier_relaxed_feasibility_tol; // Relative feasibility tolerance for barrier method f_t barrier_relaxed_optimality_tol; // Relative optimality tolerance for barrier method f_t barrier_relaxed_complementarity_tol; // Relative complementarity tolerance for barrier method + f_t barrier_relaxed_objective_gap_tol; // Relative objective gap tolerance for barrier method f_t cut_off; // If the dual objective is greater than the cutoff we stop f_t steepest_edge_ratio; // the ratio of computed steepest edge mismatch from updated steepest edge diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 80b3da2c18..1d6955a961 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -522,6 +522,7 @@ std::tuple, simplex::lp_status_t, f_t, f_t, f_t barrier_settings.barrier_relaxed_feasibility_tol = settings.tolerances.relative_primal_tolerance; barrier_settings.barrier_relaxed_optimality_tol = settings.tolerances.relative_dual_tolerance; barrier_settings.barrier_relaxed_complementarity_tol = settings.tolerances.relative_gap_tolerance; + barrier_settings.barrier_relaxed_objective_gap_tol = settings.tolerances.relative_gap_tolerance; if (barrier_settings.concurrent_halt != nullptr) { // Don't show the barrier log in concurrent mode. Show the PDLP log instead barrier_settings.log.log = false; From 269cbd97a296bc3e755783de719b9697eb97c0d0 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Fri, 14 Aug 2026 18:34:55 -0700 Subject: [PATCH 02/23] Add objective gap reporting and checks in suboptimal cases --- cpp/src/barrier/barrier.cu | 49 ++++++++++++++++++++++++++++--------- cpp/src/barrier/barrier.hpp | 1 + 2 files changed, 38 insertions(+), 12 deletions(-) diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index ef346eebc1..de11f22a16 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -4009,6 +4009,7 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( f_t& primal_residual_norm, f_t& dual_residual_norm, f_t& complementarity_residual_norm, + f_t& objective_gap, f_t& relative_primal_residual, f_t& relative_dual_residual, f_t& relative_complementarity_residual, @@ -4046,12 +4047,16 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( settings.log.printf("Complementarity gap (abs/rel): %8.2e/%8.2e\n", complementarity_residual_norm, relative_complementarity_residual); + settings.log.printf( + "Objective gap (abs/rel): %8.2e/%8.2e\n", objective_gap, relative_objective_gap); settings.log.printf("\n"); return lp_status_t::OPTIMAL; // TODO: Barrier should probably have a separate suboptimal // status } f_t primal_objective_save = data.c.inner_product(data.x_save); + f_t dual_objective_save = + data.b.inner_product(data.y_save) - data.restrict_u_.inner_product(data.v_save); if (data.Q.n > 0) { dense_vector_t Qx_save(data.Q.n); dense_vector_t x_save_host(data.Q.n); @@ -4059,11 +4064,21 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( matrix_vector_multiply(data.Q, 1.0, x_save_host, 0.0, Qx_save); f_t quad_objective = 0.5 * x_save_host.inner_product(Qx_save); primal_objective_save += quad_objective; + dual_objective_save -= quad_objective; } + f_t objective_gap_save = std::abs(primal_objective_save - dual_objective_save); + f_t user_primal_objective_save = compute_user_objective(lp, primal_objective_save); + f_t relative_objective_gap_save = + objective_gap_save / + (1.0 + std::min(std::abs(user_primal_objective_save), std::abs(primal_objective_save))); + bool small_gap_save = (!data.has_cones() && data.Q.n == 0) || + relative_objective_gap_save < settings.barrier_relaxed_objective_gap_tol; + if (data.relative_primal_residual_save < settings.barrier_relaxed_feasibility_tol && data.relative_dual_residual_save < settings.barrier_relaxed_optimality_tol && - data.relative_complementarity_residual_save < settings.barrier_relaxed_complementarity_tol) { + data.relative_complementarity_residual_save < settings.barrier_relaxed_complementarity_tol && + small_gap_save) { settings.log.printf("Restoring previous solution\n"); data.restore_saved_iterate(); data.to_solution(lp, @@ -4086,14 +4101,19 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( settings.log.printf("Complementarity gap (abs/rel): %8.2e/%8.2e\n", data.complementarity_residual_norm_save, data.relative_complementarity_residual_save); + settings.log.printf("Objective gap (abs/rel): %8.2e/%8.2e\n", + objective_gap_save, + relative_objective_gap_save); settings.log.printf("\n"); return lp_status_t::OPTIMAL; // TODO: Barrier should probably have a separate suboptimal // status } else { - settings.log.printf("Primal residual %.2e dual residual %.2e complementarity residual %.2e\n", - relative_primal_residual, - relative_dual_residual, - relative_complementarity_residual); + settings.log.printf( + "Primal residual %.2e dual residual %.2e complementarity residual %.2e objective gap %.2e\n", + relative_primal_residual, + relative_dual_residual, + relative_complementarity_residual, + relative_objective_gap); } settings.log.printf("Search direction computation failed\n"); return lp_status_t::NUMERICAL_ISSUES; @@ -4241,10 +4261,9 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t Date: Mon, 17 Aug 2026 14:50:24 -0700 Subject: [PATCH 03/23] Rename the relaxed objective tolerance --- cpp/src/barrier/barrier.cu | 9 +++++---- cpp/src/dual_simplex/simplex_solver_settings.hpp | 5 +++-- cpp/src/pdlp/solve.cu | 3 ++- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index de11f22a16..cc0c00c36b 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -4018,7 +4018,7 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( { raft::common::nvtx::range fun_scope("Barrier: check_for_suboptimal_solution"); bool small_gap = (!data.has_cones() && data.Q.n == 0) || - relative_objective_gap < settings.barrier_relaxed_objective_gap_tol; + relative_objective_gap < settings.barrier_relaxed_relative_objective_gap_tol; if (relative_primal_residual < settings.barrier_relaxed_feasibility_tol && relative_dual_residual < settings.barrier_relaxed_optimality_tol && relative_complementarity_residual < settings.barrier_relaxed_complementarity_tol && @@ -4072,8 +4072,9 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( f_t relative_objective_gap_save = objective_gap_save / (1.0 + std::min(std::abs(user_primal_objective_save), std::abs(primal_objective_save))); - bool small_gap_save = (!data.has_cones() && data.Q.n == 0) || - relative_objective_gap_save < settings.barrier_relaxed_objective_gap_tol; + bool small_gap_save = + (!data.has_cones() && data.Q.n == 0) || + relative_objective_gap_save < settings.barrier_relaxed_relative_objective_gap_tol; if (data.relative_primal_residual_save < settings.barrier_relaxed_feasibility_tol && data.relative_dual_residual_save < settings.barrier_relaxed_optimality_tol && @@ -4288,7 +4289,7 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t::infinity()), steepest_edge_ratio(0.5), steepest_edge_primal_tol(1e-9), @@ -146,7 +146,8 @@ struct simplex_solver_settings_t { f_t barrier_relaxed_feasibility_tol; // Relative feasibility tolerance for barrier method f_t barrier_relaxed_optimality_tol; // Relative optimality tolerance for barrier method f_t barrier_relaxed_complementarity_tol; // Relative complementarity tolerance for barrier method - f_t barrier_relaxed_objective_gap_tol; // Relative objective gap tolerance for barrier method + f_t barrier_relaxed_relative_objective_gap_tol; // Relative objective gap tolerance for barrier + // method f_t cut_off; // If the dual objective is greater than the cutoff we stop f_t steepest_edge_ratio; // the ratio of computed steepest edge mismatch from updated steepest edge diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 1d6955a961..2ef1ea7874 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -522,7 +522,8 @@ std::tuple, simplex::lp_status_t, f_t, f_t, f_t barrier_settings.barrier_relaxed_feasibility_tol = settings.tolerances.relative_primal_tolerance; barrier_settings.barrier_relaxed_optimality_tol = settings.tolerances.relative_dual_tolerance; barrier_settings.barrier_relaxed_complementarity_tol = settings.tolerances.relative_gap_tolerance; - barrier_settings.barrier_relaxed_objective_gap_tol = settings.tolerances.relative_gap_tolerance; + barrier_settings.barrier_relaxed_relative_objective_gap_tol = + settings.tolerances.relative_gap_tolerance; if (barrier_settings.concurrent_halt != nullptr) { // Don't show the barrier log in concurrent mode. Show the PDLP log instead barrier_settings.log.log = false; From ed622ed91f19ded53d1650ea22e0e89bb232f4bc Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Mon, 17 Aug 2026 15:59:42 -0700 Subject: [PATCH 04/23] Adjust the tolerance to 1e-8 --- cpp/src/dual_simplex/simplex_solver_settings.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index febd7298e7..469518e794 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -48,7 +48,7 @@ struct simplex_solver_settings_t { barrier_relative_feasibility_tol(1e-8), barrier_relative_optimality_tol(1e-8), barrier_relative_complementarity_tol(1e-8), - barrier_relative_objective_gap_tol(1e-6), + barrier_relative_objective_gap_tol(1e-8), barrier_relaxed_feasibility_tol(1e-4), barrier_relaxed_optimality_tol(1e-4), barrier_relaxed_complementarity_tol(1e-4), From 1fedd5aebdebb7d48e8eadb7c6cff30e791653a2 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Tue, 18 Aug 2026 14:34:41 -0700 Subject: [PATCH 05/23] Revert the tolerance --- cpp/src/dual_simplex/simplex_solver_settings.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 469518e794..febd7298e7 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -48,7 +48,7 @@ struct simplex_solver_settings_t { barrier_relative_feasibility_tol(1e-8), barrier_relative_optimality_tol(1e-8), barrier_relative_complementarity_tol(1e-8), - barrier_relative_objective_gap_tol(1e-8), + barrier_relative_objective_gap_tol(1e-6), barrier_relaxed_feasibility_tol(1e-4), barrier_relaxed_optimality_tol(1e-4), barrier_relaxed_complementarity_tol(1e-4), From 6742aae10de7f9fc221329bf417201a136ef15a2 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 26 Aug 2026 11:47:50 -0700 Subject: [PATCH 06/23] Address PR comments --- cpp/src/barrier/barrier.cu | 13 ++++++------- cpp/src/dual_simplex/solve.cpp | 24 ++++++++++++++++++++++++ cpp/src/dual_simplex/solve.hpp | 7 +++++++ 3 files changed, 37 insertions(+), 7 deletions(-) diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index 3273ba545b..0cacec016b 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -4028,7 +4028,7 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( if (relative_primal_residual < settings.barrier_relaxed_feasibility_tol && relative_dual_residual < settings.barrier_relaxed_optimality_tol && relative_complementarity_residual < settings.barrier_relaxed_complementarity_tol && - small_gap) { + small_gap && primal_objective == primal_objective) { raft::copy(data.x.data(), data.d_x_.data(), data.d_x_.size(), stream_view_); raft::copy(data.y.data(), data.d_y_.data(), data.d_y_.size(), stream_view_); raft::copy(data.z.data(), data.d_z_.data(), data.d_z_.size(), stream_view_); @@ -4278,9 +4278,9 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t& lp, f_t obj) return user_obj; } +template +void compute_objective_gap(const lp_problem_t& problem, + f_t primal_obj, + f_t dual_obj, + f_t& objective_gap, + f_t& relative_objective_gap) +{ + objective_gap = std::abs(primal_obj - dual_obj); + f_t user_primal_obj = compute_user_objective(problem, primal_obj); + f_t user_dual_obj = compute_user_objective(problem, dual_obj); + + f_t denom_1 = std::min(std::abs(user_primal_obj), std::abs(primal_obj)); + f_t denom_2 = std::min(std::abs(user_dual_obj), std::abs(dual_obj)); + f_t denom = 1.0 + std::max(denom_1, denom_2); + + relative_objective_gap = objective_gap / denom; +} + template f_t compute_presolved_objective(const lp_problem_t& lp, f_t user_obj) { @@ -819,6 +837,12 @@ template double compute_user_objective(const lp_problem_t& lp, double obj); +template void compute_objective_gap(const lp_problem_t& problem, + double primal_obj, + double dual_obj, + double& objective_gap, + double& relative_objective_gap); + template double compute_presolved_objective(const lp_problem_t& lp, double user_obj); template lp_status_t solve_linear_program_advanced( diff --git a/cpp/src/dual_simplex/solve.hpp b/cpp/src/dual_simplex/solve.hpp index 308c462de5..e150260302 100644 --- a/cpp/src/dual_simplex/solve.hpp +++ b/cpp/src/dual_simplex/solve.hpp @@ -63,6 +63,13 @@ f_t compute_user_objective(const lp_problem_t& lp, const std::vector f_t compute_user_objective(const lp_problem_t& lp, f_t obj); +template +void compute_objective_gap(const lp_problem_t& problem, + f_t primal_obj, + f_t dual_obj, + f_t& objective_gap, + f_t& relative_objective_gap); + template f_t compute_presolved_objective(const lp_problem_t& lp, f_t user_obj); From 51b7dc80a2cc21fb0c21fba849c3c31038bc35a1 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Tue, 15 Sep 2026 23:33:19 -0700 Subject: [PATCH 07/23] Eliminate zero-cost free linear variables by equality substitution in QP/SOCP --- cpp/src/dual_simplex/presolve.cpp | 301 +++++++++++++++++++++++++++ cpp/src/dual_simplex/presolve.hpp | 22 ++ cpp/src/dual_simplex/solve.cpp | 64 +++--- cpp/tests/socp/solve_barrier_socp.cu | 198 +++++++++++++++++- 4 files changed, 548 insertions(+), 37 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index c204e0c798..7a9f8a227b 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -16,6 +16,7 @@ #include #include #include +#include namespace cuopt::mathematical_optimization::simplex { @@ -95,6 +96,240 @@ static void remove_variables_from_Q(csr_matrix_t& Q, Q = std::move(Qout); } +// Eliminate zero-cost, Q-uncoupled free linear variables by sparse equality substitution. +// This is especially useful for conic formulations containing chains of auxiliary free +// variables: unlike regularizing those variables in the KKT system, substitution is exact. +template +static void eliminate_free_variables(lp_problem_t& problem, + presolve_info_t& presolve_info) +{ + const i_t old_m = problem.num_rows; + const i_t old_n = problem.num_cols; + const i_t linear_cols = linear_variable_count(problem); + if (old_m == 0 || linear_cols == 0) { return; } + + std::vector q_present(old_n, 0); + if (problem.Q.n > 0) { + for (i_t row = 0; row < problem.Q.m; ++row) { + for (i_t p = problem.Q.row_start[row]; p < problem.Q.row_start[row + 1]; ++p) { + if (problem.Q.x[p] == 0) { continue; } + q_present[row] = 1; + q_present[problem.Q.j[p]] = 1; + } + } + } + + csr_matrix_t Arow(0, 0, 0); + problem.A.to_compressed_row(Arow); + std::vector> rows(old_m); + std::vector> col_rows(old_n); + for (i_t i = 0; i < old_m; ++i) { + auto& row = rows[i]; + row.reserve(static_cast(Arow.row_start[i + 1] - Arow.row_start[i]) * 2); + for (i_t p = Arow.row_start[i]; p < Arow.row_start[i + 1]; ++p) { + row[Arow.j[p]] = Arow.x[p]; + col_rows[Arow.j[p]].push_back(i); + } + } + + std::vector active_row(old_m, 1); + std::vector active_col(old_n, 1); + auto& eliminations = presolve_info.free_variable_eliminations; + + // col_rows is append-only: fill re-adds a row that may already be listed, and eliminated + // rows are never unlisted. Compact it during the scan so the pass stays linear. + std::vector row_stamp(old_m, -1); + i_t stamp = 0; + std::vector incident; + + // One pass in column order. Peeling a chain from one end keeps each pivot row sparse, so + // revisiting columns buys almost nothing and costs a requeue storm on models with + // hundreds of thousands of free columns. + for (i_t j = 0; j < linear_cols; ++j) { + if (!active_col[j] || problem.lower[j] != -inf || problem.upper[j] != inf || + problem.objective[j] != 0 || q_present[j]) { + continue; + } + + incident.clear(); + auto& listed = col_rows[j]; + size_t keep = 0; + ++stamp; + for (size_t idx = 0; idx < listed.size(); ++idx) { + const i_t i = listed[idx]; + if (row_stamp[i] == stamp) { continue; } + row_stamp[i] = stamp; + if (!active_row[i] || rows[i].find(j) == rows[i].end()) { continue; } + listed[keep++] = i; + incident.push_back(i); + } + listed.resize(keep); + + free_variable_elimination_t elimination; + elimination.variable = j; + if (incident.empty()) { + elimination.pivot_row = -1; + elimination.pivot_coefficient = 1; + elimination.rhs = 0; + eliminations.push_back(std::move(elimination)); + active_col[j] = 0; + continue; + } + + // A small pivot row limits fill. On a path this peels sparse boundary rows instead of + // repeatedly traversing the growing aggregate row. + i_t pivot = incident.front(); + for (const i_t i : incident) { + if (rows[i].size() < rows[pivot].size() || + (rows[i].size() == rows[pivot].size() && + std::abs(rows[i].at(j)) > std::abs(rows[pivot].at(j)))) { + pivot = i; + } + } + const f_t pivot_coefficient = rows[pivot].at(j); + if (pivot_coefficient == 0) { continue; } + + // Accept only if the substitution does not add nonzeros. Dropping the pivot row and the + // a_ij entries pays for the entries the pivot row scatters into the other incident rows. + // Integrator chains have a sparse pivot and degree 2, so they clear this easily; the + // dense equalities of a converted QCQP do not, and aggregating those is what densified A + // and stalled the barrier. + i_t added = 0; + for (const i_t i : incident) { + if (i == pivot) { continue; } + for (const auto& [col, value] : rows[pivot]) { + if (col == j) { continue; } + if (rows[i].find(col) == rows[i].end()) { ++added; } + } + } + const i_t removed = + static_cast(rows[pivot].size()) + static_cast(incident.size()) - 1; + if (added > removed) { continue; } + + elimination.pivot_row = pivot; + elimination.pivot_coefficient = pivot_coefficient; + elimination.rhs = problem.rhs[pivot]; + elimination.columns.reserve(rows[pivot].size() - 1); + elimination.coefficients.reserve(rows[pivot].size() - 1); + for (const auto& [col, value] : rows[pivot]) { + if (col == j) { continue; } + elimination.columns.push_back(col); + elimination.coefficients.push_back(value); + } + + for (const i_t i : incident) { + if (i == pivot) { continue; } + auto j_it = rows[i].find(j); + if (j_it == rows[i].end()) { continue; } + const f_t factor = j_it->second / pivot_coefficient; + rows[i].erase(j_it); + for (size_t k = 0; k < elimination.columns.size(); ++k) { + const i_t col = elimination.columns[k]; + const f_t delta = factor * elimination.coefficients[k]; + auto col_it = rows[i].find(col); + const f_t old_value = col_it == rows[i].end() ? f_t{0} : col_it->second; + const f_t new_value = old_value - delta; + const f_t drop_tol = f_t{100} * std::numeric_limits::epsilon() * + std::max({f_t{1}, std::abs(old_value), std::abs(delta)}); + if (std::abs(new_value) <= drop_tol) { + if (col_it != rows[i].end()) { rows[i].erase(col_it); } + } else { + if (col_it == rows[i].end()) { col_rows[col].push_back(i); } + rows[i][col] = new_value; + } + } + problem.rhs[i] -= factor * elimination.rhs; + elimination.affected_rows.push_back(i); + elimination.factors.push_back(factor); + } + + active_row[pivot] = 0; + active_col[j] = 0; + eliminations.push_back(std::move(elimination)); + } + + if (eliminations.empty()) { return; } + + presolve_info.free_elimination_num_variables = old_n; + presolve_info.free_elimination_num_constraints = old_m; + auto& remaining_cols = presolve_info.free_elimination_remaining_variables; + auto& remaining_rows = presolve_info.free_elimination_remaining_constraints; + remaining_cols.clear(); + remaining_rows.clear(); + + std::vector old_to_new_col(old_n, -1); + for (i_t j = 0; j < old_n; ++j) { + if (!active_col[j]) { continue; } + old_to_new_col[j] = static_cast(remaining_cols.size()); + remaining_cols.push_back(j); + } + for (i_t i = 0; i < old_m; ++i) { + if (active_row[i]) { remaining_rows.push_back(i); } + } + + i_t new_n = static_cast(remaining_cols.size()); + i_t new_m = static_cast(remaining_rows.size()); + i_t new_nnz = 0; + for (const i_t i : remaining_rows) { + for (const auto& [j, value] : rows[i]) { + if (active_col[j] && value != 0) { ++new_nnz; } + } + } + + csr_matrix_t reduced_A(new_m, new_n, new_nnz); + std::vector reduced_rhs(new_m); + i_t nz = 0; + for (i_t new_i = 0; new_i < new_m; ++new_i) { + const i_t old_i = remaining_rows[new_i]; + reduced_A.row_start[new_i] = nz; + std::vector> entries; + entries.reserve(rows[old_i].size()); + for (const auto& [old_j, value] : rows[old_i]) { + if (active_col[old_j] && value != 0) { entries.emplace_back(old_to_new_col[old_j], value); } + } + std::sort(entries.begin(), entries.end()); + for (const auto& [new_j, value] : entries) { + reduced_A.j[nz] = new_j; + reduced_A.x[nz] = value; + ++nz; + } + reduced_rhs[new_i] = problem.rhs[old_i]; + } + reduced_A.row_start[new_m] = nz; + + std::vector objective(new_n); + std::vector lower(new_n); + std::vector upper(new_n); + for (i_t new_j = 0; new_j < new_n; ++new_j) { + const i_t old_j = remaining_cols[new_j]; + objective[new_j] = problem.objective[old_j]; + lower[new_j] = problem.lower[old_j]; + upper[new_j] = problem.upper[old_j]; + } + + std::vector col_marker(old_n, 0); + for (i_t j = 0; j < old_n; ++j) { + if (!active_col[j]) { col_marker[j] = 1; } + } + if (problem.Q.n > 0) { remove_variables_from_Q(problem.Q, col_marker, old_to_new_col, new_n); } + + reduced_A.to_compressed_col(problem.A); + problem.rhs = std::move(reduced_rhs); + problem.objective = std::move(objective); + problem.lower = std::move(lower); + problem.upper = std::move(upper); + problem.num_rows = new_m; + problem.num_cols = new_n; + problem.cone_var_start = old_to_new_col[problem.cone_var_start]; + + presolve_info.direct_free_variables.clear(); + for (i_t new_j = 0; new_j < problem.cone_var_start; ++new_j) { + if (problem.lower[new_j] == -inf && problem.upper[new_j] == inf) { + presolve_info.direct_free_variables.push_back(new_j); + } + } +} + template i_t remove_empty_cols(lp_problem_t& problem, i_t& num_empty_cols, @@ -1497,6 +1732,21 @@ i_t presolve(const lp_problem_t& original, } settings.log.printf("Dependent row check in %.2fs\n", toc(dependent_row_start)); } + + // LP already goes through PSLP; this substitution is for QP/SOCP only. + if (settings.barrier_presolve && (has_cones || problem.Q.n > 0)) { + const i_t old_free_count = static_cast(presolve_info.direct_free_variables.size()); + const f_t free_elimination_start = tic(); + eliminate_free_variables(problem, presolve_info); + const i_t eliminated = + old_free_count - static_cast(presolve_info.direct_free_variables.size()); + if (eliminated > 0) { + settings.log.printf("Eliminated %d free variables by equality substitution in %.2fs\n", + eliminated, + toc(free_elimination_start)); + } + } + assert(problem.num_rows == problem.A.m); assert(problem.num_cols == problem.A.n); if (settings.print_presolve_stats && problem.A.m < original.A.m) { @@ -1769,6 +2019,57 @@ void uncrush_solution(const presolve_info_t& presolve_info, std::vector input_y = crushed_y; std::vector input_z = crushed_z; std::vector free_variable_pairs = presolve_info.free_variable_pairs; + + // Free-variable substitution is the last presolve transformation, so undo it first. + if (!presolve_info.free_variable_eliminations.empty()) { + if (settings.postsolve_info == 1) { + settings.log.printf("Post-solve: Reconstructing %d eliminated free variables\n", + static_cast(presolve_info.free_variable_eliminations.size())); + } + assert(static_cast(input_x.size()) == + static_cast(presolve_info.free_elimination_remaining_variables.size())); + assert(static_cast(input_y.size()) == + static_cast(presolve_info.free_elimination_remaining_constraints.size())); + std::vector expanded_x(presolve_info.free_elimination_num_variables, 0); + std::vector expanded_z(presolve_info.free_elimination_num_variables, 0); + for (i_t k = 0; k < static_cast(presolve_info.free_elimination_remaining_variables.size()); + ++k) { + const i_t j = presolve_info.free_elimination_remaining_variables[k]; + expanded_x[j] = input_x[k]; + expanded_z[j] = input_z[k]; + } + input_x = std::move(expanded_x); + input_z = std::move(expanded_z); + + std::vector expanded_y(presolve_info.free_elimination_num_constraints, 0); + for (i_t k = 0; + k < static_cast(presolve_info.free_elimination_remaining_constraints.size()); + ++k) { + expanded_y[presolve_info.free_elimination_remaining_constraints[k]] = input_y[k]; + } + input_y = std::move(expanded_y); + + for (auto it = presolve_info.free_variable_eliminations.rbegin(); + it != presolve_info.free_variable_eliminations.rend(); + ++it) { + const auto& elimination = *it; + f_t value = elimination.rhs; + for (size_t k = 0; k < elimination.columns.size(); ++k) { + value -= elimination.coefficients[k] * input_x[elimination.columns[k]]; + } + input_x[elimination.variable] = value / elimination.pivot_coefficient; + input_z[elimination.variable] = 0; + + if (elimination.pivot_row >= 0) { + f_t pivot_dual = 0; + for (size_t k = 0; k < elimination.affected_rows.size(); ++k) { + pivot_dual -= elimination.factors[k] * input_y[elimination.affected_rows[k]]; + } + input_y[elimination.pivot_row] = pivot_dual; + } + } + } + if (presolve_info.folding_info.is_folded) { // We solved a foled problem in the form // minimize c_prime^T x_prime diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 4c0eb31b3f..63d69aa42b 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -187,6 +187,21 @@ struct bounded_free_var_t { f_t coefficient; // a_{i*,j}: the coefficient of x_j in constraint i* }; +// Algebraic elimination of a free variable using an equality pivot row. +// All indices refer to the problem immediately before the free-variable elimination pass. +template +struct free_variable_elimination_t { + i_t variable; + i_t pivot_row; + f_t pivot_coefficient; + f_t rhs; + std::vector columns; + std::vector coefficients; + // Each affected row was replaced by row - factor * pivot_row. + std::vector affected_rows; + std::vector factors; +}; + template struct presolve_info_t { // indices of variables in the original problem that remain in the presolved problem @@ -218,6 +233,13 @@ struct presolve_info_t { // Originally-free variables that received implied bounds, with the constraint used std::vector> bounded_free_variables; + + // Free variables and pivot rows removed algebraically at the end of presolve. + std::vector> free_variable_eliminations; + std::vector free_elimination_remaining_variables; + std::vector free_elimination_remaining_constraints; + i_t free_elimination_num_variables{0}; + i_t free_elimination_num_constraints{0}; }; template diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 6f6ecc88f9..3b79db5abf 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -115,6 +115,28 @@ void write_matlab(const std::string& filename, const simplex::lp_problem_t +void compute_stationarity_residual(const lp_problem_t& lp, + const std::vector& x, + const std::vector& y, + const std::vector& z, + std::vector& residual) +{ + residual = z; + for (i_t j = 0; j < lp.num_cols; ++j) { + residual[j] -= lp.objective[j]; + } + if (lp.Q.n > 0) { + for (i_t i = 0; i < lp.Q.m; ++i) { + for (i_t p = lp.Q.row_start[i]; p < lp.Q.row_start[i + 1]; ++p) { + residual[i] -= lp.Q.x[p] * x[lp.Q.j[p]]; + } + } + } + matrix_transpose_vector_multiply(lp.A, 1.0, y, 1.0, residual); +} + } // namespace template @@ -569,19 +591,14 @@ lp_status_t solve_linear_program_with_barrier( settings.log.printf("Unscaled Primal infeasibility (abs/rel): %.2e/%.2e\n", primal_residual, primal_residual / (1.0 + vector_norm_inf(presolved_lp.rhs))); - if (barrier_lp.Q.n == 0) { - std::vector unscaled_dual_residual = unscaled_z; - for (i_t j = 0; j < unscaled_dual_residual.size(); ++j) { - unscaled_dual_residual[j] -= presolved_lp.objective[j]; - } - matrix_transpose_vector_multiply( - presolved_lp.A, 1.0, unscaled_y, 1.0, unscaled_dual_residual); - f_t unscaled_dual_residual_norm = vector_norm_inf(unscaled_dual_residual); - settings.log.printf( - "Unscaled Dual infeasibility (abs/rel): %.2e/%.2e\n", - unscaled_dual_residual_norm, - unscaled_dual_residual_norm / (1.0 + vector_norm_inf(presolved_lp.objective))); - } + std::vector unscaled_dual_residual; + compute_stationarity_residual( + presolved_lp, unscaled_x, unscaled_y, unscaled_z, unscaled_dual_residual); + f_t unscaled_dual_residual_norm = vector_norm_inf(unscaled_dual_residual); + settings.log.printf( + "Unscaled Dual infeasibility (abs/rel): %.2e/%.2e\n", + unscaled_dual_residual_norm, + unscaled_dual_residual_norm / (1.0 + vector_norm_inf(presolved_lp.objective))); } // Undo presolve @@ -604,19 +621,14 @@ lp_status_t solve_linear_program_with_barrier( post_solve_primal_residual, post_solve_primal_residual / (1.0 + vector_norm_inf(original_lp.rhs))); - if (barrier_lp.Q.n == 0) { - std::vector post_solve_dual_residual = lp_solution.z; - for (i_t j = 0; j < post_solve_dual_residual.size(); ++j) { - post_solve_dual_residual[j] -= original_lp.objective[j]; - } - matrix_transpose_vector_multiply( - original_lp.A, 1.0, lp_solution.y, 1.0, post_solve_dual_residual); - f_t post_solve_dual_residual_norm = vector_norm_inf(post_solve_dual_residual); - settings.log.printf( - "Post-solve Dual infeasibility (abs/rel): %.2e/%.2e\n", - post_solve_dual_residual_norm, - post_solve_dual_residual_norm / (1.0 + vector_norm_inf(original_lp.objective))); - } + std::vector post_solve_dual_residual; + compute_stationarity_residual( + original_lp, lp_solution.x, lp_solution.y, lp_solution.z, post_solve_dual_residual); + f_t post_solve_dual_residual_norm = vector_norm_inf(post_solve_dual_residual); + settings.log.printf( + "Post-solve Dual infeasibility (abs/rel): %.2e/%.2e\n", + post_solve_dual_residual_norm, + post_solve_dual_residual_norm / (1.0 + vector_norm_inf(original_lp.objective))); } if (dualize_info.solving_dual) { diff --git a/cpp/tests/socp/solve_barrier_socp.cu b/cpp/tests/socp/solve_barrier_socp.cu index 59fa339904..0009471974 100644 --- a/cpp/tests/socp/solve_barrier_socp.cu +++ b/cpp/tests/socp/solve_barrier_socp.cu @@ -13,6 +13,7 @@ #include #include #include +#include #include #include @@ -33,6 +34,85 @@ static void init_handler(const raft::handle_t* handle_ptr) handle_ptr->get_stream().get())); } +static double inf_norm(const std::vector& v) +{ + double nrm = 0.0; + for (double val : v) { + nrm = std::max(nrm, std::abs(val)); + } + return nrm; +} + +// Hub-and-spoke style chain: two free integrator variables plus a quadratic on w. +// +// minimize 0.5 w^2 +// s.t. y1 + w = 1 +// -y1 + y2 = 0 +// -y2 + w = 1 +// t - w = 0 +// (t, u) in Q^2 +// +// Unique primal: w = t = 1, y1 = y2 = u = 0. +static user_problem_t make_free_substitution_qp(raft::handle_t* handle) +{ + user_problem_t user_problem(handle); + + constexpr int m = 4; + constexpr int n = 5; + constexpr int nz = 8; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective.assign(n, 0.0); + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + // Columns: y1, y2, w, t, u + user_problem.A.col_start = {0, 2, 4, 7, 8, 8}; + user_problem.A.i = {0, 1, 1, 2, 0, 2, 3, 3}; + user_problem.A.x = {1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; + + user_problem.rhs = {1.0, 0.0, 1.0, 0.0}; + user_problem.row_sense = {'E', 'E', 'E', 'E'}; + // Keep y1, y2, and w free so bound strengthening cannot pin the integrator + // chain before substitution. w is skipped later because it appears in Q. + user_problem.lower = {-inf, -inf, -inf, 0.0, 0.0}; + user_problem.upper.assign(n, inf); + + user_problem.Q_offsets = {0, 0, 0, 1, 1, 1}; + user_problem.Q_indices = {2}; + user_problem.Q_values = {1.0}; + + user_problem.num_range_rows = 0; + user_problem.problem_name = "free_substitution_qp"; + user_problem.cone_var_start = 3; + user_problem.second_order_cone_dims = {2}; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + return user_problem; +} + +static void stationarity_residual(const lp_problem_t& lp, + const std::vector& x, + const std::vector& y, + const std::vector& z, + std::vector& residual) +{ + residual = z; + for (int j = 0; j < lp.num_cols; ++j) { + residual[j] -= lp.objective[j]; + } + if (lp.Q.n > 0) { + for (int i = 0; i < lp.Q.m; ++i) { + for (int p = lp.Q.row_start[i]; p < lp.Q.row_start[i + 1]; ++p) { + residual[i] -= lp.Q.x[p] * x[lp.Q.j[p]]; + } + } + } + matrix_transpose_vector_multiply(lp.A, 1.0, y, 1.0, residual); +} + TEST(barrier, cone_metadata_reindexed_when_slack_is_inserted_before_cones) { raft::handle_t handle{}; @@ -161,8 +241,9 @@ TEST(barrier, presolve_reindexes_cone_start_after_empty_column_removal) TEST(barrier, presolve_keeps_direct_free_variables_before_cones) { // Layout: [x0, x1 | cone x2, x3, x4] with x0, x1 free and a 3-dimensional SOC block. - // SOCP barrier presolve keeps direct free variables (no x = v - w split); cone_var_start - // and column count stay unchanged. + // Free linear columns are not split into v - w. Zero-cost free x0 is substituted from + // the singleton equality (the pivot may contain cone columns); the leftover free x1 + // then has an empty column and is fixed at 0. The cone block stays trailing. raft::handle_t handle{}; init_handler(&handle); @@ -210,17 +291,21 @@ TEST(barrier, presolve_keeps_direct_free_variables_before_cones) lp_problem_t presolved_lp(user_problem.handle_ptr, 1, 1, 1); ASSERT_EQ(presolve(original_lp, settings, presolved_lp, presolve_info), 0); - EXPECT_EQ(presolved_lp.num_cols, 5); - EXPECT_EQ(presolved_lp.cone_var_start, 2); + EXPECT_EQ(presolved_lp.num_rows, 0); + EXPECT_EQ(presolved_lp.num_cols, 3); + EXPECT_EQ(presolved_lp.cone_var_start, 0); EXPECT_EQ(presolved_lp.second_order_cone_dims, std::vector({3})); EXPECT_TRUE(presolve_info.free_variable_pairs.empty()); - ASSERT_EQ(presolve_info.direct_free_variables.size(), 2); - EXPECT_EQ(presolve_info.direct_free_variables[0], 0); - EXPECT_EQ(presolve_info.direct_free_variables[1], 1); - EXPECT_EQ(presolved_lp.lower[0], -inf); - EXPECT_EQ(presolved_lp.lower[1], -inf); - EXPECT_EQ(presolved_lp.upper[0], inf); - EXPECT_EQ(presolved_lp.upper[1], inf); + EXPECT_TRUE(presolve_info.direct_free_variables.empty()); + ASSERT_EQ(presolve_info.free_variable_eliminations.size(), 2u); + EXPECT_EQ(presolve_info.free_variable_eliminations[0].variable, 0); + EXPECT_EQ(presolve_info.free_variable_eliminations[0].pivot_row, 0); + EXPECT_EQ(presolve_info.free_variable_eliminations[1].variable, 1); + EXPECT_EQ(presolve_info.free_variable_eliminations[1].pivot_row, -1); + ASSERT_EQ(presolve_info.free_elimination_remaining_variables.size(), 3u); + EXPECT_EQ(presolve_info.free_elimination_remaining_variables[0], 2); + EXPECT_EQ(presolve_info.free_elimination_remaining_variables[1], 3); + EXPECT_EQ(presolve_info.free_elimination_remaining_variables[2], 4); } TEST(barrier, rejects_middle_cone_input_before_barrier) @@ -1020,4 +1105,95 @@ TEST(barrier, sparse_soc_expansion_solves_dim_500_cone) } } +TEST(barrier, free_variable_substitution_postsolve_kkt) +{ + raft::handle_t handle{}; + init_handler(&handle); + + user_problem_t user_problem = make_free_substitution_qp(&handle); + + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + settings.scale_columns = false; + settings.postsolve_info = 1; + + std::vector new_slacks; + dualize_info_t dualize_info; + lp_problem_t original_lp(user_problem.handle_ptr, 1, 1, 1); + convert_user_problem(user_problem, settings, original_lp, new_slacks, dualize_info); + + presolve_info_t presolve_info; + lp_problem_t presolved_lp(user_problem.handle_ptr, 1, 1, 1); + ASSERT_EQ(presolve(original_lp, settings, presolved_lp, presolve_info), 0); + ASSERT_EQ(presolve_info.free_variable_eliminations.size(), 2u); + ASSERT_EQ(presolved_lp.num_cols, original_lp.num_cols - 2); + ASSERT_EQ(presolved_lp.num_rows, original_lp.num_rows - 2); + + // Known feasible point on the original problem, restricted to remaining columns. + const std::vector original_x = {0.0, 0.0, 1.0, 1.0, 0.0}; + std::vector crushed_x(presolved_lp.num_cols); + for (int k = 0; k < presolved_lp.num_cols; ++k) { + crushed_x[k] = original_x[presolve_info.free_elimination_remaining_variables[k]]; + } + std::vector crushed_y(presolved_lp.num_rows, 0.25); + std::vector crushed_z(presolved_lp.num_cols, 0.0); + std::vector reduced_stationarity; + stationarity_residual(presolved_lp, crushed_x, crushed_y, crushed_z, reduced_stationarity); + for (int j = 0; j < presolved_lp.num_cols; ++j) { + crushed_z[j] -= reduced_stationarity[j]; + } + stationarity_residual(presolved_lp, crushed_x, crushed_y, crushed_z, reduced_stationarity); + ASSERT_NEAR(inf_norm(reduced_stationarity), 0.0, 1e-12); + + std::vector uncrushed_x; + std::vector uncrushed_y; + std::vector uncrushed_z; + uncrush_solution(presolve_info, + settings, + original_lp, + crushed_x, + crushed_y, + crushed_z, + uncrushed_x, + uncrushed_y, + uncrushed_z); + + ASSERT_EQ(uncrushed_x.size(), static_cast(original_lp.num_cols)); + EXPECT_NEAR(uncrushed_x[0], 0.0, 1e-12); + EXPECT_NEAR(uncrushed_x[1], 0.0, 1e-12); + EXPECT_NEAR(uncrushed_x[2], 1.0, 1e-12); + EXPECT_NEAR(uncrushed_x[3], 1.0, 1e-12); + EXPECT_NEAR(std::abs(uncrushed_z[0]), 0.0, 1e-12); + EXPECT_NEAR(std::abs(uncrushed_z[1]), 0.0, 1e-12); + + std::vector primal_residual = original_lp.rhs; + matrix_vector_multiply(original_lp.A, 1.0, uncrushed_x, -1.0, primal_residual); + EXPECT_NEAR(inf_norm(primal_residual), 0.0, 1e-12); + + std::vector dual_residual; + stationarity_residual(original_lp, uncrushed_x, uncrushed_y, uncrushed_z, dual_residual); + EXPECT_NEAR(inf_norm(dual_residual), 0.0, 1e-10); + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + auto status = solve_linear_program_with_barrier(user_problem, settings, solution); + EXPECT_EQ(status, lp_status_t::OPTIMAL); + EXPECT_NEAR(solution.objective, 0.5, 1e-4); + EXPECT_NEAR(solution.x[0], 0.0, 1e-4); + EXPECT_NEAR(solution.x[1], 0.0, 1e-4); + EXPECT_NEAR(solution.x[2], 1.0, 1e-4); + EXPECT_NEAR(solution.x[3], 1.0, 1e-4); + EXPECT_NEAR(std::abs(solution.z[0]), 0.0, 1e-4); + EXPECT_NEAR(std::abs(solution.z[1]), 0.0, 1e-4); + + std::vector solved_primal = original_lp.rhs; + matrix_vector_multiply(original_lp.A, 1.0, solution.x, -1.0, solved_primal); + EXPECT_NEAR(inf_norm(solved_primal), 0.0, 1e-5); + + std::vector solved_dual; + stationarity_residual(original_lp, solution.x, solution.y, solution.z, solved_dual); + EXPECT_NEAR(inf_norm(solved_dual), 0.0, 1e-5); +} + } // namespace cuopt::mathematical_optimization::simplex::test From cd9e70c7e39842cfee7ce0e16b2dcc5bad581ed6 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 16 Sep 2026 11:51:14 -0700 Subject: [PATCH 08/23] Do not sort row entries before the CSR to CSC conversion --- cpp/src/dual_simplex/presolve.cpp | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 7a9f8a227b..f6aaa5dd8a 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -282,14 +282,11 @@ static void eliminate_free_variables(lp_problem_t& problem, for (i_t new_i = 0; new_i < new_m; ++new_i) { const i_t old_i = remaining_rows[new_i]; reduced_A.row_start[new_i] = nz; - std::vector> entries; - entries.reserve(rows[old_i].size()); + // Column order within a row is irrelevant here: to_compressed_col below buckets the + // entries by column in linear time, so sorting each row would only add O(nnz log nnz). for (const auto& [old_j, value] : rows[old_i]) { - if (active_col[old_j] && value != 0) { entries.emplace_back(old_to_new_col[old_j], value); } - } - std::sort(entries.begin(), entries.end()); - for (const auto& [new_j, value] : entries) { - reduced_A.j[nz] = new_j; + if (!active_col[old_j] || value == 0) { continue; } + reduced_A.j[nz] = old_to_new_col[old_j]; reduced_A.x[nz] = value; ++nz; } From 168a1b447be4f1edb4b2f27e82c010cba0620ea0 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 16 Sep 2026 15:03:46 -0700 Subject: [PATCH 09/23] Replace per-row hash maps with flat CSR/CSC arenas for free-variable substitution --- cpp/src/dual_simplex/presolve.cpp | 378 +++++++++++++++++++++++++----- 1 file changed, 317 insertions(+), 61 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index f6aaa5dd8a..038a625ba8 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -16,7 +16,6 @@ #include #include #include -#include namespace cuopt::mathematical_optimization::simplex { @@ -96,6 +95,256 @@ static void remove_variables_from_Q(csr_matrix_t& Q, Q = std::move(Qout); } +// Sparse matrix used while substituting free variables. +// +// Rows and columns live in two flat arenas (row-major and column-major, each with a little +// slack per line) rather than a hash map per row and a vector per column. Random access +// into a row goes through scatter, which indexes one row at a time; substitution keeps +// updating the same target row, so that row stays resident and each update costs +// O(pivot nonzeros). Re-scattering per elimination is quadratic when one row is dense. +template +struct substitution_matrix_t { + // Rows at or below this length are searched directly, so glancing at a short row does not + // evict a long resident one. + static constexpr i_t scan_limit = 16; + + i_t num_rows = 0; + i_t num_cols = 0; + std::vector row_col; + std::vector row_val; + std::vector row_start; + std::vector row_len; + std::vector row_cap; + // Row indices per column, append only: an entry whose row no longer holds the column is + // dropped when that column is scanned. + std::vector col_row; + std::vector col_start; + std::vector col_len; + std::vector col_cap; + std::vector scatter; // column -> offset within the resident row, -1 when absent + i_t scattered = -1; + // Arena high-water marks. A line that outgrows its slot is relocated to the tail with + // twice the capacity, so appends are amortized O(1); the holes left behind are reclaimed + // by a compaction once they outweigh the live entries. + i_t row_used = 0; + i_t col_used = 0; + + static i_t slack(i_t len) { return std::max(4, len / 8); } + + void build(const lp_problem_t& problem) + { + num_rows = problem.num_rows; + num_cols = problem.num_cols; + const i_t nnz = problem.A.col_start[num_cols]; + + // Counting pass over the CSC input gives the row lengths, then a prefix sum lays out + // the row arena and a single scatter fills it. + row_len.assign(num_rows, 0); + for (i_t p = 0; p < nnz; ++p) { + ++row_len[problem.A.i[p]]; + } + row_start.resize(num_rows); + row_cap.resize(num_rows); + i_t used = 0; + for (i_t i = 0; i < num_rows; ++i) { + row_start[i] = used; + row_cap[i] = row_len[i] + slack(row_len[i]); + used += row_cap[i]; + } + row_used = used; + row_col.assign(used, 0); + row_val.assign(used, 0); + std::vector cursor(row_start); + for (i_t j = 0; j < num_cols; ++j) { + for (i_t p = problem.A.col_start[j]; p < problem.A.col_start[j + 1]; ++p) { + const i_t q = cursor[problem.A.i[p]]++; + row_col[q] = j; + row_val[q] = problem.A.x[p]; + } + } + + col_start.resize(num_cols); + col_len.resize(num_cols); + col_cap.resize(num_cols); + used = 0; + for (i_t j = 0; j < num_cols; ++j) { + col_len[j] = problem.A.col_start[j + 1] - problem.A.col_start[j]; + col_start[j] = used; + col_cap[j] = col_len[j] + slack(col_len[j]); + used += col_cap[j]; + } + col_used = used; + col_row.assign(used, 0); + for (i_t j = 0; j < num_cols; ++j) { + i_t q = col_start[j]; + for (i_t p = problem.A.col_start[j]; p < problem.A.col_start[j + 1]; ++p, ++q) { + col_row[q] = problem.A.i[p]; + } + } + + scatter.assign(num_cols, -1); + } + + i_t column(i_t i, i_t offset) const { return row_col[row_start[i] + offset]; } + f_t value(i_t i, i_t offset) const { return row_val[row_start[i] + offset]; } + void set_value(i_t i, i_t offset, f_t value) { row_val[row_start[i] + offset] = value; } + + void unload() + { + if (scattered == -1) { return; } + const i_t start = row_start[scattered]; + for (i_t k = 0; k < row_len[scattered]; ++k) { + scatter[row_col[start + k]] = -1; + } + scattered = -1; + } + + void load(i_t i) + { + if (scattered == i) { return; } + unload(); + const i_t start = row_start[i]; + for (i_t k = 0; k < row_len[i]; ++k) { + scatter[row_col[start + k]] = k; + } + scattered = i; + } + + i_t scan(i_t i, i_t col) const + { + const i_t start = row_start[i]; + for (i_t k = 0; k < row_len[i]; ++k) { + if (row_col[start + k] == col) { return k; } + } + return -1; + } + + // Offset of (i, col) within row i, or -1 when absent. + i_t find(i_t i, i_t col) + { + if (scattered == i) { return scatter[col]; } + if (row_len[i] <= scan_limit) { return scan(i, col); } + load(i); + return scatter[col]; + } + + // erase and insert both require row i to be resident. + void erase(i_t i, i_t offset) + { + const i_t start = row_start[i]; + const i_t last = row_len[i] - 1; + scatter[row_col[start + offset]] = -1; + if (offset != last) { + row_col[start + offset] = row_col[start + last]; + row_val[start + offset] = row_val[start + last]; + scatter[row_col[start + offset]] = offset; + } + row_len[i] = last; + } + + void insert(i_t i, i_t col, f_t value) + { + reserve_row(i); + const i_t offset = row_len[i]; + row_col[row_start[i] + offset] = col; + row_val[row_start[i] + offset] = value; + scatter[col] = offset; + row_len[i] = offset + 1; + reserve_col(col); + col_row[col_start[col] + col_len[col]] = i; + ++col_len[col]; + } + + // Relocation and compaction both preserve the order within a line, so resident offsets + // stay valid across either. + void reserve_row(i_t i) + { + if (row_len[i] < row_cap[i]) { return; } + const i_t need = std::max(row_len[i] + 1, 2 * row_cap[i]); + if (row_used + need > static_cast(row_col.size())) { + i_t live = 0; + for (i_t r = 0; r < num_rows; ++r) { + live += row_len[r]; + } + if (row_used > 2 * (live + num_rows)) { compact_rows(); } + if (row_used + need > static_cast(row_col.size())) { + const i_t size = + std::max(row_used + need, 2 * static_cast(row_col.size()) + 1); + row_col.resize(size, 0); + row_val.resize(size, 0); + } + if (row_len[i] < row_cap[i]) { return; } + } + std::copy_n(row_col.begin() + row_start[i], row_len[i], row_col.begin() + row_used); + std::copy_n(row_val.begin() + row_start[i], row_len[i], row_val.begin() + row_used); + row_start[i] = row_used; + row_cap[i] = need; + row_used += need; + } + + void reserve_col(i_t j) + { + if (col_len[j] < col_cap[j]) { return; } + const i_t need = std::max(col_len[j] + 1, 2 * col_cap[j]); + if (col_used + need > static_cast(col_row.size())) { + i_t live = 0; + for (i_t c = 0; c < num_cols; ++c) { + live += col_len[c]; + } + if (col_used > 2 * (live + num_cols)) { compact_cols(); } + if (col_used + need > static_cast(col_row.size())) { + const i_t size = + std::max(col_used + need, 2 * static_cast(col_row.size()) + 1); + col_row.resize(size, 0); + } + if (col_len[j] < col_cap[j]) { return; } + } + std::copy_n(col_row.begin() + col_start[j], col_len[j], col_row.begin() + col_used); + col_start[j] = col_used; + col_cap[j] = need; + col_used += need; + } + + void compact_rows() + { + std::vector new_start(num_rows); + i_t used = 0; + for (i_t i = 0; i < num_rows; ++i) { + new_start[i] = used; + used += row_len[i] + slack(row_len[i]); + } + std::vector new_col(std::max(used, 1), 0); + std::vector new_val(std::max(used, 1), 0); + for (i_t i = 0; i < num_rows; ++i) { + std::copy_n(row_col.begin() + row_start[i], row_len[i], new_col.begin() + new_start[i]); + std::copy_n(row_val.begin() + row_start[i], row_len[i], new_val.begin() + new_start[i]); + row_cap[i] = row_len[i] + slack(row_len[i]); + } + row_col = std::move(new_col); + row_val = std::move(new_val); + row_start = std::move(new_start); + row_used = used; + } + + void compact_cols() + { + std::vector new_start(num_cols); + i_t used = 0; + for (i_t j = 0; j < num_cols; ++j) { + new_start[j] = used; + used += col_len[j] + slack(col_len[j]); + } + std::vector new_row(std::max(used, 1), 0); + for (i_t j = 0; j < num_cols; ++j) { + std::copy_n(col_row.begin() + col_start[j], col_len[j], new_row.begin() + new_start[j]); + col_cap[j] = col_len[j] + slack(col_len[j]); + } + col_row = std::move(new_row); + col_start = std::move(new_start); + col_used = used; + } +}; + // Eliminate zero-cost, Q-uncoupled free linear variables by sparse equality substitution. // This is especially useful for conic formulations containing chains of auxiliary free // variables: unlike regularizing those variables in the KKT system, substitution is exact. @@ -110,37 +359,26 @@ static void eliminate_free_variables(lp_problem_t& problem, std::vector q_present(old_n, 0); if (problem.Q.n > 0) { - for (i_t row = 0; row < problem.Q.m; ++row) { - for (i_t p = problem.Q.row_start[row]; p < problem.Q.row_start[row + 1]; ++p) { - if (problem.Q.x[p] == 0) { continue; } - q_present[row] = 1; - q_present[problem.Q.j[p]] = 1; - } + // Q is square and symmetric, so a nonempty row is exactly the Q-coupled variables. + const i_t q_n = std::min(problem.Q.m, old_n); + for (i_t row = 0; row < q_n; ++row) { + if (problem.Q.row_start[row + 1] > problem.Q.row_start[row]) { q_present[row] = 1; } } } - csr_matrix_t Arow(0, 0, 0); - problem.A.to_compressed_row(Arow); - std::vector> rows(old_m); - std::vector> col_rows(old_n); - for (i_t i = 0; i < old_m; ++i) { - auto& row = rows[i]; - row.reserve(static_cast(Arow.row_start[i + 1] - Arow.row_start[i]) * 2); - for (i_t p = Arow.row_start[i]; p < Arow.row_start[i + 1]; ++p) { - row[Arow.j[p]] = Arow.x[p]; - col_rows[Arow.j[p]].push_back(i); - } - } + substitution_matrix_t matrix; + matrix.build(problem); std::vector active_row(old_m, 1); std::vector active_col(old_n, 1); auto& eliminations = presolve_info.free_variable_eliminations; - // col_rows is append-only: fill re-adds a row that may already be listed, and eliminated - // rows are never unlisted. Compact it during the scan so the pass stays linear. + // A column list is append-only: fill re-adds a row that may already be listed, and + // eliminated rows are never unlisted. Compact it during the scan so the pass stays linear. std::vector row_stamp(old_m, -1); i_t stamp = 0; std::vector incident; + std::vector incident_value; // One pass in column order. Peeling a chain from one end keeps each pivot row sparse, so // revisiting columns buys almost nothing and costs a requeue storm on models with @@ -152,18 +390,24 @@ static void eliminate_free_variables(lp_problem_t& problem, } incident.clear(); - auto& listed = col_rows[j]; - size_t keep = 0; + incident_value.clear(); + const i_t listed = matrix.col_start[j]; + i_t keep = 0; ++stamp; - for (size_t idx = 0; idx < listed.size(); ++idx) { - const i_t i = listed[idx]; + for (i_t idx = 0; idx < matrix.col_len[j]; ++idx) { + const i_t i = matrix.col_row[listed + idx]; if (row_stamp[i] == stamp) { continue; } row_stamp[i] = stamp; - if (!active_row[i] || rows[i].find(j) == rows[i].end()) { continue; } - listed[keep++] = i; + if (!active_row[i]) { continue; } + // Keep a_ij from this probe; pivot selection and the update reuse it. + const i_t offset = matrix.find(i, j); + if (offset == -1 || matrix.value(i, offset) == 0) { continue; } + matrix.col_row[listed + keep] = i; + ++keep; incident.push_back(i); + incident_value.push_back(matrix.value(i, offset)); } - listed.resize(keep); + matrix.col_len[j] = keep; free_variable_elimination_t elimination; elimination.variable = j; @@ -177,17 +421,32 @@ static void eliminate_free_variables(lp_problem_t& problem, } // A small pivot row limits fill. On a path this peels sparse boundary rows instead of - // repeatedly traversing the growing aggregate row. - i_t pivot = incident.front(); - for (const i_t i : incident) { - if (rows[i].size() < rows[pivot].size() || - (rows[i].size() == rows[pivot].size() && - std::abs(rows[i].at(j)) > std::abs(rows[pivot].at(j)))) { - pivot = i; + // repeatedly traversing the growing aggregate row. Row lengths come from the arena, so + // the dense row loses on length alone and is never scanned here. + size_t pivot_slot = 0; + for (size_t slot = 1; slot < incident.size(); ++slot) { + const i_t len = matrix.row_len[incident[slot]]; + const i_t best = matrix.row_len[incident[pivot_slot]]; + if (len < best || + (len == best && + std::abs(incident_value[slot]) > std::abs(incident_value[pivot_slot]))) { + pivot_slot = slot; } } - const f_t pivot_coefficient = rows[pivot].at(j); - if (pivot_coefficient == 0) { continue; } + const i_t pivot = incident[pivot_slot]; + const f_t pivot_coefficient = incident_value[pivot_slot]; + const i_t pivot_len = matrix.row_len[pivot]; + + // Snapshot the pivot row before touching the matrix: the arena can relocate rows when a + // substitution inserts, and postsolve needs these coefficients anyway. + elimination.columns.reserve(pivot_len - 1); + elimination.coefficients.reserve(pivot_len - 1); + for (i_t k = 0; k < pivot_len; ++k) { + const i_t col = matrix.column(pivot, k); + if (col == j) { continue; } + elimination.columns.push_back(col); + elimination.coefficients.push_back(matrix.value(pivot, k)); + } // Accept only if the substitution does not add nonzeros. Dropping the pivot row and the // a_ij entries pays for the entries the pivot row scatters into the other incident rows. @@ -197,45 +456,38 @@ static void eliminate_free_variables(lp_problem_t& problem, i_t added = 0; for (const i_t i : incident) { if (i == pivot) { continue; } - for (const auto& [col, value] : rows[pivot]) { - if (col == j) { continue; } - if (rows[i].find(col) == rows[i].end()) { ++added; } + for (const i_t col : elimination.columns) { + if (matrix.find(i, col) == -1) { ++added; } } } - const i_t removed = - static_cast(rows[pivot].size()) + static_cast(incident.size()) - 1; + const i_t removed = pivot_len + static_cast(incident.size()) - 1; if (added > removed) { continue; } elimination.pivot_row = pivot; elimination.pivot_coefficient = pivot_coefficient; elimination.rhs = problem.rhs[pivot]; - elimination.columns.reserve(rows[pivot].size() - 1); - elimination.coefficients.reserve(rows[pivot].size() - 1); - for (const auto& [col, value] : rows[pivot]) { - if (col == j) { continue; } - elimination.columns.push_back(col); - elimination.coefficients.push_back(value); - } for (const i_t i : incident) { if (i == pivot) { continue; } - auto j_it = rows[i].find(j); - if (j_it == rows[i].end()) { continue; } - const f_t factor = j_it->second / pivot_coefficient; - rows[i].erase(j_it); + matrix.load(i); + const i_t j_offset = matrix.scatter[j]; + if (j_offset == -1) { continue; } + const f_t factor = matrix.value(i, j_offset) / pivot_coefficient; + matrix.erase(i, j_offset); for (size_t k = 0; k < elimination.columns.size(); ++k) { const i_t col = elimination.columns[k]; const f_t delta = factor * elimination.coefficients[k]; - auto col_it = rows[i].find(col); - const f_t old_value = col_it == rows[i].end() ? f_t{0} : col_it->second; + const i_t offset = matrix.scatter[col]; + const f_t old_value = offset == -1 ? f_t{0} : matrix.value(i, offset); const f_t new_value = old_value - delta; const f_t drop_tol = f_t{100} * std::numeric_limits::epsilon() * std::max({f_t{1}, std::abs(old_value), std::abs(delta)}); if (std::abs(new_value) <= drop_tol) { - if (col_it != rows[i].end()) { rows[i].erase(col_it); } + if (offset != -1) { matrix.erase(i, offset); } + } else if (offset == -1) { + matrix.insert(i, col, new_value); } else { - if (col_it == rows[i].end()) { col_rows[col].push_back(i); } - rows[i][col] = new_value; + matrix.set_value(i, offset, new_value); } } problem.rhs[i] -= factor * elimination.rhs; @@ -243,6 +495,7 @@ static void eliminate_free_variables(lp_problem_t& problem, elimination.factors.push_back(factor); } + if (matrix.scattered == pivot) { matrix.unload(); } active_row[pivot] = 0; active_col[j] = 0; eliminations.push_back(std::move(elimination)); @@ -269,10 +522,11 @@ static void eliminate_free_variables(lp_problem_t& problem, i_t new_n = static_cast(remaining_cols.size()); i_t new_m = static_cast(remaining_rows.size()); + matrix.unload(); i_t new_nnz = 0; for (const i_t i : remaining_rows) { - for (const auto& [j, value] : rows[i]) { - if (active_col[j] && value != 0) { ++new_nnz; } + for (i_t k = 0; k < matrix.row_len[i]; ++k) { + if (active_col[matrix.column(i, k)] && matrix.value(i, k) != 0) { ++new_nnz; } } } @@ -284,7 +538,9 @@ static void eliminate_free_variables(lp_problem_t& problem, reduced_A.row_start[new_i] = nz; // Column order within a row is irrelevant here: to_compressed_col below buckets the // entries by column in linear time, so sorting each row would only add O(nnz log nnz). - for (const auto& [old_j, value] : rows[old_i]) { + for (i_t k = 0; k < matrix.row_len[old_i]; ++k) { + const i_t old_j = matrix.column(old_i, k); + const f_t value = matrix.value(old_i, k); if (!active_col[old_j] || value == 0) { continue; } reduced_A.j[nz] = old_to_new_col[old_j]; reduced_A.x[nz] = value; From 6b0f33230e928c5a3e1a6c81531677720e8d00cc Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Tue, 22 Sep 2026 11:20:00 -0700 Subject: [PATCH 10/23] Use threshold pivoting when eliminating free variables by substitution. A pivot is accepted only if it is large relative to both its row and its column, otherwise the variable stays free in the augmented system. --- .../mathematical_optimization/constants.h | 5 ++ .../pdlp/solver_settings.hpp | 7 ++ cpp/src/dual_simplex/presolve.cpp | 85 ++++++++++++++----- .../dual_simplex/simplex_solver_settings.hpp | 9 ++ cpp/src/math_optimization/solver_settings.cu | 3 + cpp/src/pdlp/solve.cu | 5 ++ cpp/tests/socp/solve_barrier_socp.cu | 75 ++++++++++++++++ 7 files changed, 169 insertions(+), 20 deletions(-) diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 3656791a98..ae1d6153c1 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -161,6 +161,11 @@ /* @brief Barrier initial point safeguard */ #define CUOPT_BARRIER_INITIAL_POINT_SAFEGUARD "barrier_initial_point_safeguard" +/* @brief Free-variable equality substitution in barrier presolve */ +#define CUOPT_BARRIER_PRESOLVE_FREE_ELIMINATION "barrier_presolve_free_elimination" +#define CUOPT_BARRIER_FREE_ELIMINATION_ROW_PIVOT_TOL "barrier_free_elimination_row_pivot_tol" +#define CUOPT_BARRIER_FREE_ELIMINATION_COL_PIVOT_TOL "barrier_free_elimination_col_pivot_tol" + /* @brief MIP determinism mode constants */ #define CUOPT_MODE_OPPORTUNISTIC 0 #define CUOPT_MODE_DETERMINISTIC 1 diff --git a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp index fd6c497850..bc1040c3f3 100644 --- a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp @@ -303,6 +303,13 @@ class pdlp_solver_settings_t { barrier_dual_initial_point_t barrier_dual_initial_point{barrier_dual_initial_point_t::Automatic}; i_t postsolve_info{-1}; i_t barrier_presolve_bound_free_variables{-1}; // -1 automatic, 0 disabled, 1 enabled + i_t barrier_presolve_free_elimination{-1}; // -1 automatic, 0 disabled, 1 enabled + // Threshold pivoting tolerances for the free-variable equality substitution. The pivot + // a_pj must be at least this fraction of the largest other entry in its row and of the + // largest entry in its column; 0 disables the corresponding test. At 1.0 no substitution + // can amplify a coefficient. + f_t barrier_free_elimination_row_pivot_tol{1.0}; + f_t barrier_free_elimination_col_pivot_tol{1.0}; // Ruiz equilibration for QCQP (barrier) scaling: -1 automatic (row/column // imbalance heuristic), 0 disabled, 1 enabled. Distinct from PDLP's own Ruiz // scaling in pdlp_hyper_params_t. diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index cbbeb18229..b96e7492de 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -17,6 +17,7 @@ #include #include #include +#include namespace cuopt::mathematical_optimization::simplex { @@ -349,14 +350,16 @@ struct substitution_matrix_t { // Eliminate zero-cost, Q-uncoupled free linear variables by sparse equality substitution. // This is especially useful for conic formulations containing chains of auxiliary free // variables: unlike regularizing those variables in the KKT system, substitution is exact. +// Returns the number of columns left alone because no incident row passed the pivot threshold. template -static void eliminate_free_variables(lp_problem_t& problem, - presolve_info_t& presolve_info) +static i_t eliminate_free_variables(lp_problem_t& problem, + presolve_info_t& presolve_info, + const simplex_solver_settings_t& settings) { const i_t old_m = problem.num_rows; const i_t old_n = problem.num_cols; const i_t linear_cols = linear_variable_count(problem); - if (old_m == 0 || linear_cols == 0) { return; } + if (old_m == 0 || linear_cols == 0) { return 0; } std::vector q_present(old_n, 0); if (problem.Q.n > 0) { @@ -380,6 +383,11 @@ static void eliminate_free_variables(lp_problem_t& problem, i_t stamp = 0; std::vector incident; std::vector incident_value; + std::vector candidate_order; + + const f_t row_pivot_tol = settings.barrier_free_elimination_row_pivot_tol; + const f_t col_pivot_tol = settings.barrier_free_elimination_col_pivot_tol; + i_t pivot_rejected = 0; // One pass in column order. Peeling a chain from one end keeps each pivot row sparse, so // revisiting columns buys almost nothing and costs a requeue storm on models with @@ -422,17 +430,47 @@ static void eliminate_free_variables(lp_problem_t& problem, } // A small pivot row limits fill. On a path this peels sparse boundary rows instead of - // repeatedly traversing the growing aggregate row. Row lengths come from the arena, so - // the dense row loses on length alone and is never scanned here. - size_t pivot_slot = 0; - for (size_t slot = 1; slot < incident.size(); ++slot) { - const i_t len = matrix.row_len[incident[slot]]; - const i_t best = matrix.row_len[incident[pivot_slot]]; - if (len < best || - (len == best && - std::abs(incident_value[slot]) > std::abs(incident_value[pivot_slot]))) { - pivot_slot = slot; + // repeatedly traversing the growing aggregate row. Candidates are ranked shortest first + // and the first one passing the threshold tests wins, so a long aggregate row is only + // scanned once every shorter row has failed. + candidate_order.resize(incident.size()); + std::iota(candidate_order.begin(), candidate_order.end(), 0); + std::sort(candidate_order.begin(), candidate_order.end(), [&](i_t a, i_t b) { + const i_t len_a = matrix.row_len[incident[a]]; + const i_t len_b = matrix.row_len[incident[b]]; + if (len_a != len_b) { return len_a < len_b; } + return std::abs(incident_value[a]) > std::abs(incident_value[b]); + }); + + // Threshold pivoting, as in the LU: a pivot that is tiny relative to its own row or to the + // rest of its column makes the multiplier a_ij / a_pj huge and inflates every surviving row + // by that factor, which the barrier's KKT solve cannot recover from. The column test bounds + // the multiplier, the row test bounds the coefficients the pivot row scatters. At tolerance + // 1 the pivot is the largest entry of both its row and its column, so no coefficient grows. + f_t column_max = 0; + for (const f_t value : incident_value) { + column_max = std::max(column_max, std::abs(value)); + } + i_t pivot_slot = -1; + for (const i_t slot : candidate_order) { + const f_t a_ij = std::abs(incident_value[slot]); + if (a_ij < col_pivot_tol * column_max) { continue; } + if (row_pivot_tol > 0) { + const i_t row = incident[slot]; + f_t row_max = 0; + for (i_t k = 0; k < matrix.row_len[row]; ++k) { + if (matrix.column(row, k) == j) { continue; } + row_max = std::max(row_max, std::abs(matrix.value(row, k))); + } + if (a_ij < row_pivot_tol * row_max) { continue; } } + pivot_slot = slot; + break; + } + // No stable pivot: leave the column as a free variable handled directly in the KKT system. + if (pivot_slot == -1) { + ++pivot_rejected; + continue; } const i_t pivot = incident[pivot_slot]; const f_t pivot_coefficient = incident_value[pivot_slot]; @@ -502,7 +540,7 @@ static void eliminate_free_variables(lp_problem_t& problem, eliminations.push_back(std::move(elimination)); } - if (eliminations.empty()) { return; } + if (eliminations.empty()) { return pivot_rejected; } presolve_info.free_elimination_num_variables = old_n; presolve_info.free_elimination_num_constraints = old_m; @@ -582,6 +620,7 @@ static void eliminate_free_variables(lp_problem_t& problem, presolve_info.direct_free_variables.push_back(new_j); } } + return pivot_rejected; } template @@ -1990,16 +2029,22 @@ i_t presolve(const lp_problem_t& original, } // LP already goes through PSLP; this substitution is for QP/SOCP only. - if (settings.barrier_presolve && (has_cones || problem.Q.n > 0)) { + if (settings.barrier_presolve && settings.barrier_presolve_free_elimination != 0 && + (has_cones || problem.Q.n > 0)) { const i_t old_free_count = static_cast(presolve_info.direct_free_variables.size()); const f_t free_elimination_start = tic(); - eliminate_free_variables(problem, presolve_info); + const i_t pivot_rejected = eliminate_free_variables(problem, presolve_info, settings); const i_t eliminated = old_free_count - static_cast(presolve_info.direct_free_variables.size()); - if (eliminated > 0) { - settings.log.printf("Eliminated %d free variables by equality substitution in %.2fs\n", - eliminated, - toc(free_elimination_start)); + if (eliminated > 0 || pivot_rejected > 0) { + settings.log.printf( + "Eliminated %d free variables by equality substitution in %.2fs (%d skipped by pivot " + "threshold, row tol %g, col tol %g)\n", + eliminated, + toc(free_elimination_start), + pivot_rejected, + settings.barrier_free_elimination_row_pivot_tol, + settings.barrier_free_elimination_col_pivot_tol); } } diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 8b3eba56d3..45e5b28fa9 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -84,6 +84,9 @@ struct simplex_solver_settings_t { barrier_dual_initial_point(barrier_dual_initial_point_t::Automatic), postsolve_info(-1), barrier_presolve_bound_free_variables(-1), + barrier_presolve_free_elimination(-1), + barrier_free_elimination_row_pivot_tol(1.0), + barrier_free_elimination_col_pivot_tol(1.0), qcqp_ruiz_equilibration(-1), barrier_initial_point_safeguard(10.0), check_Q(false), @@ -193,6 +196,12 @@ struct simplex_solver_settings_t { // 1 dual least squares, 2 SeDuMi mu-based i_t postsolve_info; // -1 automatic (disabled), 0 disabled, 1 enabled i_t barrier_presolve_bound_free_variables; // -1 automatic, 0 disabled, 1 enabled + i_t barrier_presolve_free_elimination; // -1 automatic, 0 disabled, 1 enabled + // Threshold pivoting tolerances for the free-variable equality substitution: the pivot must + // be at least this fraction of the largest other entry of its row / of the largest entry of + // its column. 0 disables a test, 1 forbids any coefficient growth. + f_t barrier_free_elimination_row_pivot_tol; + f_t barrier_free_elimination_col_pivot_tol; i_t qcqp_ruiz_equilibration; // -1 automatic (imbalance heuristic), 0 disabled, 1 enabled f_t barrier_initial_point_safeguard; // margin pushing the barrier initial iterate into // the interior of the nonnegative orthant / SOC diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 5a4ab72c32..2e92514bc1 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -145,6 +145,8 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_BARRIER_DUAL_REGULARIZATION, &pdlp_settings.barrier_dual_regularization, f_t(-1.0), std::numeric_limits::infinity(), f_t(-1.0), "initial dual regularization for the augmented system; -1 automatic"}, {CUOPT_BARRIER_STEP_SCALE, &pdlp_settings.barrier_step_scale, f_t(0.5), f_t(0.9999), f_t(0.9)}, {CUOPT_BARRIER_INITIAL_POINT_SAFEGUARD, &pdlp_settings.barrier_initial_point_safeguard, f_t(0.0), std::numeric_limits::infinity(), f_t(10.0), "margin pushing the barrier initial iterate into the interior of the nonnegative orthant / SOC"}, + {CUOPT_BARRIER_FREE_ELIMINATION_ROW_PIVOT_TOL, &pdlp_settings.barrier_free_elimination_row_pivot_tol, f_t(0.0), f_t(1.0), f_t(1.0), "free-variable elimination: pivot must be this fraction of the largest other entry in its row; 0 disables the test"}, + {CUOPT_BARRIER_FREE_ELIMINATION_COL_PIVOT_TOL, &pdlp_settings.barrier_free_elimination_col_pivot_tol, f_t(0.0), f_t(1.0), f_t(1.0), "free-variable elimination: pivot must be this fraction of the largest entry in its column; 0 disables the test"}, // MIP heuristic hyper-parameters (hidden from default --help: name contains "hyper_") {CUOPT_MIP_HYPER_HEURISTIC_ROOT_LP_TIME_RATIO, &mip_settings.heuristic_params.root_lp_time_ratio, f_t(0.0), f_t(1.0), f_t(0.1), "fraction of total time for root LP"}, {CUOPT_MIP_HYPER_HEURISTIC_ROOT_LP_MAX_TIME, &mip_settings.heuristic_params.root_lp_max_time, f_t(0.0), std::numeric_limits::infinity(), f_t(15.0), "hard cap on root LP seconds"}, @@ -232,6 +234,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_HYPER_SUBMIP_ITERATION_LIMIT_OFFSET, &mip_settings.submip_params.iteration_limit_offset, 0, std::numeric_limits::max(), 10000, "base sub-MIP simplex-iteration limit for root heuristics"}, {CUOPT_MIP_HYPER_SUBMIP_MAX_LEVEL, &mip_settings.submip_params.max_level, 0, std::numeric_limits::max(), 10, "maximum sub-MIP recursion level"}, {CUOPT_BARRIER_PRESOLVE_BOUND_FREE_VARIABLES, &pdlp_settings.barrier_presolve_bound_free_variables, -1, 1, -1, "Bound free variables during barrier presolve: -1 automatic (default behavior), 0 disabled, 1 enabled"}, + {CUOPT_BARRIER_PRESOLVE_FREE_ELIMINATION, &pdlp_settings.barrier_presolve_free_elimination, -1, 1, -1, "Eliminate zero-cost free variables by equality substitution during barrier presolve: -1 automatic (default behavior), 0 disabled, 1 enabled"}, {CUOPT_BARRIER_ADAPTIVE_REGULARIZATION, &pdlp_settings.barrier_adaptive_regularization, -1, 1, -1, "Adaptive regularization for barrier method: -1 automatic (default behavior), 0 disabled, 1 enabled"}, // QCQP (barrier) scaling hyper-parameter {CUOPT_QCQP_HYPER_RUIZ_EQUILIBRATION, &pdlp_settings.qcqp_ruiz_equilibration, -1, 1, -1, "Ruiz equilibration for QCQP barrier scaling: -1 automatic (row/column imbalance heuristic), 0 disabled, 1 enabled"}, diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 38202c6b51..46eaa36f8c 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -543,6 +543,11 @@ std::tuple, simplex::lp_status_t, f_t, f_t, f_t barrier_settings.postsolve_info = settings.postsolve_info; barrier_settings.barrier_presolve_bound_free_variables = settings.barrier_presolve_bound_free_variables; + barrier_settings.barrier_presolve_free_elimination = settings.barrier_presolve_free_elimination; + barrier_settings.barrier_free_elimination_row_pivot_tol = + settings.barrier_free_elimination_row_pivot_tol; + barrier_settings.barrier_free_elimination_col_pivot_tol = + settings.barrier_free_elimination_col_pivot_tol; barrier_settings.barrier_initial_point_safeguard = settings.barrier_initial_point_safeguard; barrier_settings.barrier = true; barrier_settings.barrier_presolve = true; diff --git a/cpp/tests/socp/solve_barrier_socp.cu b/cpp/tests/socp/solve_barrier_socp.cu index 0009471974..2bcab7d666 100644 --- a/cpp/tests/socp/solve_barrier_socp.cu +++ b/cpp/tests/socp/solve_barrier_socp.cu @@ -308,6 +308,81 @@ TEST(barrier, presolve_keeps_direct_free_variables_before_cones) EXPECT_EQ(presolve_info.free_elimination_remaining_variables[2], 4); } +TEST(barrier, presolve_skips_free_elimination_on_unstable_pivot) +{ + // Layout: [x0, x1 | cone x2, x3, x4] with only x1 free and zero-cost. Both rows holding x1 + // carry it with a 1e-10 coefficient next to O(1) entries, so substituting it would scatter + // the pivot row amplified by 1e10. Every candidate pivot must fail the row threshold test + // and x1 must survive as a direct free variable. + raft::handle_t handle{}; + init_handler(&handle); + + user_problem_t user_problem(&handle); + + constexpr int m = 2; + constexpr int n = 5; + constexpr int nz = 6; + constexpr double tiny = 1e-10; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {0.0, 0.0, 0.0, 0.0, 0.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + // x0 + tiny*x1 + x2 = 1, tiny*x1 + x3 + x4 = 1 + user_problem.A.col_start = {0, 1, 3, 4, 5, 6}; + const std::vector rows_of_entries = {0, 0, 1, 0, 1, 1}; + const std::vector entry_values = {1.0, tiny, tiny, 1.0, 1.0, 1.0}; + for (int p = 0; p < nz; ++p) { + user_problem.A.i[p] = rows_of_entries[p]; + user_problem.A.x[p] = entry_values[p]; + } + + user_problem.rhs = {1.0, 1.0}; + user_problem.row_sense = {'E', 'E'}; + user_problem.lower = {0.0, -inf, 0.0, 0.0, 0.0}; + user_problem.upper.assign(n, inf); + user_problem.num_range_rows = 0; + user_problem.cone_var_start = 2; + user_problem.second_order_cone_dims = {3}; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + settings.scale_columns = false; + + std::vector new_slacks; + dualize_info_t dualize_info; + lp_problem_t original_lp(user_problem.handle_ptr, 1, 1, 1); + convert_user_problem(user_problem, settings, original_lp, new_slacks, dualize_info); + + presolve_info_t presolve_info; + lp_problem_t presolved_lp(user_problem.handle_ptr, 1, 1, 1); + ASSERT_EQ(presolve(original_lp, settings, presolved_lp, presolve_info), 0); + + EXPECT_EQ(presolved_lp.num_rows, m); + EXPECT_EQ(presolved_lp.num_cols, n); + EXPECT_EQ(presolved_lp.cone_var_start, 2); + EXPECT_TRUE(presolve_info.free_variable_eliminations.empty()); + ASSERT_EQ(presolve_info.direct_free_variables.size(), 1u); + EXPECT_EQ(presolve_info.direct_free_variables[0], 1); + + // Without the thresholds the same column is substituted, which is what the default rejects. + settings.barrier_free_elimination_row_pivot_tol = 0.0; + settings.barrier_free_elimination_col_pivot_tol = 0.0; + presolve_info_t unguarded_info; + lp_problem_t unguarded_lp(user_problem.handle_ptr, 1, 1, 1); + ASSERT_EQ(presolve(original_lp, settings, unguarded_lp, unguarded_info), 0); + ASSERT_EQ(unguarded_info.free_variable_eliminations.size(), 1u); + EXPECT_EQ(unguarded_info.free_variable_eliminations[0].variable, 1); + EXPECT_TRUE(unguarded_info.direct_free_variables.empty()); +} + TEST(barrier, rejects_middle_cone_input_before_barrier) { raft::handle_t handle{}; From d397bbdc893f18cb69aebeab4de61692f778a158 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Tue, 22 Sep 2026 11:38:11 -0700 Subject: [PATCH 11/23] Hardcode the free-variable elimination pivot thresholds. Keep the largest-in-row and largest-in-column tests without exposing them as solver settings. --- .../mathematical_optimization/constants.h | 5 --- .../pdlp/solver_settings.hpp | 7 ---- cpp/src/dual_simplex/presolve.cpp | 36 +++++++++---------- .../dual_simplex/simplex_solver_settings.hpp | 9 ----- cpp/src/math_optimization/solver_settings.cu | 3 -- cpp/src/pdlp/solve.cu | 5 --- cpp/tests/socp/solve_barrier_socp.cu | 10 ------ 7 files changed, 16 insertions(+), 59 deletions(-) diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index ae1d6153c1..3656791a98 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -161,11 +161,6 @@ /* @brief Barrier initial point safeguard */ #define CUOPT_BARRIER_INITIAL_POINT_SAFEGUARD "barrier_initial_point_safeguard" -/* @brief Free-variable equality substitution in barrier presolve */ -#define CUOPT_BARRIER_PRESOLVE_FREE_ELIMINATION "barrier_presolve_free_elimination" -#define CUOPT_BARRIER_FREE_ELIMINATION_ROW_PIVOT_TOL "barrier_free_elimination_row_pivot_tol" -#define CUOPT_BARRIER_FREE_ELIMINATION_COL_PIVOT_TOL "barrier_free_elimination_col_pivot_tol" - /* @brief MIP determinism mode constants */ #define CUOPT_MODE_OPPORTUNISTIC 0 #define CUOPT_MODE_DETERMINISTIC 1 diff --git a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp index bc1040c3f3..fd6c497850 100644 --- a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp @@ -303,13 +303,6 @@ class pdlp_solver_settings_t { barrier_dual_initial_point_t barrier_dual_initial_point{barrier_dual_initial_point_t::Automatic}; i_t postsolve_info{-1}; i_t barrier_presolve_bound_free_variables{-1}; // -1 automatic, 0 disabled, 1 enabled - i_t barrier_presolve_free_elimination{-1}; // -1 automatic, 0 disabled, 1 enabled - // Threshold pivoting tolerances for the free-variable equality substitution. The pivot - // a_pj must be at least this fraction of the largest other entry in its row and of the - // largest entry in its column; 0 disables the corresponding test. At 1.0 no substitution - // can amplify a coefficient. - f_t barrier_free_elimination_row_pivot_tol{1.0}; - f_t barrier_free_elimination_col_pivot_tol{1.0}; // Ruiz equilibration for QCQP (barrier) scaling: -1 automatic (row/column // imbalance heuristic), 0 disabled, 1 enabled. Distinct from PDLP's own Ruiz // scaling in pdlp_hyper_params_t. diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index b96e7492de..30e58afa3a 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -353,8 +353,7 @@ struct substitution_matrix_t { // Returns the number of columns left alone because no incident row passed the pivot threshold. template static i_t eliminate_free_variables(lp_problem_t& problem, - presolve_info_t& presolve_info, - const simplex_solver_settings_t& settings) + presolve_info_t& presolve_info) { const i_t old_m = problem.num_rows; const i_t old_n = problem.num_cols; @@ -385,9 +384,11 @@ static i_t eliminate_free_variables(lp_problem_t& problem, std::vector incident_value; std::vector candidate_order; - const f_t row_pivot_tol = settings.barrier_free_elimination_row_pivot_tol; - const f_t col_pivot_tol = settings.barrier_free_elimination_col_pivot_tol; - i_t pivot_rejected = 0; + // The pivot must be the largest entry of both its row and its column, so the + // substitution cannot amplify a coefficient. + constexpr f_t row_pivot_tol = 1.0; + constexpr f_t col_pivot_tol = 1.0; + i_t pivot_rejected = 0; // One pass in column order. Peeling a chain from one end keeps each pivot row sparse, so // revisiting columns buys almost nothing and costs a requeue storm on models with @@ -455,15 +456,13 @@ static i_t eliminate_free_variables(lp_problem_t& problem, for (const i_t slot : candidate_order) { const f_t a_ij = std::abs(incident_value[slot]); if (a_ij < col_pivot_tol * column_max) { continue; } - if (row_pivot_tol > 0) { - const i_t row = incident[slot]; - f_t row_max = 0; - for (i_t k = 0; k < matrix.row_len[row]; ++k) { - if (matrix.column(row, k) == j) { continue; } - row_max = std::max(row_max, std::abs(matrix.value(row, k))); - } - if (a_ij < row_pivot_tol * row_max) { continue; } + const i_t row = incident[slot]; + f_t row_max = 0; + for (i_t k = 0; k < matrix.row_len[row]; ++k) { + if (matrix.column(row, k) == j) { continue; } + row_max = std::max(row_max, std::abs(matrix.value(row, k))); } + if (a_ij < row_pivot_tol * row_max) { continue; } pivot_slot = slot; break; } @@ -2029,22 +2028,19 @@ i_t presolve(const lp_problem_t& original, } // LP already goes through PSLP; this substitution is for QP/SOCP only. - if (settings.barrier_presolve && settings.barrier_presolve_free_elimination != 0 && - (has_cones || problem.Q.n > 0)) { + if (settings.barrier_presolve && (has_cones || problem.Q.n > 0)) { const i_t old_free_count = static_cast(presolve_info.direct_free_variables.size()); const f_t free_elimination_start = tic(); - const i_t pivot_rejected = eliminate_free_variables(problem, presolve_info, settings); + const i_t pivot_rejected = eliminate_free_variables(problem, presolve_info); const i_t eliminated = old_free_count - static_cast(presolve_info.direct_free_variables.size()); if (eliminated > 0 || pivot_rejected > 0) { settings.log.printf( "Eliminated %d free variables by equality substitution in %.2fs (%d skipped by pivot " - "threshold, row tol %g, col tol %g)\n", + "threshold)\n", eliminated, toc(free_elimination_start), - pivot_rejected, - settings.barrier_free_elimination_row_pivot_tol, - settings.barrier_free_elimination_col_pivot_tol); + pivot_rejected); } } diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 45e5b28fa9..8b3eba56d3 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -84,9 +84,6 @@ struct simplex_solver_settings_t { barrier_dual_initial_point(barrier_dual_initial_point_t::Automatic), postsolve_info(-1), barrier_presolve_bound_free_variables(-1), - barrier_presolve_free_elimination(-1), - barrier_free_elimination_row_pivot_tol(1.0), - barrier_free_elimination_col_pivot_tol(1.0), qcqp_ruiz_equilibration(-1), barrier_initial_point_safeguard(10.0), check_Q(false), @@ -196,12 +193,6 @@ struct simplex_solver_settings_t { // 1 dual least squares, 2 SeDuMi mu-based i_t postsolve_info; // -1 automatic (disabled), 0 disabled, 1 enabled i_t barrier_presolve_bound_free_variables; // -1 automatic, 0 disabled, 1 enabled - i_t barrier_presolve_free_elimination; // -1 automatic, 0 disabled, 1 enabled - // Threshold pivoting tolerances for the free-variable equality substitution: the pivot must - // be at least this fraction of the largest other entry of its row / of the largest entry of - // its column. 0 disables a test, 1 forbids any coefficient growth. - f_t barrier_free_elimination_row_pivot_tol; - f_t barrier_free_elimination_col_pivot_tol; i_t qcqp_ruiz_equilibration; // -1 automatic (imbalance heuristic), 0 disabled, 1 enabled f_t barrier_initial_point_safeguard; // margin pushing the barrier initial iterate into // the interior of the nonnegative orthant / SOC diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 2e92514bc1..5a4ab72c32 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -145,8 +145,6 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_BARRIER_DUAL_REGULARIZATION, &pdlp_settings.barrier_dual_regularization, f_t(-1.0), std::numeric_limits::infinity(), f_t(-1.0), "initial dual regularization for the augmented system; -1 automatic"}, {CUOPT_BARRIER_STEP_SCALE, &pdlp_settings.barrier_step_scale, f_t(0.5), f_t(0.9999), f_t(0.9)}, {CUOPT_BARRIER_INITIAL_POINT_SAFEGUARD, &pdlp_settings.barrier_initial_point_safeguard, f_t(0.0), std::numeric_limits::infinity(), f_t(10.0), "margin pushing the barrier initial iterate into the interior of the nonnegative orthant / SOC"}, - {CUOPT_BARRIER_FREE_ELIMINATION_ROW_PIVOT_TOL, &pdlp_settings.barrier_free_elimination_row_pivot_tol, f_t(0.0), f_t(1.0), f_t(1.0), "free-variable elimination: pivot must be this fraction of the largest other entry in its row; 0 disables the test"}, - {CUOPT_BARRIER_FREE_ELIMINATION_COL_PIVOT_TOL, &pdlp_settings.barrier_free_elimination_col_pivot_tol, f_t(0.0), f_t(1.0), f_t(1.0), "free-variable elimination: pivot must be this fraction of the largest entry in its column; 0 disables the test"}, // MIP heuristic hyper-parameters (hidden from default --help: name contains "hyper_") {CUOPT_MIP_HYPER_HEURISTIC_ROOT_LP_TIME_RATIO, &mip_settings.heuristic_params.root_lp_time_ratio, f_t(0.0), f_t(1.0), f_t(0.1), "fraction of total time for root LP"}, {CUOPT_MIP_HYPER_HEURISTIC_ROOT_LP_MAX_TIME, &mip_settings.heuristic_params.root_lp_max_time, f_t(0.0), std::numeric_limits::infinity(), f_t(15.0), "hard cap on root LP seconds"}, @@ -234,7 +232,6 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_HYPER_SUBMIP_ITERATION_LIMIT_OFFSET, &mip_settings.submip_params.iteration_limit_offset, 0, std::numeric_limits::max(), 10000, "base sub-MIP simplex-iteration limit for root heuristics"}, {CUOPT_MIP_HYPER_SUBMIP_MAX_LEVEL, &mip_settings.submip_params.max_level, 0, std::numeric_limits::max(), 10, "maximum sub-MIP recursion level"}, {CUOPT_BARRIER_PRESOLVE_BOUND_FREE_VARIABLES, &pdlp_settings.barrier_presolve_bound_free_variables, -1, 1, -1, "Bound free variables during barrier presolve: -1 automatic (default behavior), 0 disabled, 1 enabled"}, - {CUOPT_BARRIER_PRESOLVE_FREE_ELIMINATION, &pdlp_settings.barrier_presolve_free_elimination, -1, 1, -1, "Eliminate zero-cost free variables by equality substitution during barrier presolve: -1 automatic (default behavior), 0 disabled, 1 enabled"}, {CUOPT_BARRIER_ADAPTIVE_REGULARIZATION, &pdlp_settings.barrier_adaptive_regularization, -1, 1, -1, "Adaptive regularization for barrier method: -1 automatic (default behavior), 0 disabled, 1 enabled"}, // QCQP (barrier) scaling hyper-parameter {CUOPT_QCQP_HYPER_RUIZ_EQUILIBRATION, &pdlp_settings.qcqp_ruiz_equilibration, -1, 1, -1, "Ruiz equilibration for QCQP barrier scaling: -1 automatic (row/column imbalance heuristic), 0 disabled, 1 enabled"}, diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 46eaa36f8c..38202c6b51 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -543,11 +543,6 @@ std::tuple, simplex::lp_status_t, f_t, f_t, f_t barrier_settings.postsolve_info = settings.postsolve_info; barrier_settings.barrier_presolve_bound_free_variables = settings.barrier_presolve_bound_free_variables; - barrier_settings.barrier_presolve_free_elimination = settings.barrier_presolve_free_elimination; - barrier_settings.barrier_free_elimination_row_pivot_tol = - settings.barrier_free_elimination_row_pivot_tol; - barrier_settings.barrier_free_elimination_col_pivot_tol = - settings.barrier_free_elimination_col_pivot_tol; barrier_settings.barrier_initial_point_safeguard = settings.barrier_initial_point_safeguard; barrier_settings.barrier = true; barrier_settings.barrier_presolve = true; diff --git a/cpp/tests/socp/solve_barrier_socp.cu b/cpp/tests/socp/solve_barrier_socp.cu index 2bcab7d666..9f0c3bd620 100644 --- a/cpp/tests/socp/solve_barrier_socp.cu +++ b/cpp/tests/socp/solve_barrier_socp.cu @@ -371,16 +371,6 @@ TEST(barrier, presolve_skips_free_elimination_on_unstable_pivot) EXPECT_TRUE(presolve_info.free_variable_eliminations.empty()); ASSERT_EQ(presolve_info.direct_free_variables.size(), 1u); EXPECT_EQ(presolve_info.direct_free_variables[0], 1); - - // Without the thresholds the same column is substituted, which is what the default rejects. - settings.barrier_free_elimination_row_pivot_tol = 0.0; - settings.barrier_free_elimination_col_pivot_tol = 0.0; - presolve_info_t unguarded_info; - lp_problem_t unguarded_lp(user_problem.handle_ptr, 1, 1, 1); - ASSERT_EQ(presolve(original_lp, settings, unguarded_lp, unguarded_info), 0); - ASSERT_EQ(unguarded_info.free_variable_eliminations.size(), 1u); - EXPECT_EQ(unguarded_info.free_variable_eliminations[0].variable, 1); - EXPECT_TRUE(unguarded_info.direct_free_variables.empty()); } TEST(barrier, rejects_middle_cone_input_before_barrier) From 79a2898d0eda318012281fdac7e5cd9355b7fc37 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Tue, 22 Sep 2026 15:34:21 -0700 Subject: [PATCH 12/23] Fix an indexing bug --- cpp/src/dual_simplex/presolve.cpp | 37 ++++++++++++++++++------------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 30e58afa3a..0564816581 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -233,8 +233,8 @@ struct substitution_matrix_t { // erase and insert both require row i to be resident. void erase(i_t i, i_t offset) { - const i_t start = row_start[i]; - const i_t last = row_len[i] - 1; + const i_t start = row_start[i]; + const i_t last = row_len[i] - 1; scatter[row_col[start + offset]] = -1; if (offset != last) { row_col[start + offset] = row_col[start + last]; @@ -270,8 +270,7 @@ struct substitution_matrix_t { } if (row_used > 2 * (live + num_rows)) { compact_rows(); } if (row_used + need > static_cast(row_col.size())) { - const i_t size = - std::max(row_used + need, 2 * static_cast(row_col.size()) + 1); + const i_t size = std::max(row_used + need, 2 * static_cast(row_col.size()) + 1); row_col.resize(size, 0); row_val.resize(size, 0); } @@ -295,8 +294,7 @@ struct substitution_matrix_t { } if (col_used > 2 * (live + num_cols)) { compact_cols(); } if (col_used + need > static_cast(col_row.size())) { - const i_t size = - std::max(col_used + need, 2 * static_cast(col_row.size()) + 1); + const i_t size = std::max(col_used + need, 2 * static_cast(col_row.size()) + 1); col_row.resize(size, 0); } if (col_len[j] < col_cap[j]) { return; } @@ -558,8 +556,8 @@ static i_t eliminate_free_variables(lp_problem_t& problem, if (active_row[i]) { remaining_rows.push_back(i); } } - i_t new_n = static_cast(remaining_cols.size()); - i_t new_m = static_cast(remaining_rows.size()); + i_t new_n = static_cast(remaining_cols.size()); + i_t new_m = static_cast(remaining_rows.size()); matrix.unload(); i_t new_nnz = 0; for (const i_t i : remaining_rows) { @@ -605,16 +603,23 @@ static i_t eliminate_free_variables(lp_problem_t& problem, if (problem.Q.n > 0) { remove_variables_from_Q(problem.Q, col_marker, old_to_new_col, new_n); } reduced_A.to_compressed_col(problem.A); - problem.rhs = std::move(reduced_rhs); - problem.objective = std::move(objective); - problem.lower = std::move(lower); - problem.upper = std::move(upper); - problem.num_rows = new_m; - problem.num_cols = new_n; - problem.cone_var_start = old_to_new_col[problem.cone_var_start]; + problem.rhs = std::move(reduced_rhs); + problem.objective = std::move(objective); + problem.lower = std::move(lower); + problem.upper = std::move(upper); + problem.num_rows = new_m; + problem.num_cols = new_n; + // Cone columns are never eliminated, so the cone block stays trailing and its new start is + // just the remapped old one. Without cones cone_var_start is 0 and must be left alone. + if (!problem.second_order_cone_dims.empty()) { + const i_t new_cone_start = old_to_new_col[problem.cone_var_start]; + assert(new_cone_start != -1); + problem.cone_var_start = new_cone_start; + } presolve_info.direct_free_variables.clear(); - for (i_t new_j = 0; new_j < problem.cone_var_start; ++new_j) { + const i_t new_linear_cols = linear_variable_count(problem); + for (i_t new_j = 0; new_j < new_linear_cols; ++new_j) { if (problem.lower[new_j] == -inf && problem.upper[new_j] == inf) { presolve_info.direct_free_variables.push_back(new_j); } From f8668f6dd398f130bea8fcf60027da318ee41665 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Tue, 22 Sep 2026 15:50:04 -0700 Subject: [PATCH 13/23] Use TMPDIR instead of explicit tmp --- benchmarks/linear_programming/run_mps_files.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/linear_programming/run_mps_files.sh b/benchmarks/linear_programming/run_mps_files.sh index d8625a0f2a..5e411c77e9 100755 --- a/benchmarks/linear_programming/run_mps_files.sh +++ b/benchmarks/linear_programming/run_mps_files.sh @@ -383,7 +383,7 @@ mps_files=("${mps_files[@]:$start_idx:$((end_idx-start_idx))}") file_count=${#mps_files[@]} # Initialize the index file for locking mechanism -INDEX_FILE="/tmp/mps_file_index.$$" +INDEX_FILE="${TMPDIR:-/tmp}/mps_file_index.$$" # Remove the index file if it exists rm -f "$INDEX_FILE" From 3cd2e17d10a0d5bac0ce8bcb47ef90592b04e1f9 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Tue, 22 Sep 2026 21:27:29 -0700 Subject: [PATCH 14/23] Resize uncrush vectors appropriately --- cpp/tests/socp/solve_barrier_socp.cu | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/tests/socp/solve_barrier_socp.cu b/cpp/tests/socp/solve_barrier_socp.cu index 9f0c3bd620..d098658528 100644 --- a/cpp/tests/socp/solve_barrier_socp.cu +++ b/cpp/tests/socp/solve_barrier_socp.cu @@ -1212,9 +1212,9 @@ TEST(barrier, free_variable_substitution_postsolve_kkt) stationarity_residual(presolved_lp, crushed_x, crushed_y, crushed_z, reduced_stationarity); ASSERT_NEAR(inf_norm(reduced_stationarity), 0.0, 1e-12); - std::vector uncrushed_x; - std::vector uncrushed_y; - std::vector uncrushed_z; + std::vector uncrushed_x(original_lp.num_cols); + std::vector uncrushed_y(original_lp.num_rows); + std::vector uncrushed_z(original_lp.num_cols); uncrush_solution(presolve_info, settings, original_lp, From 8dd5d19e6c9026d2d3ca9dfb4c1a97faeb321e01 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Tue, 22 Sep 2026 22:33:51 -0700 Subject: [PATCH 15/23] Fix CI failures --- cpp/src/grpc/codegen/generated/cuopt_mcp_schema.json | 4 ++++ cpp/src/grpc/codegen/generated/cuopt_remote_data.proto | 2 -- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/cpp/src/grpc/codegen/generated/cuopt_mcp_schema.json b/cpp/src/grpc/codegen/generated/cuopt_mcp_schema.json index ca36b5c8e7..d1d13a537b 100644 --- a/cpp/src/grpc/codegen/generated/cuopt_mcp_schema.json +++ b/cpp/src/grpc/codegen/generated/cuopt_mcp_schema.json @@ -170,6 +170,10 @@ "pdlp_precision": { "type": "integer", "description": "Precision mode for the PDLP solver. DefaultPrecision uses the problem's native precision; SinglePrecision runs PDHG in FP32 (half the memory, roughly 2x faster iterations, possibly more of them); DoublePrecision forces FP64; MixedPrecision stores the constraint matrix in FP32 for faster SpMV while keeping vectors and compute in FP64 (convergence checks still use the FP64 matrix, so memory is not reduced). Default: DefaultPrecision." + }, + "do_curtis_reid_scaling": { + "type": "boolean", + "description": "Whether Curtis-Reid prescaling runs before Ruiz/Pock-Chambolle matrix scaling. Default: true." } }, "additionalProperties": false diff --git a/cpp/src/grpc/codegen/generated/cuopt_remote_data.proto b/cpp/src/grpc/codegen/generated/cuopt_remote_data.proto index 3c6338f718..0b48826dcc 100644 --- a/cpp/src/grpc/codegen/generated/cuopt_remote_data.proto +++ b/cpp/src/grpc/codegen/generated/cuopt_remote_data.proto @@ -198,8 +198,6 @@ message PDLPSolverSettings { optional double barrier_step_scale = 32; optional int32 postsolve_info = 33; optional int32 barrier_adaptive_regularization = 34; - // Whether Curtis-Reid prescaling runs before Ruiz/Pock-Chambolle matrix - // scaling. (default: true) optional bool do_curtis_reid_scaling = 35; PDLPWarmStartData warm_start_data = 50; } From e3e1e0f9d573eeb70af33f0d8a3917b3bdfa6cc0 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 23 Sep 2026 07:51:39 -0700 Subject: [PATCH 16/23] Disable CR scaling for PDLP test --- cpp/tests/linear_programming/pdlp_test.cu | 3 +++ 1 file changed, 3 insertions(+) diff --git a/cpp/tests/linear_programming/pdlp_test.cu b/cpp/tests/linear_programming/pdlp_test.cu index b9d3e72fb0..7e2b74db3b 100644 --- a/cpp/tests/linear_programming/pdlp_test.cu +++ b/cpp/tests/linear_programming/pdlp_test.cu @@ -3773,6 +3773,8 @@ TEST(pdlp_class, run_batch_pdlp_many_different_bounds) regular_pdlp_settings.method = cuopt::mathematical_optimization::method_t::PDLP; regular_pdlp_settings.pdlp_solver_mode = pdlp_solver_mode_t::Stable3; regular_pdlp_settings.presolver = presolver_t::None; + // Known issue with Curtis-Reid scaling on batch PDLP. + regular_pdlp_settings.hyper_params.do_curtis_reid_scaling = false; regular_pdlp_settings.set_optimality_tolerance(result_tolerance); const std::vector>> bound_offsets_by_climber = { @@ -3836,6 +3838,7 @@ TEST(pdlp_class, run_batch_pdlp_many_different_bounds) auto batch_settings = regular_pdlp_settings; batch_settings.generate_batch_primal_dual_solution = true; + batch_settings.hyper_params.do_curtis_reid_scaling = false; for (int i = 0; i < batch_size; ++i) { for (const auto& bounds : custom_bounds_by_climber[i]) { batch_settings.new_bounds.push_back( From 47e106dd9bda03567ddcb6b3341692838d9e77e7 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 23 Sep 2026 08:41:54 -0700 Subject: [PATCH 17/23] Address PR review --- cpp/src/dual_simplex/solve.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 84121b7fb4..fa9778511b 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -117,11 +117,11 @@ void write_matlab(const std::string& filename, const simplex::lp_problem_t -void compute_stationarity_residual(const lp_problem_t& lp, - const std::vector& x, - const std::vector& y, - const std::vector& z, - std::vector& residual) +void compute_dual_residual(const lp_problem_t& lp, + const std::vector& x, + const std::vector& y, + const std::vector& z, + std::vector& residual) { residual = z; for (i_t j = 0; j < lp.num_cols; ++j) { @@ -604,7 +604,7 @@ lp_status_t solve_linear_program_with_barrier( primal_residual, primal_residual / (1.0 + vector_norm_inf(presolved_lp.rhs))); std::vector unscaled_dual_residual; - compute_stationarity_residual( + compute_dual_residual( presolved_lp, unscaled_x, unscaled_y, unscaled_z, unscaled_dual_residual); f_t unscaled_dual_residual_norm = vector_norm_inf(unscaled_dual_residual); settings.log.printf( @@ -634,7 +634,7 @@ lp_status_t solve_linear_program_with_barrier( post_solve_primal_residual / (1.0 + vector_norm_inf(original_lp.rhs))); std::vector post_solve_dual_residual; - compute_stationarity_residual( + compute_dual_residual( original_lp, lp_solution.x, lp_solution.y, lp_solution.z, post_solve_dual_residual); f_t post_solve_dual_residual_norm = vector_norm_inf(post_solve_dual_residual); settings.log.printf( From 9625c54aabef7cf0e080eafa98419fe402fbd97b Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 23 Sep 2026 12:04:53 -0700 Subject: [PATCH 18/23] disable CR scaling for warm start test --- cpp/tests/linear_programming/unit_tests/presolve_test.cu | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cpp/tests/linear_programming/unit_tests/presolve_test.cu b/cpp/tests/linear_programming/unit_tests/presolve_test.cu index f8ba6432fe..57c61b2aac 100644 --- a/cpp/tests/linear_programming/unit_tests/presolve_test.cu +++ b/cpp/tests/linear_programming/unit_tests/presolve_test.cu @@ -929,6 +929,8 @@ TEST_P(crush_warmstart, round_trip) settings.dual_postsolve = true; settings.method = cuopt::mathematical_optimization::method_t::PDLP; settings.time_limit = 60.0; + // Known issue with Curtis-Reid scaling on PDLP warm starts. + settings.hyper_params.do_curtis_reid_scaling = false; auto cold_solution = solve_lp(result.reduced_problem, settings); ASSERT_EQ(cold_solution.get_termination_status(), pdlp_termination_status_t::Optimal); From 97065f7dba6ab39e2d024c730ca36916ca3ab0b1 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 23 Sep 2026 12:42:55 -0700 Subject: [PATCH 19/23] Address remaining PR review comments. --- cpp/src/dual_simplex/presolve.cpp | 60 +++++++++++------------- cpp/src/dual_simplex/solve.cpp | 8 +--- cpp/src/linear_algebra/sparse_matrix.hpp | 29 ++++++++++++ cpp/tests/socp/solve_barrier_socp.cu | 52 ++++++++------------ 4 files changed, 77 insertions(+), 72 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 0564816581..1b27e53ca0 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -17,7 +17,7 @@ #include #include #include -#include +#include namespace cuopt::mathematical_optimization::simplex { @@ -380,7 +380,6 @@ static i_t eliminate_free_variables(lp_problem_t& problem, i_t stamp = 0; std::vector incident; std::vector incident_value; - std::vector candidate_order; // The pivot must be the largest entry of both its row and its column, so the // substitution cannot amplify a coefficient. @@ -428,41 +427,38 @@ static i_t eliminate_free_variables(lp_problem_t& problem, continue; } - // A small pivot row limits fill. On a path this peels sparse boundary rows instead of - // repeatedly traversing the growing aggregate row. Candidates are ranked shortest first - // and the first one passing the threshold tests wins, so a long aggregate row is only - // scanned once every shorter row has failed. - candidate_order.resize(incident.size()); - std::iota(candidate_order.begin(), candidate_order.end(), 0); - std::sort(candidate_order.begin(), candidate_order.end(), [&](i_t a, i_t b) { - const i_t len_a = matrix.row_len[incident[a]]; - const i_t len_b = matrix.row_len[incident[b]]; - if (len_a != len_b) { return len_a < len_b; } - return std::abs(incident_value[a]) > std::abs(incident_value[b]); - }); - // Threshold pivoting, as in the LU: a pivot that is tiny relative to its own row or to the // rest of its column makes the multiplier a_ij / a_pj huge and inflates every surviving row // by that factor, which the barrier's KKT solve cannot recover from. The column test bounds // the multiplier, the row test bounds the coefficients the pivot row scatters. At tolerance // 1 the pivot is the largest entry of both its row and its column, so no coefficient grows. + // Walk incident rows by increasing length (as in the LU degree walk) so a chain peels from + // a sparse end instead of aggregating into a dense row. f_t column_max = 0; - for (const f_t value : incident_value) { - column_max = std::max(column_max, std::abs(value)); + i_t min_len = std::numeric_limits::max(); + i_t max_len = 0; + for (size_t slot = 0; slot < incident.size(); ++slot) { + column_max = std::max(column_max, std::abs(incident_value[slot])); + const i_t len = matrix.row_len[incident[slot]]; + min_len = std::min(min_len, len); + max_len = std::max(max_len, len); } i_t pivot_slot = -1; - for (const i_t slot : candidate_order) { - const f_t a_ij = std::abs(incident_value[slot]); - if (a_ij < col_pivot_tol * column_max) { continue; } - const i_t row = incident[slot]; - f_t row_max = 0; - for (i_t k = 0; k < matrix.row_len[row]; ++k) { - if (matrix.column(row, k) == j) { continue; } - row_max = std::max(row_max, std::abs(matrix.value(row, k))); + for (i_t len = min_len; len <= max_len && pivot_slot == -1; ++len) { + for (size_t slot = 0; slot < incident.size(); ++slot) { + if (matrix.row_len[incident[slot]] != len) { continue; } + const f_t a_ij = std::abs(incident_value[slot]); + if (a_ij < col_pivot_tol * column_max) { continue; } + const i_t row = incident[slot]; + f_t row_max = 0; + for (i_t k = 0; k < matrix.row_len[row]; ++k) { + if (matrix.column(row, k) == j) { continue; } + row_max = std::max(row_max, std::abs(matrix.value(row, k))); + } + if (a_ij < row_pivot_tol * row_max) { continue; } + pivot_slot = static_cast(slot); + break; } - if (a_ij < row_pivot_tol * row_max) { continue; } - pivot_slot = slot; - break; } // No stable pivot: leave the column as a free variable handled directly in the KKT system. if (pivot_slot == -1) { @@ -2041,11 +2037,11 @@ i_t presolve(const lp_problem_t& original, old_free_count - static_cast(presolve_info.direct_free_variables.size()); if (eliminated > 0 || pivot_rejected > 0) { settings.log.printf( - "Eliminated %d free variables by equality substitution in %.2fs (%d skipped by pivot " - "threshold)\n", + "Eliminated %d free variables by equality substitution (%d skipped by pivot " + "threshold) in %.2fs\n", eliminated, - toc(free_elimination_start), - pivot_rejected); + pivot_rejected, + toc(free_elimination_start)); } } diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index fa9778511b..39931bb42a 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -127,13 +127,7 @@ void compute_dual_residual(const lp_problem_t& lp, for (i_t j = 0; j < lp.num_cols; ++j) { residual[j] -= lp.objective[j]; } - if (lp.Q.n > 0) { - for (i_t i = 0; i < lp.Q.m; ++i) { - for (i_t p = lp.Q.row_start[i]; p < lp.Q.row_start[i + 1]; ++p) { - residual[i] -= lp.Q.x[p] * x[lp.Q.j[p]]; - } - } - } + if (lp.Q.n > 0) { matrix_vector_multiply(lp.Q, -1.0, x, 1.0, residual); } matrix_transpose_vector_multiply(lp.A, 1.0, y, 1.0, residual); } diff --git a/cpp/src/linear_algebra/sparse_matrix.hpp b/cpp/src/linear_algebra/sparse_matrix.hpp index 95ee0f0f32..33937d435d 100644 --- a/cpp/src/linear_algebra/sparse_matrix.hpp +++ b/cpp/src/linear_algebra/sparse_matrix.hpp @@ -311,4 +311,33 @@ i_t matrix_vector_multiply( return 0; } +// y <- alpha*A*x + beta*y +template +i_t matrix_vector_multiply( + const csr_matrix_t& A, f_t alpha, const VectorX& x, f_t beta, VectorY& y) +{ + const i_t m = A.m; + const i_t n = A.n; + assert(y.size() == static_cast(m)); + assert(x.size() == static_cast(n)); + + if (beta != 1.0) { + for (i_t i = 0; i < m; ++i) { + y[i] *= beta; + } + } + + for (i_t i = 0; i < m; ++i) { + const i_t row_start = A.row_start[i]; + const i_t row_end = A.row_start[i + 1]; + f_t dot = 0.0; + for (i_t p = row_start; p < row_end; ++p) { + dot += A.x[p] * x[A.j[p]]; + } + y[i] += alpha * dot; + } + + return 0; +} + } // namespace cuopt::mathematical_optimization diff --git a/cpp/tests/socp/solve_barrier_socp.cu b/cpp/tests/socp/solve_barrier_socp.cu index d098658528..550a7e3021 100644 --- a/cpp/tests/socp/solve_barrier_socp.cu +++ b/cpp/tests/socp/solve_barrier_socp.cu @@ -14,6 +14,7 @@ #include #include #include +#include #include #include @@ -34,15 +35,6 @@ static void init_handler(const raft::handle_t* handle_ptr) handle_ptr->get_stream().get())); } -static double inf_norm(const std::vector& v) -{ - double nrm = 0.0; - for (double val : v) { - nrm = std::max(nrm, std::abs(val)); - } - return nrm; -} - // Hub-and-spoke style chain: two free integrator variables plus a quadratic on w. // // minimize 0.5 w^2 @@ -93,23 +85,17 @@ static user_problem_t make_free_substitution_qp(raft::handle_t* han return user_problem; } -static void stationarity_residual(const lp_problem_t& lp, - const std::vector& x, - const std::vector& y, - const std::vector& z, - std::vector& residual) +static void dual_residual(const lp_problem_t& lp, + const std::vector& x, + const std::vector& y, + const std::vector& z, + std::vector& residual) { residual = z; for (int j = 0; j < lp.num_cols; ++j) { residual[j] -= lp.objective[j]; } - if (lp.Q.n > 0) { - for (int i = 0; i < lp.Q.m; ++i) { - for (int p = lp.Q.row_start[i]; p < lp.Q.row_start[i + 1]; ++p) { - residual[i] -= lp.Q.x[p] * x[lp.Q.j[p]]; - } - } - } + if (lp.Q.n > 0) { matrix_vector_multiply(lp.Q, -1.0, x, 1.0, residual); } matrix_transpose_vector_multiply(lp.A, 1.0, y, 1.0, residual); } @@ -1204,13 +1190,13 @@ TEST(barrier, free_variable_substitution_postsolve_kkt) } std::vector crushed_y(presolved_lp.num_rows, 0.25); std::vector crushed_z(presolved_lp.num_cols, 0.0); - std::vector reduced_stationarity; - stationarity_residual(presolved_lp, crushed_x, crushed_y, crushed_z, reduced_stationarity); + std::vector reduced_dual; + dual_residual(presolved_lp, crushed_x, crushed_y, crushed_z, reduced_dual); for (int j = 0; j < presolved_lp.num_cols; ++j) { - crushed_z[j] -= reduced_stationarity[j]; + crushed_z[j] -= reduced_dual[j]; } - stationarity_residual(presolved_lp, crushed_x, crushed_y, crushed_z, reduced_stationarity); - ASSERT_NEAR(inf_norm(reduced_stationarity), 0.0, 1e-12); + dual_residual(presolved_lp, crushed_x, crushed_y, crushed_z, reduced_dual); + ASSERT_NEAR((vector_norm_inf(reduced_dual)), 0.0, 1e-12); std::vector uncrushed_x(original_lp.num_cols); std::vector uncrushed_y(original_lp.num_rows); @@ -1235,11 +1221,11 @@ TEST(barrier, free_variable_substitution_postsolve_kkt) std::vector primal_residual = original_lp.rhs; matrix_vector_multiply(original_lp.A, 1.0, uncrushed_x, -1.0, primal_residual); - EXPECT_NEAR(inf_norm(primal_residual), 0.0, 1e-12); + EXPECT_NEAR((vector_norm_inf(primal_residual)), 0.0, 1e-12); - std::vector dual_residual; - stationarity_residual(original_lp, uncrushed_x, uncrushed_y, uncrushed_z, dual_residual); - EXPECT_NEAR(inf_norm(dual_residual), 0.0, 1e-10); + std::vector uncrushed_dual; + dual_residual(original_lp, uncrushed_x, uncrushed_y, uncrushed_z, uncrushed_dual); + EXPECT_NEAR((vector_norm_inf(uncrushed_dual)), 0.0, 1e-10); lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); auto status = solve_linear_program_with_barrier(user_problem, settings, solution); @@ -1254,11 +1240,11 @@ TEST(barrier, free_variable_substitution_postsolve_kkt) std::vector solved_primal = original_lp.rhs; matrix_vector_multiply(original_lp.A, 1.0, solution.x, -1.0, solved_primal); - EXPECT_NEAR(inf_norm(solved_primal), 0.0, 1e-5); + EXPECT_NEAR((vector_norm_inf(solved_primal)), 0.0, 1e-5); std::vector solved_dual; - stationarity_residual(original_lp, solution.x, solution.y, solution.z, solved_dual); - EXPECT_NEAR(inf_norm(solved_dual), 0.0, 1e-5); + dual_residual(original_lp, solution.x, solution.y, solution.z, solved_dual); + EXPECT_NEAR((vector_norm_inf(solved_dual)), 0.0, 1e-5); } } // namespace cuopt::mathematical_optimization::simplex::test From 04e8c462818e59482e2592415e24219a5cb97180 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 23 Sep 2026 14:01:27 -0700 Subject: [PATCH 20/23] Cache row maxima when choosing a substitution pivot. --- cpp/src/dual_simplex/presolve.cpp | 66 ++++++++++++++++++++----------- 1 file changed, 43 insertions(+), 23 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 1b27e53ca0..deac1d7109 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -124,6 +124,7 @@ struct substitution_matrix_t { std::vector col_len; std::vector col_cap; std::vector scatter; // column -> offset within the resident row, -1 when absent + std::vector row_max; // inf-norm of each row, kept up to date on insert/erase/set i_t scattered = -1; // Arena high-water marks. A line that outgrows its slot is relocated to the tail with // twice the capacity, so appends are amortized O(1); the holes left behind are reclaimed @@ -185,11 +186,37 @@ struct substitution_matrix_t { } scatter.assign(num_cols, -1); + row_max.assign(num_rows, 0); + for (i_t i = 0; i < num_rows; ++i) { + recompute_row_max(i); + } } i_t column(i_t i, i_t offset) const { return row_col[row_start[i] + offset]; } f_t value(i_t i, i_t offset) const { return row_val[row_start[i] + offset]; } - void set_value(i_t i, i_t offset, f_t value) { row_val[row_start[i] + offset] = value; } + + void recompute_row_max(i_t i) + { + f_t max_abs = 0; + const i_t start = row_start[i]; + for (i_t k = 0; k < row_len[i]; ++k) { + max_abs = std::max(max_abs, std::abs(row_val[start + k])); + } + row_max[i] = max_abs; + } + + void set_value(i_t i, i_t offset, f_t value) + { + f_t& slot = row_val[row_start[i] + offset]; + const f_t old_abs = std::abs(slot); + slot = value; + const f_t new_abs = std::abs(value); + if (new_abs >= row_max[i]) { + row_max[i] = new_abs; + } else if (old_abs == row_max[i]) { + recompute_row_max(i); + } + } void unload() { @@ -235,6 +262,7 @@ struct substitution_matrix_t { { const i_t start = row_start[i]; const i_t last = row_len[i] - 1; + const f_t old_abs = std::abs(row_val[start + offset]); scatter[row_col[start + offset]] = -1; if (offset != last) { row_col[start + offset] = row_col[start + last]; @@ -242,6 +270,7 @@ struct substitution_matrix_t { scatter[row_col[start + offset]] = offset; } row_len[i] = last; + if (old_abs == row_max[i]) { recompute_row_max(i); } } void insert(i_t i, i_t col, f_t value) @@ -252,6 +281,7 @@ struct substitution_matrix_t { row_val[row_start[i] + offset] = value; scatter[col] = offset; row_len[i] = offset + 1; + row_max[i] = std::max(row_max[i], std::abs(value)); reserve_col(col); col_row[col_start[col] + col_len[col]] = i; ++col_len[col]; @@ -432,32 +462,22 @@ static i_t eliminate_free_variables(lp_problem_t& problem, // by that factor, which the barrier's KKT solve cannot recover from. The column test bounds // the multiplier, the row test bounds the coefficients the pivot row scatters. At tolerance // 1 the pivot is the largest entry of both its row and its column, so no coefficient grows. - // Walk incident rows by increasing length (as in the LU degree walk) so a chain peels from - // a sparse end instead of aggregating into a dense row. + // Among those, pick the shortest row so a chain peels from a sparse end. f_t column_max = 0; - i_t min_len = std::numeric_limits::max(); - i_t max_len = 0; - for (size_t slot = 0; slot < incident.size(); ++slot) { - column_max = std::max(column_max, std::abs(incident_value[slot])); - const i_t len = matrix.row_len[incident[slot]]; - min_len = std::min(min_len, len); - max_len = std::max(max_len, len); + for (const f_t value : incident_value) { + column_max = std::max(column_max, std::abs(value)); } i_t pivot_slot = -1; - for (i_t len = min_len; len <= max_len && pivot_slot == -1; ++len) { - for (size_t slot = 0; slot < incident.size(); ++slot) { - if (matrix.row_len[incident[slot]] != len) { continue; } - const f_t a_ij = std::abs(incident_value[slot]); - if (a_ij < col_pivot_tol * column_max) { continue; } - const i_t row = incident[slot]; - f_t row_max = 0; - for (i_t k = 0; k < matrix.row_len[row]; ++k) { - if (matrix.column(row, k) == j) { continue; } - row_max = std::max(row_max, std::abs(matrix.value(row, k))); - } - if (a_ij < row_pivot_tol * row_max) { continue; } + i_t best_len = std::numeric_limits::max(); + for (size_t slot = 0; slot < incident.size(); ++slot) { + const f_t a_ij = std::abs(incident_value[slot]); + if (a_ij < col_pivot_tol * column_max) { continue; } + const i_t row = incident[slot]; + if (a_ij < row_pivot_tol * matrix.row_max[row]) { continue; } + const i_t len = matrix.row_len[row]; + if (len < best_len) { + best_len = len; pivot_slot = static_cast(slot); - break; } } // No stable pivot: leave the column as a free variable handled directly in the KKT system. From 3c17339620ddbbddd1f73cc1cf563aaa7449f9ab Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 23 Sep 2026 16:32:22 -0700 Subject: [PATCH 21/23] fix style checks --- cpp/src/dual_simplex/presolve.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index deac1d7109..47d879ff42 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -198,6 +198,7 @@ struct substitution_matrix_t { void recompute_row_max(i_t i) { f_t max_abs = 0; + const i_t start = row_start[i]; for (i_t k = 0; k < row_len[i]; ++k) { max_abs = std::max(max_abs, std::abs(row_val[start + k])); @@ -207,10 +208,10 @@ struct substitution_matrix_t { void set_value(i_t i, i_t offset, f_t value) { - f_t& slot = row_val[row_start[i] + offset]; - const f_t old_abs = std::abs(slot); - slot = value; - const f_t new_abs = std::abs(value); + f_t& slot = row_val[row_start[i] + offset]; + const f_t old_abs = std::abs(slot); + slot = value; + const f_t new_abs = std::abs(value); if (new_abs >= row_max[i]) { row_max[i] = new_abs; } else if (old_abs == row_max[i]) { From 1ffe501fc51f58213b2b47d132c649844d2ce7fd Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Mon, 28 Sep 2026 11:33:50 -0700 Subject: [PATCH 22/23] Remove default handling of IR tolerance --- cpp/src/barrier/barrier.cu | 9 +++++---- cpp/src/barrier/iterative_refinement.hpp | 6 +++--- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index dc56fe78e8..d60251f29f 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -3263,8 +3263,9 @@ i_t barrier_solver_t::gpu_compute_search_direction(iteration_data_t(adat_op, data.d_h_, data.d_dy_); + iterative_refinement(adat_op, data.d_h_, data.d_dy_, ir_tol); if (adat_solve_err > 1e-1) { settings.log.debug("||ADAT*dy - h|| %e after IR\n", adat_solve_err); } @@ -4473,15 +4474,15 @@ lp_status_t barrier_solver_t::barrier_advanced_solve(f_t start_time, mu, primal_objective, dual_objective); - f_t user_primal_objective = compute_user_objective(lp, primal_objective); - f_t user_dual_objective = compute_user_objective(lp, dual_objective); + f_t user_primal_objective = compute_user_objective(lp, primal_objective); + f_t user_dual_objective = compute_user_objective(lp, dual_objective); f_t relative_primal_residual = primal_residual_norm / (1.0 + norm_b); f_t relative_dual_residual = dual_residual_norm / (1.0 + norm_c); f_t relative_complementarity_residual = complementarity_residual_norm / (1.0 + std::min(std::abs(user_primal_objective), std::abs(primal_objective))); - + f_t objective_gap, relative_objective_gap; compute_objective_gap( lp, primal_objective, dual_objective, objective_gap, relative_objective_gap); diff --git a/cpp/src/barrier/iterative_refinement.hpp b/cpp/src/barrier/iterative_refinement.hpp index 9391dd0b6d..7680dfe610 100644 --- a/cpp/src/barrier/iterative_refinement.hpp +++ b/cpp/src/barrier/iterative_refinement.hpp @@ -58,7 +58,7 @@ template f_t iterative_refinement_simple(T& op, const rmm::device_uvector& b, rmm::device_uvector& x, - f_t tol = 1e-8) + f_t tol) { rmm::device_uvector x_sav(x, x.stream()); @@ -125,7 +125,7 @@ template f_t iterative_refinement_gmres(T& op, const rmm::device_uvector& b, rmm::device_uvector& x, - f_t tol = 1e-8) + f_t tol) { // Parameters // Ideally, we do not need to restart here. But having restarts helps as a checkpoint to get @@ -383,7 +383,7 @@ template f_t iterative_refinement(T& op, const rmm::device_uvector& b, rmm::device_uvector& x, - f_t tol = 1e-8) + f_t tol) { return iterative_refinement_gmres(op, b, x, tol); } From 6f2eb4933e2c9b743c262facf2a29d20ea93de64 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Mon, 28 Sep 2026 14:03:37 -0700 Subject: [PATCH 23/23] Implement adaptive IR --- cpp/src/barrier/barrier.cu | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index d60251f29f..93de531290 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -2203,6 +2203,8 @@ class iteration_data_t { f_t dual_residual_norm_save; f_t complementarity_residual_norm_save; + bool use_high_accuracy_ir = false; + dense_vector_t diag; pinned_dense_vector_t inv_diag; dense_vector_t inv_sqrt_diag; @@ -2581,7 +2583,8 @@ int barrier_solver_t::initial_point(iteration_data_t& data) } op(data); if (settings.barrier_iterative_refinement) { - const f_t ir_tol = data.has_sparse_cones() ? f_t(1e-12) : f_t(1e-8); + const f_t ir_tol = + (data.has_sparse_cones() || data.use_high_accuracy_ir) ? f_t(1e-12) : f_t(1e-8); iterative_refinement(op, rhs, soln, ir_tol); } @@ -3175,7 +3178,8 @@ i_t barrier_solver_t::gpu_compute_search_direction(iteration_data_t( op, data.d_augmented_rhs_, data.d_augmented_soln_, ir_tol); if (solve_err > 1e-1) { @@ -4647,6 +4651,11 @@ lp_status_t barrier_solver_t::barrier_advanced_solve(f_t start_time, compute_objective_gap( lp, primal_objective, dual_objective, objective_gap, relative_objective_gap); + if (data.has_cones() && !data.use_high_accuracy_ir && relative_primal_residual < 1e-6 && + relative_dual_residual < 1e-6 && relative_complementarity_residual < 1e-6) { + data.use_high_accuracy_ir = true; + } + if (relative_primal_residual < settings.barrier_relaxed_feasibility_tol && relative_dual_residual < settings.barrier_relaxed_optimality_tol && relative_complementarity_residual < settings.barrier_relaxed_complementarity_tol) {