From c62b025f0572ba3a556606b80dedd7e880f410c0 Mon Sep 17 00:00:00 2001 From: yuwenchen95 Date: Mon, 14 Sep 2026 06:51:19 -0700 Subject: [PATCH] Barrier presolve: remove fixed variables Variables with equal lower and upper bounds contribute a constant to the objective and a constant to each row activity, so they can be folded into the objective constant and the right-hand side and dropped. Postsolve restores their value from the recorded bound and recovers their reduced cost from the final duals as z_j = c_j - a_j^T y, which is what dual feasibility A^T y + z = c requires; it cannot be recorded during presolve because y is not known yet. Columns carrying a quadratic objective term are left in place, since a fixed variable's contribution to x'Qx is not a pure constant. Cone columns are also left alone: they must stay a contiguous trailing block and the barrier rejects explicit bounds on them. Two existing problems this depends on, both fixed here: - The number of leading linear columns is now re-measured after the column reductions. It was captured before them and then used to index lower/upper, so it already read out of bounds whenever remove_empty_cols dropped a column, and would have read 167184 entries past the end on the instance below. - The barrier is skipped when presolve leaves no columns at all. Reading the correct column count surfaces that case rather than hiding it behind the out-of-bounds read above. The barrier cannot form a KKT system with no columns, and the solution is already fully determined by postsolve. Empty rows are now removed after the column reductions rather than before, so rows left empty by them are removed too. Removing an empty row cannot change any column's nonzero count, so the previous order gained nothing. Controlled by barrier_presolve_fixed_variables (-1 automatic/enabled, 0 disabled, 1 enabled). On db-joint-soerensen this removes 167184 of 2477529 columns and 2.5% of the nonzeros. Signed-off-by: yuwenchen95 Co-Authored-By: Claude Opus 5 (1M context) --- .../mathematical_optimization/constants.h | 1 + .../pdlp/solver_settings.hpp | 3 + cpp/src/dual_simplex/presolve.cpp | 185 +++++++++++++----- cpp/src/dual_simplex/presolve.hpp | 5 + .../dual_simplex/simplex_solver_settings.hpp | 2 + cpp/src/dual_simplex/solve.cpp | 14 +- cpp/src/math_optimization/solver_settings.cu | 1 + cpp/src/pdlp/solve.cu | 31 +-- 8 files changed, 181 insertions(+), 61 deletions(-) diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 64da263a4d..498ce94eb8 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -49,6 +49,7 @@ #define CUOPT_BARRIER_DUAL_INITIAL_POINT "barrier_dual_initial_point" #define CUOPT_POSTSOLVE_INFO "postsolve_info" #define CUOPT_BARRIER_PRESOLVE_BOUND_FREE_VARIABLES "barrier_presolve_bound_free_variables" +#define CUOPT_BARRIER_PRESOLVE_FIXED_VARIABLES "barrier_presolve_fixed_variables" #define CUOPT_BARRIER_ITERATIVE_REFINEMENT "barrier_iterative_refinement" #define CUOPT_BARRIER_ADAPTIVE_REGULARIZATION "barrier_adaptive_regularization" #define CUOPT_BARRIER_PRIMAL_REGULARIZATION "barrier_primal_regularization" diff --git a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp index 98f2190852..2ca5b3b11a 100644 --- a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp @@ -301,6 +301,9 @@ 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 + // Remove variables whose lower and upper bounds are equal during barrier presolve. + // -1 automatic (enabled), 0 disabled, 1 enabled + i_t barrier_presolve_fixed_variables{-1}; // 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 c204e0c798..9cc819a7a7 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -243,10 +243,9 @@ i_t remove_empty_rows(lp_problem_t& problem, template i_t remove_fixed_variables(f_t fixed_tolerance, lp_problem_t& problem, + presolve_info_t& presolve_info, i_t& fixed_variables) { - constexpr bool verbose = false; - if (verbose) { printf("Removing %d fixed variables\n", fixed_variables); } // We have a variable with l_j = x_j = u_j // Constraints of the form // @@ -258,45 +257,108 @@ i_t remove_fixed_variables(f_t fixed_tolerance, // sum_{k != j} c_k * x_k + c_j * x_j // becomes // sum_{k != j} c_k * x_k + c_j l_j + const i_t linear_cols = linear_variable_count(problem); - std::vector col_marker(problem.num_cols); - for (i_t j = 0; j < problem.num_cols; ++j) { - if (std::abs(problem.upper[j] - problem.lower[j]) < fixed_tolerance) { - col_marker[j] = 1; - for (i_t p = problem.A.col_start[j]; p < problem.A.col_start[j + 1]; ++p) { - const i_t i = problem.A.i[p]; - const f_t aij = problem.A.x[p]; - problem.rhs[i] -= aij * problem.lower[j]; - } - problem.obj_constant += problem.objective[j] * problem.lower[j]; - } else { - col_marker[j] = 0; + // A fixed variable's contribution to (1/2) x^T Q x is not a pure constant, so leave columns + // participating in the quadratic objective in place. + std::vector has_quadratic_term(problem.num_cols, false); + if (problem.Q.n > 0) { + for (i_t j = 0; j < linear_cols; ++j) { + has_quadratic_term[j] = problem.Q.row_start[j + 1] > problem.Q.row_start[j]; } } - problem.A.remove_columns(col_marker); - - // Clean up objective, lower, upper, and col_names - i_t new_cols = problem.A.n; - if (verbose) { printf("new cols %d\n", new_cols); } + // Cone columns must stay a contiguous trailing block, and the barrier rejects explicit bounds + // on them, so only linear columns are eligible. + std::vector col_marker(problem.num_cols, 0); + i_t num_removed = 0; + std::vector kahan_compensation(problem.num_rows, 0.0); + for (i_t j = 0; j < linear_cols; ++j) { + if (has_quadratic_term[j] || problem.lower[j] == -inf || problem.upper[j] == inf) { continue; } + if (std::abs(problem.upper[j] - problem.lower[j]) > fixed_tolerance) { continue; } + col_marker[j] = 1; + num_removed++; + for (i_t p = problem.A.col_start[j]; p < problem.A.col_start[j + 1]; ++p) { + const i_t i = problem.A.i[p]; + const f_t aij = problem.A.x[p]; + const f_t val = -aij * problem.lower[j]; + const f_t y = val - kahan_compensation[i]; + const f_t t = problem.rhs[i] + y; + kahan_compensation[i] = (t - problem.rhs[i]) - y; + problem.rhs[i] = t; + } + problem.obj_constant += problem.objective[j] * problem.lower[j]; + } + fixed_variables = num_removed; + if (num_removed == 0) { return 0; } + + // An earlier column reduction may already have renumbered the columns. remaining_variables maps + // the current numbering back to the original one; compose with it so the bookkeeping this + // records stays in the original index space that postsolve expects. + const bool previously_reduced = !presolve_info.remaining_variables.empty(); + auto to_original = [&](i_t j) { + return previously_reduced ? presolve_info.remaining_variables[j] : j; + }; + + const i_t new_cols = problem.num_cols - num_removed; std::vector objective(new_cols); std::vector lower(new_cols); std::vector upper(new_cols); + std::vector remaining; + remaining.reserve(new_cols); + std::vector col_old_to_new(problem.num_cols, -1); + i_t new_j = 0; for (i_t j = 0; j < problem.num_cols; ++j) { - if (!col_marker[j]) { - objective[new_j] = problem.objective[j]; - lower[new_j] = problem.lower[j]; - upper[new_j] = problem.upper[j]; - new_j++; - fixed_variables--; + if (col_marker[j]) { + // The value restored here is the current lower bound; if this variable also had its lower + // bound shifted to zero earlier, removed_lower_bounds adds the original offset back. + presolve_info.removed_variables.push_back(to_original(j)); + presolve_info.removed_values.push_back(problem.lower[j]); + presolve_info.removed_reduced_costs.push_back(0.0); + presolve_info.fixed_variables.push_back(to_original(j)); + continue; } + objective[new_j] = problem.objective[j]; + lower[new_j] = problem.lower[j]; + upper[new_j] = problem.upper[j]; + remaining.push_back(to_original(j)); + col_old_to_new[j] = new_j; + new_j++; } - problem.objective = objective; - problem.lower = lower; - problem.upper = upper; - problem.num_cols = problem.A.n; - if (verbose) { printf("Finishing fixed columns\n"); } + + problem.A.remove_columns(col_marker); + assert(new_cols == problem.A.n); + + if (problem.Q.n > 0) { + // Removed columns have no Q entries, so the rows only need compacting and reindexing. + for (i_t j = 0; j < problem.num_cols; ++j) { + const i_t mapped = col_old_to_new[j]; + if (mapped != -1) { problem.Q.row_start[mapped] = problem.Q.row_start[j]; } + } + problem.Q.row_start[new_cols] = problem.Q.row_start[problem.num_cols]; + problem.Q.row_start.resize(new_cols + 1); + for (size_t p = 0; p < problem.Q.j.size(); ++p) { + const i_t mapped = col_old_to_new[problem.Q.j[p]]; + assert(mapped != -1); + problem.Q.j[p] = mapped; + } + problem.Q.m = new_cols; + problem.Q.n = new_cols; + problem.Q.check_matrix("After removing fixed columns"); + } + + if (!problem.second_order_cone_dims.empty()) { + const i_t new_cone_start = col_old_to_new[problem.cone_var_start]; + assert(new_cone_start != -1); + problem.cone_var_start = new_cone_start; + } + + presolve_info.remaining_variables = remaining; + problem.objective = objective; + problem.lower = lower; + problem.upper = upper; + problem.num_cols = new_cols; return 0; } @@ -1363,7 +1425,30 @@ i_t presolve(const lp_problem_t& original, } } - // Check for empty rows + // Check for empty cols + i_t num_empty_cols = 0; + { + for (i_t j = 0; j < linear_cols; ++j) { + if ((problem.A.col_start[j + 1] - problem.A.col_start[j]) == 0) { num_empty_cols++; } + } + } + if (num_empty_cols > 0) { + settings.log.printf("Presolve attempt to remove %d empty cols\n", num_empty_cols); + remove_empty_cols(problem, num_empty_cols, presolve_info); + } + + // Check for fixed variables + if (settings.barrier_presolve && settings.barrier_presolve_fixed_variables != 0) { + i_t num_fixed_variables = 0; + // Exact bound equality: a tolerance here would silently perturb the primal solution. + remove_fixed_variables(static_cast(0.0), problem, presolve_info, num_fixed_variables); + if (num_fixed_variables > 0) { + settings.log.printf("Presolve removed %d fixed variables\n", num_fixed_variables); + } + } + + // Check for empty rows. This runs after the column reductions above so that rows left empty by + // them are removed too (an empty row would make A*A^T singular). i_t num_empty_rows = 0; { csr_matrix_t Arow(0, 0, 0); @@ -1378,21 +1463,13 @@ i_t presolve(const lp_problem_t& original, if (i != 0) { return -1; } } - // Check for empty cols - i_t num_empty_cols = 0; - { - for (i_t j = 0; j < linear_cols; ++j) { - if ((problem.A.col_start[j + 1] - problem.A.col_start[j]) == 0) { num_empty_cols++; } - } - } - if (num_empty_cols > 0) { - settings.log.printf("Presolve attempt to remove %d empty cols\n", num_empty_cols); - remove_empty_cols(problem, num_empty_cols, presolve_info); - } + // The column reductions above renumber the columns, so the leading linear block must be + // re-measured before it is indexed again. + const i_t reduced_linear_cols = linear_variable_count(problem); // Check for free variables (exclude cone variables — they are naturally unbounded) free_variables = 0; - for (i_t j = 0; j < linear_cols; j++) { + for (i_t j = 0; j < reduced_linear_cols; j++) { if (problem.lower[j] == -inf && problem.upper[j] == inf) { free_variables++; } } problem.Q.check_matrix("Before free variable expansion"); @@ -1406,7 +1483,7 @@ i_t presolve(const lp_problem_t& original, // Only free linear decision variables need to be handled; cone/stack columns // are unbounded by construction and must not be counted here. i_t direct_free_count = 0; - for (i_t j = 0; j < linear_cols; j++) { + for (i_t j = 0; j < reduced_linear_cols; j++) { if (problem.lower[j] == -inf && problem.upper[j] == inf) { presolve_info.direct_free_variables.push_back(j); direct_free_count++; @@ -1958,6 +2035,26 @@ void uncrush_solution(const presolve_info_t& presolve_info, } } + // Variables removed because their bounds were equal were dropped before y was known, so their + // reduced cost could not be recorded at presolve time. Recover it from the final duals as + // z_j = c_j - a_j^T y, which is what dual feasibility A^T y + z = c requires. A fixed variable's + // reduced cost is sign-unrestricted, so this is always dual feasible. Columns participating in + // the quadratic objective are never removed, so the Q x term is zero here. + if (!presolve_info.fixed_variables.empty()) { + if (settings.postsolve_info == 1) { + settings.log.printf("Post-solve: Recovering reduced costs for %d fixed variables\n", + static_cast(presolve_info.fixed_variables.size())); + } + const csc_matrix_t& A = original_problem.A; + for (const i_t j : presolve_info.fixed_variables) { + f_t zj = original_problem.objective[j]; + for (i_t p = A.col_start[j]; p < A.col_start[j + 1]; ++p) { + zj -= A.x[p] * input_y[A.i[p]]; + } + input_z[j] = zj; + } + } + assert(uncrushed_x.size() == input_x.size()); assert(uncrushed_y.size() == input_y.size()); assert(uncrushed_z.size() == input_z.size()); diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 4c0eb31b3f..e393be5ae4 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -218,6 +218,11 @@ struct presolve_info_t { // Originally-free variables that received implied bounds, with the constraint used std::vector> bounded_free_variables; + + // Original indices of variables removed because their bounds were equal. Their reduced costs + // depend on the final duals, so postsolve recomputes them rather than reading + // removed_reduced_costs. Subset of removed_variables. + std::vector fixed_variables; }; template diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 8b3eba56d3..eecdc14617 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -84,6 +84,7 @@ 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_fixed_variables(-1), qcqp_ruiz_equilibration(-1), barrier_initial_point_safeguard(10.0), check_Q(false), @@ -193,6 +194,7 @@ 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_fixed_variables; // -1 automatic (enabled), 0 disabled, 1 enabled 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/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 388bb43b35..572fa25971 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -394,8 +394,18 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us // Solve using barrier lp_solution_t barrier_solution(barrier_lp.num_rows, barrier_lp.num_cols); - barrier::barrier_solver_t barrier_solver(barrier_lp, presolve_info, barrier_settings); - lp_status_t barrier_status = barrier_solver.solve(start_time, barrier_solution); + lp_status_t barrier_status; + if (barrier_lp.num_cols == 0) { + // Presolve determined every variable, so there is no KKT system left for the barrier to form. + // What remains of the objective is the constant presolve folded the removed columns into, and + // postsolve reconstructs the variables from presolve_info. + barrier_settings.log.printf("Presolve solved the problem, skipping barrier\n"); + barrier_solution.user_objective = compute_user_objective(barrier_lp, static_cast(0.0)); + barrier_status = lp_status_t::OPTIMAL; + } else { + barrier::barrier_solver_t barrier_solver(barrier_lp, presolve_info, barrier_settings); + barrier_status = barrier_solver.solve(start_time, barrier_solution); + } if (barrier_status == lp_status_t::OPTIMAL) { #ifdef COMPUTE_SCALED_RESIDUALS std::vector scaled_residual = barrier_lp.rhs; diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 2704efd51a..4e74a0324b 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -232,6 +232,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_FIXED_VARIABLES, &pdlp_settings.barrier_presolve_fixed_variables, -1, 1, -1, "Remove variables with equal lower and upper bounds during barrier presolve: -1 automatic (enabled), 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 40b8d80378..48cf8575b7 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -512,21 +512,22 @@ 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_initial_point_safeguard = settings.barrier_initial_point_safeguard; - barrier_settings.barrier = true; - barrier_settings.barrier_presolve = true; - barrier_settings.crossover = settings.crossover; - barrier_settings.eliminate_dense_columns = settings.eliminate_dense_columns; - barrier_settings.barrier_iterative_refinement = settings.barrier_iterative_refinement; - barrier_settings.barrier_adaptive_regularization = settings.barrier_adaptive_regularization; - barrier_settings.barrier_primal_regularization = settings.barrier_primal_regularization; - barrier_settings.barrier_dual_regularization = settings.barrier_dual_regularization; - barrier_settings.barrier_soc_threshold = settings.barrier_soc_threshold; - barrier_settings.barrier_step_scale = settings.barrier_step_scale; - barrier_settings.qcqp_ruiz_equilibration = settings.qcqp_ruiz_equilibration; - barrier_settings.cudss_deterministic = settings.cudss_deterministic; - 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_presolve_fixed_variables = settings.barrier_presolve_fixed_variables; + barrier_settings.barrier_initial_point_safeguard = settings.barrier_initial_point_safeguard; + barrier_settings.barrier = true; + barrier_settings.barrier_presolve = true; + barrier_settings.crossover = settings.crossover; + barrier_settings.eliminate_dense_columns = settings.eliminate_dense_columns; + barrier_settings.barrier_iterative_refinement = settings.barrier_iterative_refinement; + barrier_settings.barrier_adaptive_regularization = settings.barrier_adaptive_regularization; + barrier_settings.barrier_primal_regularization = settings.barrier_primal_regularization; + barrier_settings.barrier_dual_regularization = settings.barrier_dual_regularization; + barrier_settings.barrier_soc_threshold = settings.barrier_soc_threshold; + barrier_settings.barrier_step_scale = settings.barrier_step_scale; + barrier_settings.qcqp_ruiz_equilibration = settings.qcqp_ruiz_equilibration; + barrier_settings.cudss_deterministic = settings.cudss_deterministic; + 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; if (barrier_settings.concurrent_halt != nullptr) { // Don't show the barrier log in concurrent mode. Show the PDLP log instead