diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 3656791a98..8eddcb9d36 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 fd6c497850..ec23582398 100644 --- a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp @@ -303,6 +303,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 61d824aeae..1ed4f37620 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -244,10 +244,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 // @@ -259,45 +258,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; } @@ -1364,7 +1426,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); @@ -1379,21 +1464,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"); @@ -1407,7 +1484,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++; @@ -1959,6 +2036,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 6f6ecc88f9..db8cb83483 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -485,51 +485,63 @@ lp_status_t solve_linear_program_with_barrier( // Solve using barrier lp_solution_t barrier_solution(barrier_lp.num_rows, barrier_lp.num_cols); - lp_problem_t const* solver_lp = &barrier_lp; - if (cache != nullptr) { - cache->clear(); - auto xf = std::make_unique(); - xf->user_num_cols = user_problem.num_cols; - xf->user_num_rows = user_problem.num_rows; - xf->original_num_cols = original_lp.num_cols; - xf->original_num_rows = original_lp.num_rows; - xf->obj_scale = user_problem.obj_scale; - xf->obj_constant = user_problem.obj_constant; - xf->row_sense = user_problem.row_sense; - xf->cone_var_start = user_problem.cone_var_start; - xf->second_order_cone_dims = user_problem.second_order_cone_dims; - xf->expanded_original_num_cols = user_problem.original_num_cols; - xf->original_col_to_expanded_col = user_problem.original_col_to_expanded_col; - xf->presolve_info = presolve_info; - xf->column_scales = column_scales; - xf->row_scales = row_scales; - xf->barrier_lp = std::make_unique>(barrier_lp); - solver_lp = xf->barrier_lp.get(); - cache->store_transform(std::move(xf)); - } - - barrier::barrier_solver_t barrier_solver(*solver_lp, presolve_info, barrier_settings); - lp_status_t barrier_status = barrier_solver.solve(start_time, barrier_solution, cache); + 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; + // No factorization was ever formed, so there is nothing a later solve could warm start from. + if (cache != nullptr) { cache->clear(); } + } else { + lp_problem_t const* solver_lp = &barrier_lp; + if (cache != nullptr) { + cache->clear(); + auto xf = std::make_unique(); + xf->user_num_cols = user_problem.num_cols; + xf->user_num_rows = user_problem.num_rows; + xf->original_num_cols = original_lp.num_cols; + xf->original_num_rows = original_lp.num_rows; + xf->obj_scale = user_problem.obj_scale; + xf->obj_constant = user_problem.obj_constant; + xf->row_sense = user_problem.row_sense; + xf->cone_var_start = user_problem.cone_var_start; + xf->second_order_cone_dims = user_problem.second_order_cone_dims; + xf->expanded_original_num_cols = user_problem.original_num_cols; + xf->original_col_to_expanded_col = user_problem.original_col_to_expanded_col; + xf->presolve_info = presolve_info; + xf->column_scales = column_scales; + xf->row_scales = row_scales; + xf->barrier_lp = std::make_unique>(barrier_lp); + solver_lp = xf->barrier_lp.get(); + cache->store_transform(std::move(xf)); + } - if (cache != nullptr) { - if (barrier_status == lp_status_t::OPTIMAL) { - auto* xf = cache->transform(); - try { - auto crushed = cuopt::mathematical_optimization::crush_user_linear_objective( - *xf, user_problem.objective.data(), user_problem.num_cols); - xf->linear_obj_shift.resize(static_cast(solver_lp->num_cols), 0.0); - if (static_cast(crushed.size()) == solver_lp->num_cols) { - for (int j = 0; j < solver_lp->num_cols; ++j) { - xf->linear_obj_shift[static_cast(j)] = - solver_lp->objective[static_cast(j)] - - crushed[static_cast(j)]; + barrier::barrier_solver_t barrier_solver(*solver_lp, presolve_info, barrier_settings); + barrier_status = barrier_solver.solve(start_time, barrier_solution, cache); + + if (cache != nullptr) { + if (barrier_status == lp_status_t::OPTIMAL) { + auto* xf = cache->transform(); + try { + auto crushed = cuopt::mathematical_optimization::crush_user_linear_objective( + *xf, user_problem.objective.data(), user_problem.num_cols); + xf->linear_obj_shift.resize(static_cast(solver_lp->num_cols), 0.0); + if (static_cast(crushed.size()) == solver_lp->num_cols) { + for (int j = 0; j < solver_lp->num_cols; ++j) { + xf->linear_obj_shift[static_cast(j)] = + solver_lp->objective[static_cast(j)] - + crushed[static_cast(j)]; + } } + } catch (std::exception const&) { + xf->linear_obj_shift.assign(static_cast(solver_lp->num_cols), 0.0); } - } catch (std::exception const&) { - xf->linear_obj_shift.assign(static_cast(solver_lp->num_cols), 0.0); + } else { + cache->clear(); } - } else { - cache->clear(); } } diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 5a4ab72c32..1b9d8e2e58 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 38202c6b51..2069b316f7 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -543,21 +543,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