From 2ec4a227c79f98e478e3e42869d19b8681ee1e30 Mon Sep 17 00:00:00 2001 From: yuwenchen95 Date: Wed, 14 Oct 2026 03:13:28 -0700 Subject: [PATCH 1/7] Add fixed-variables presolve Signed-off-by: yuwenchen95 --- .../mathematical_optimization/constants.h | 1 + .../pdlp/solver_settings.hpp | 2 + cpp/src/dual_simplex/presolve.cpp | 192 ++++++++++++++---- cpp/src/dual_simplex/presolve.hpp | 4 + .../dual_simplex/simplex_solver_settings.hpp | 2 + cpp/src/math_optimization/solver_settings.cu | 1 + cpp/src/pdlp/solve.cu | 9 + 7 files changed, 170 insertions(+), 41 deletions(-) diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 8b2c0958b9..a896d25cbc 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -50,6 +50,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 3554122963..b80f81e357 100644 --- a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp @@ -311,6 +311,8 @@ 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 + // -1 automatic (on, off when the solve feeds a barrier cache), 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 b053645e74..d3b2a60752 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -794,10 +794,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 // @@ -809,45 +808,124 @@ 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]; + // Identify columns that participate in the quadratic objective to avoid removing them as fixed + // variables. + std::vector has_quadratic_term(problem.num_cols, false); + 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) { + has_quadratic_term[row] = true; + has_quadratic_term[problem.Q.j[p]] = true; } - problem.obj_constant += problem.objective[j] * problem.lower[j]; - } else { - col_marker[j] = 0; } } - problem.A.remove_columns(col_marker); + // Map the current column index to the original problem's column index, accounting for any + // previous reductions. + const bool previously_reduced = !presolve_info.remaining_variables.empty(); + auto to_original = [&](i_t j) { + return previously_reduced ? presolve_info.remaining_variables[j] : j; + }; + // Compute the full value of a fixed variable, including any earlier lower-bound shift. + auto full_value = [&](i_t j) { + const auto& shifted = presolve_info.removed_lower_bounds; + return shifted.empty() ? problem.lower[j] : problem.lower[j] + shifted[to_original(j)]; + }; - // Clean up objective, lower, upper, and col_names - i_t new_cols = problem.A.n; - if (verbose) { printf("new cols %d\n", new_cols); } + std::vector col_marker(problem.num_cols, 0); + i_t num_removed = 0; + std::vector kahan_compensation(problem.num_rows, 0.0); + // Per row: whether it was touched by a removed fixed column, and the magnitude of the removed + // contribution. + std::vector touched(problem.num_rows, 0); + std::vector removed_magnitude(problem.num_rows, 0.0); + // Loop over all linear columns to identify and remove fixed variables. + 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++; + const f_t value = full_value(j); + 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; + touched[i] = 1; + removed_magnitude[i] += std::abs(aij * value); + } + problem.obj_constant += problem.objective[j] * problem.lower[j]; + } + fixed_variables = num_removed; + if (num_removed == 0) { return 0; } + + // Cut off rounding-level residuals of empty rows to zero. + constexpr f_t empty_row_tolerance = 1e-12; + std::vector row_nnz_left(problem.num_rows, 0); + for (i_t j = 0; j < problem.num_cols; ++j) { + if (col_marker[j]) { continue; } + for (i_t p = problem.A.col_start[j]; p < problem.A.col_start[j + 1]; ++p) { + row_nnz_left[problem.A.i[p]]++; + } + } + for (i_t i = 0; i < problem.num_rows; ++i) { + if (!touched[i] || row_nnz_left[i] != 0) { continue; } + const f_t scale = std::max(f_t(1), std::abs(problem.rhs[i]) + removed_magnitude[i]); + if (std::abs(problem.rhs[i]) <= empty_row_tolerance * scale) { problem.rhs[i] = 0.0; } + } + + 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 this only compacts and reindexes Q. + remove_variables_from_Q(problem.Q, col_marker, col_old_to_new, 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; } @@ -1828,7 +1906,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, linear_cols); + } + + // Check for fixed variables + if (settings.barrier_presolve && settings.barrier_presolve_fixed_variables != 0) { + i_t num_fixed_variables = 0; + 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); + linear_cols = linear_variable_count(problem); + } + } + + // Check for empty rows. This runs after the column reductions above so that rows left empty by + // them are removed too. i_t num_empty_rows = 0; { csr_matrix_t Arow(0, 0, 0); @@ -1843,18 +1944,6 @@ 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, linear_cols); - } - // Check for free variables (exclude cone variables — they are naturally unbounded) free_variables = 0; for (i_t j = 0; j < linear_cols; j++) { @@ -2495,6 +2584,27 @@ 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 and removed_reduced_costs holds no valid + // value for them. 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 123e87efc3..db4300b374 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -241,6 +241,10 @@ struct presolve_info_t { std::vector free_elimination_remaining_constraints; i_t free_elimination_num_variables{0}; i_t free_elimination_num_constraints{0}; + + // Original indices of variables removed because their bounds were equal. 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 d863be9804..963cf8a114 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -90,6 +90,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), gpu_ruiz_nnz_threshold(500000), barrier_initial_point_safeguard(10.0), @@ -210,6 +211,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 i_t gpu_ruiz_nnz_threshold; // nnz(A)+nnz(Q) at or above which barrier Ruiz runs on // GPU instead of CPU diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 141fb86f74..3d73691885 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -242,6 +242,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 in barrier presolve: -1 automatic (on), 0 disabled, 1 enabled; never applied when the solve feeds or reuses a barrier cache"}, {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 2d674dbf44..a92163a381 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -552,6 +552,14 @@ bool is_cached_sequence_resolve(pdlp_solver_settings_t const& settings settings.barrier_cache->transform() != nullptr; } +// Fixed variables are removed during barrier presolve. Automatic (-1) choice resolves to 0 for a +// sequence solve. +template +i_t effective_fixed_variables(pdlp_solver_settings_t const& settings) +{ + return settings.sequence_solve ? 0 : settings.barrier_presolve_fixed_variables; +} + template std::tuple, simplex::lp_status_t, f_t, f_t, f_t> run_barrier( const simplex::user_problem_t& user_problem, @@ -584,6 +592,7 @@ 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 = effective_bound_free_variables(settings); barrier_settings.barrier_eliminate_free_variables = effective_eliminate_free_variables(settings); + barrier_settings.barrier_presolve_fixed_variables = effective_fixed_variables(settings); barrier_settings.barrier_initial_point_safeguard = settings.barrier_initial_point_safeguard; barrier_settings.barrier = true; barrier_settings.barrier_presolve = true; From 2cba682162a208677eb65aa11dbd07bba5dd06e5 Mon Sep 17 00:00:00 2001 From: yuwenchen95 Date: Wed, 21 Oct 2026 03:10:21 -0700 Subject: [PATCH 2/7] Regularize the ADAT factorization for QPs Signed-off-by: yuwenchen95 --- cpp/src/barrier/barrier.cu | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index ceda102376..80efb79005 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -1252,6 +1252,9 @@ class iteration_data_t { raft::common::nvtx::range fun_scope("Barrier: Form ADAT"); float64_t start_form_adat = tic(); const i_t m = AD.m; + // QPs factor A (D + dual_perturb)^{-1} A^T + primal_perturb I, the ADAT form of the regularized + // augmented system. As there, it only preconditions IR, which solves the unregularized system. + const bool regularize = Q.n > 0; { raft::common::nvtx::range scope("Barrier: Form ADAT: restore A"); @@ -1281,6 +1284,16 @@ class iteration_data_t { d_inv_diag_prime.resize(inv_diag.size(), stream_view_); raft::copy(d_inv_diag_prime.data(), d_inv_diag.data(), inv_diag.size(), stream_view_); } + if (regularize && dual_perturb > 0) { + thrust::transform(rmm::exec_policy_nosync(stream_view_), + d_inv_diag_prime.begin(), + d_inv_diag_prime.end(), + d_inv_diag_prime.begin(), + [dual_perturb_value = dual_perturb] __device__(f_t inv_d) { + return inv_d / (f_t(1) + dual_perturb_value * inv_d); + }); + RAFT_CHECK_CUDA(stream_view_.get()); + } } cuopt_assert(static_cast(d_inv_diag_prime.size()) == AD.n, @@ -1336,6 +1349,28 @@ class iteration_data_t { multiply_kernels(handle_ptr, device_ADAT, spgemm_info()); handle_ptr->sync_stream(); } + if (regularize && primal_perturb > 0) { + raft::common::nvtx::range scope("Barrier: Form ADAT: regularize"); + // Scaling with ADAT_ii keeps the shift above pivot roundoff on dependent rows as D spreads. + constexpr f_t relative_perturb = 1e-14; + thrust::for_each_n(rmm::exec_policy_nosync(stream_view_), + thrust::make_counting_iterator(0), + m, + [row_start = cuopt::make_span(device_ADAT.row_start), + col_index = cuopt::make_span(device_ADAT.j), + values = cuopt::make_span(device_ADAT.x), + primal_perturb_value = primal_perturb, + relative_perturb] __device__(i_t i) { + for (i_t p = row_start[i]; p < row_start[i + 1]; ++p) { + if (col_index[p] == i) { + values[p] += + cuda::std::max(primal_perturb_value, relative_perturb * values[p]); + break; + } + } + }); + RAFT_CHECK_CUDA(stream_view_.get()); + } auto adat_nnz = device_ADAT.row_start.element(device_ADAT.m, handle_ptr->get_stream()); float64_t adat_time = toc(start_form_adat); From 7b08288550b7ec38986257ef15a60fd99edb655c Mon Sep 17 00:00:00 2001 From: yuwenchen95 Date: Wed, 21 Oct 2026 09:08:32 -0700 Subject: [PATCH 3/7] BARRIER_PRESOLVE_LOG prints detailed presolve info Signed-off-by: yuwenchen95 --- cpp/src/dual_simplex/presolve.cpp | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index d3b2a60752..5cabe98bff 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -1837,7 +1837,9 @@ i_t presolve(const lp_problem_t& original, if (problem.lower[j] != 0.0 && problem.lower[j] > -inf) { nonzero_lower_bounds++; } } if (settings.barrier_presolve && nonzero_lower_bounds > 0) { +#ifdef BARRIER_PRESOLVE_LOG settings.log.printf("Transforming %ld nonzero lower bound\n", nonzero_lower_bounds); +#endif presolve_info.removed_lower_bounds.resize(problem.num_cols); // We can construct a new variable: x'_j = x_j - l_j or x_j = x'_j + l_j // than we have 0 <= x'_j <= u_j - l_j @@ -1914,7 +1916,9 @@ i_t presolve(const lp_problem_t& original, } } if (num_empty_cols > 0) { +#ifdef BARRIER_PRESOLVE_LOG settings.log.printf("Presolve attempt to remove %d empty cols\n", num_empty_cols); +#endif remove_empty_cols(problem, num_empty_cols, presolve_info, linear_cols); } @@ -1923,7 +1927,9 @@ i_t presolve(const lp_problem_t& original, i_t num_fixed_variables = 0; remove_fixed_variables(static_cast(0.0), problem, presolve_info, num_fixed_variables); if (num_fixed_variables > 0) { +#ifdef BARRIER_PRESOLVE_LOG settings.log.printf("Presolve removed %d fixed variables\n", num_fixed_variables); +#endif linear_cols = linear_variable_count(problem); } } @@ -1939,7 +1945,9 @@ i_t presolve(const lp_problem_t& original, } } if (num_empty_rows > 0) { +#ifdef BARRIER_PRESOLVE_LOG settings.log.printf("Presolve removing %d empty rows\n", num_empty_rows); +#endif i_t i = remove_empty_rows(problem, row_sense, num_empty_rows, presolve_info); if (i != 0) { return -1; } } @@ -1961,15 +1969,17 @@ i_t presolve(const lp_problem_t& original, presolve_info.direct_free_variables.clear(); // 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; + [[maybe_unused]] i_t direct_free_count = 0; for (i_t j = 0; j < linear_cols; j++) { if (problem.lower[j] == -inf && problem.upper[j] == inf) { presolve_info.direct_free_variables.push_back(j); direct_free_count++; } } +#ifdef BARRIER_PRESOLVE_LOG settings.log.printf("Handling %d free variables directly in augmented system\n", direct_free_count); +#endif } else if (settings.barrier_presolve && !has_cones && free_variables > 0) { // For pure LP problems (Q is empty and there are no cones in this branch) // We have a variable x_j: with -inf < x_j < inf @@ -2057,9 +2067,11 @@ i_t presolve(const lp_problem_t& original, // LP already goes through PSLP; this substitution is for QP/SOCP only. A cache fill skips it. if (settings.barrier_presolve && settings.barrier_eliminate_free_variables && (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); + [[maybe_unused]] const i_t old_free_count = + static_cast(presolve_info.direct_free_variables.size()); + [[maybe_unused]] const f_t free_elimination_start = tic(); + [[maybe_unused]] const i_t pivot_rejected = eliminate_free_variables(problem, presolve_info); +#ifdef BARRIER_PRESOLVE_LOG const i_t eliminated = old_free_count - static_cast(presolve_info.direct_free_variables.size()); if (eliminated > 0 || pivot_rejected > 0) { @@ -2070,6 +2082,7 @@ i_t presolve(const lp_problem_t& original, pivot_rejected, toc(free_elimination_start)); } +#endif } assert(problem.num_rows == problem.A.m); From 74b48c4ff978ec4fe0f4e5388bce7f25e2791aca Mon Sep 17 00:00:00 2001 From: yuwenchen95 Date: Wed, 21 Oct 2026 03:17:55 -0700 Subject: [PATCH 4/7] Turn singleton rows into bounds or fixed values in barrier presolve Signed-off-by: yuwenchen95 --- .../mathematical_optimization/constants.h | 1 + .../pdlp/solver_settings.hpp | 2 + cpp/src/dual_simplex/presolve.cpp | 448 ++++++++++++++---- cpp/src/dual_simplex/presolve.hpp | 17 + .../dual_simplex/simplex_solver_settings.hpp | 2 + cpp/src/math_optimization/solver_settings.cu | 1 + cpp/src/pdlp/solve.cu | 9 + 7 files changed, 393 insertions(+), 87 deletions(-) diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index a896d25cbc..12009a7df9 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -51,6 +51,7 @@ #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_PRESOLVE_SINGLETON_ROWS "barrier_presolve_singleton_rows" #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 b80f81e357..29bfb1a2e4 100644 --- a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp @@ -313,6 +313,8 @@ class pdlp_solver_settings_t { i_t barrier_presolve_bound_free_variables{-1}; // -1 automatic, 0 disabled, 1 enabled // -1 automatic (on, off when the solve feeds a barrier cache), 0 disabled, 1 enabled. i_t barrier_presolve_fixed_variables{-1}; + // -1 automatic (on, off when the solve feeds a barrier cache), 0 disabled, 1 enabled. + i_t barrier_presolve_singleton_rows{-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 5cabe98bff..356a3cac7f 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -929,6 +929,279 @@ i_t remove_fixed_variables(f_t fixed_tolerance, return 0; } +/** + * @brief Shifts every finite nonzero lower bound of a linear column to zero. + * + * The shift is accumulated in presolve_info.removed_lower_bounds, indexed by original column. + */ +template +static void shift_lower_bounds_to_zero(lp_problem_t& problem, + i_t linear_cols, + i_t original_num_cols, + presolve_info_t& presolve_info) +{ + // We can construct a new variable: x'_j = x_j - l_j or x_j = x'_j + l_j + // than we have 0 <= x'_j <= u_j - l_j + // Constraints in the form: + // sum_{k != j} a_ik x_k + a_ij * x_j {=, <=} beta_i + // become + // sum_{k != j} a_ik x_k + a_ij * (x'_j + l_j) {=, <=} beta_i + // or + // sum_{k != j} a_ik x_k + a_ij * x'_j {=, <=} beta_i - a_{ij} l_j + // + // the cost function + // sum_{k != j} c_k x_k + c_j * x_j + // becomes + // sum_{k != j} c_k x_k + c_j (x'_j + l_j) + // + // so we get the constant term c_j * l_j + + std::vector lower_bounds_removed(problem.num_cols, false); + bool any_removed = false; + for (i_t j = 0; j < linear_cols; j++) { + if (problem.lower[j] != 0.0 && problem.lower[j] > -inf) { + lower_bounds_removed[j] = true; + any_removed = true; + } + } + if (!any_removed) { return; } + + if (presolve_info.removed_lower_bounds.empty()) { + presolve_info.removed_lower_bounds.resize(original_num_cols, 0.0); + } + const bool previously_reduced = !presolve_info.remaining_variables.empty(); + for (i_t j = 0; j < linear_cols; j++) { + if (!lower_bounds_removed[j]) { continue; } + const i_t original_j = previously_reduced ? presolve_info.remaining_variables[j] : j; + presolve_info.removed_lower_bounds[original_j] += problem.lower[j]; + } + + auto old_objective = problem.objective; + if (problem.Q.n > 0) { + for (i_t row = 0; row < linear_cols; row++) { + i_t row_start = problem.Q.row_start[row]; + i_t row_end = problem.Q.row_start[row + 1]; + for (i_t p = row_start; p < row_end; p++) { + i_t col = problem.Q.j[p]; + f_t qij = problem.Q.x[p]; + + if (lower_bounds_removed[row]) { problem.objective[col] += 0.5 * qij * problem.lower[row]; } + if (lower_bounds_removed[col]) { problem.objective[row] += 0.5 * qij * problem.lower[col]; } + if (lower_bounds_removed[row] && lower_bounds_removed[col]) { + problem.obj_constant += 0.5 * qij * problem.lower[row] * problem.lower[col]; + } + } + } + } + + std::vector kahan_compensation(problem.num_rows, 0.0); + for (i_t j = 0; j < linear_cols; j++) { + if (lower_bounds_removed[j]) { + i_t col_start = problem.A.col_start[j]; + i_t col_end = problem.A.col_start[j + 1]; + for (i_t p = col_start; p < col_end; p++) { + i_t i = problem.A.i[p]; + f_t aij = problem.A.x[p]; + f_t val = -aij * problem.lower[j]; + f_t y = val - kahan_compensation[i]; + f_t t = problem.rhs[i] + y; + kahan_compensation[i] = (t - problem.rhs[i]) - y; + problem.rhs[i] = t; + } + problem.obj_constant += old_objective[j] * problem.lower[j]; + problem.upper[j] -= problem.lower[j]; + problem.lower[j] = 0.0; + } + } +} + +/** + * @brief Turns rows that fix or bound a single variable into bounds and empties those rows. + * + * A row qualifies if it has one entry, or two entries one of which is a zero-cost column singleton + * with at least one finite bound (typically an inequality slack). The emptied rows are left for + * remove_empty_rows and their slacks are fixed for remove_fixed_variables. + * + * @return The number of rows emptied. + */ +template +static i_t remove_singleton_rows(lp_problem_t& problem, + i_t original_num_cols, + presolve_info_t& presolve_info) +{ + const i_t m = problem.num_rows; + const i_t n = problem.num_cols; + const i_t linear_cols = linear_variable_count(problem); + + std::vector has_quadratic_term(n, 0); + if (problem.Q.n > 0) { + for (i_t row = 0; row < std::min(problem.Q.m, n); ++row) { + for (i_t p = problem.Q.row_start[row]; p < problem.Q.row_start[row + 1]; ++p) { + has_quadratic_term[row] = 1; + has_quadratic_term[problem.Q.j[p]] = 1; + } + } + } + + const bool previously_reduced = !presolve_info.remaining_variables.empty(); + auto to_original = [&](i_t j) { + return previously_reduced ? presolve_info.remaining_variables[j] : j; + }; + std::vector negated(original_num_cols, 0); + for (const i_t j : presolve_info.negated_variables) { + negated[j] = 1; + } + + // A bound within feas_tol of the opposite bound is snapped to it; further away the row is kept. + constexpr f_t feas_tol = 1e-9; + auto close = [&](f_t a, f_t b) { + return std::abs(a - b) <= feas_tol * std::max({f_t(1), std::abs(a), std::abs(b)}); + }; + auto is_slack = [&](i_t j) { + return j < linear_cols && problem.A.col_length(j) == 1 && problem.objective[j] == 0 && + !has_quadratic_term[j] && (problem.lower[j] > -inf || problem.upper[j] < inf); + }; + + csr_matrix_t Arow(0, 0, 0); + problem.A.to_compressed_row(Arow); + std::vector remove_row(m, 0); + // Columns fixed by a row in this pass, which other rows must not bound again. + std::vector fixed_here(n, 0); + i_t removed = 0; + + for (i_t i = 0; i < m; ++i) { + const i_t row_start = Arow.row_start[i]; + const i_t row_len = Arow.row_start[i + 1] - row_start; + if (row_len == 0 || row_len > 2) { continue; } + bool touches_cone = false; + for (i_t p = row_start; p < row_start + row_len; ++p) { + if (Arow.j[p] >= linear_cols) { touches_cone = true; } + } + if (touches_cone) { continue; } + + singleton_row_t record{i, -1, -1, false, false}; + if (row_len == 1) { + // a_ij x_j = b_i fixes x_j. + const i_t j = Arow.j[row_start]; + if (has_quadratic_term[j] || Arow.x[row_start] == 0) { continue; } + f_t value = problem.rhs[i] / Arow.x[row_start]; + if (fixed_here[j]) { + // A second row fixing the same value is redundant and keeps a zero dual. + if (!close(value, problem.lower[j])) { continue; } + remove_row[i] = 1; + ++removed; + continue; + } + if (value < problem.lower[j]) { + if (!close(value, problem.lower[j])) { continue; } + value = problem.lower[j]; + } + if (value > problem.upper[j]) { + if (!close(value, problem.upper[j])) { continue; } + value = problem.upper[j]; + } + problem.lower[j] = value; + problem.upper[j] = value; + fixed_here[j] = 1; + record.variable = to_original(j); + record.sets_lower = true; + record.sets_upper = true; + } else { + // a_ij x_j + a_is s = b_i with s in [l_s, u_s] bounds a_ij x_j to [b_i - a_is s] over s. + i_t s = -1; + i_t j = -1; + for (i_t k = 1; k >= 0; --k) { + const i_t cand_s = Arow.j[row_start + k]; + const i_t cand_j = Arow.j[row_start + 1 - k]; + if (is_slack(cand_s) && problem.lower[cand_j] > -inf && !fixed_here[cand_j] && + !fixed_here[cand_s]) { + s = cand_s; + j = cand_j; + break; + } + } + if (s == -1) { continue; } + const f_t a_ij = s == Arow.j[row_start] ? Arow.x[row_start + 1] : Arow.x[row_start]; + const f_t a_is = s == Arow.j[row_start] ? Arow.x[row_start] : Arow.x[row_start + 1]; + if (a_ij == 0 || a_is == 0) { continue; } + const f_t b = problem.rhs[i]; + // Range of a_is * s, then of a_ij * x_j. + const f_t slack_lo = a_is > 0 ? a_is * problem.lower[s] : a_is * problem.upper[s]; + const f_t slack_hi = a_is > 0 ? a_is * problem.upper[s] : a_is * problem.lower[s]; + const f_t activity_lo = b - slack_hi; + const f_t activity_hi = b - slack_lo; + f_t new_lower = a_ij > 0 ? activity_lo / a_ij : activity_hi / a_ij; + f_t new_upper = a_ij > 0 ? activity_hi / a_ij : activity_lo / a_ij; + if (std::isnan(new_lower)) { new_lower = -inf; } + if (std::isnan(new_upper)) { new_upper = inf; } + + f_t lower = problem.lower[j]; + f_t upper = problem.upper[j]; + if (new_lower > lower) { + lower = new_lower; + record.sets_lower = true; + } + if (new_upper < upper) { + upper = new_upper; + record.sets_upper = true; + } + if (lower > upper) { + if (!close(lower, upper)) { continue; } + // Keep the bound that was already there, so the row only fixes the variable. + if (record.sets_lower) { + lower = upper; + } else { + upper = lower; + } + } + // Fixed columns with a quadratic term stay in the barrier with no interior, so leave them. + if (has_quadratic_term[j] && lower == upper) { continue; } + problem.lower[j] = lower; + problem.upper[j] = upper; + if (lower == upper) { fixed_here[j] = 1; } + + // The slack's value is recovered from the row in postsolve, so fix it anywhere. + problem.lower[s] = 0.0; + problem.upper[s] = 0.0; + fixed_here[s] = 1; + record.variable = to_original(j); + record.slack = to_original(s); + } + // Bounds are recorded in the original orientation of the variable. + if (negated[record.variable]) { std::swap(record.sets_lower, record.sets_upper); } + if (record.slack != -1 || record.sets_lower || record.sets_upper) { + presolve_info.singleton_rows.push_back(record); + } + remove_row[i] = 1; + ++removed; + } + if (removed == 0) { return 0; } + + // Empty the removed rows. Their entries no longer constrain anything once the bounds are moved. + i_t nz = 0; + for (i_t j = 0; j < n; ++j) { + const i_t col_start = problem.A.col_start[j]; + const i_t col_end = problem.A.col_start[j + 1]; + problem.A.col_start[j] = nz; + for (i_t p = col_start; p < col_end; ++p) { + if (remove_row[problem.A.i[p]]) { continue; } + problem.A.i[nz] = problem.A.i[p]; + problem.A.x[nz] = problem.A.x[p]; + ++nz; + } + } + problem.A.col_start[n] = nz; + problem.A.i.resize(nz); + problem.A.x.resize(nz); + problem.A.nz_max = nz; + for (i_t i = 0; i < m; ++i) { + if (remove_row[i]) { problem.rhs[i] = 0.0; } + } + + shift_lower_bounds_to_zero(problem, linear_cols, original_num_cols, presolve_info); + return removed; +} + // Adds one slack column s_k per entry of `rows`; slack k has a single nonzero in row i = rows[k]: // "<=" row (is_range = false): a_i^T x + s_k = rhs_i, s_k >= 0 // range row (is_range = true): a_i^T x - s_k = 0, lower[k] <= s_k <= upper[k] @@ -1840,72 +2113,7 @@ i_t presolve(const lp_problem_t& original, #ifdef BARRIER_PRESOLVE_LOG settings.log.printf("Transforming %ld nonzero lower bound\n", nonzero_lower_bounds); #endif - presolve_info.removed_lower_bounds.resize(problem.num_cols); - // We can construct a new variable: x'_j = x_j - l_j or x_j = x'_j + l_j - // than we have 0 <= x'_j <= u_j - l_j - // Constraints in the form: - // sum_{k != j} a_ik x_k + a_ij * x_j {=, <=} beta_i - // become - // sum_{k != j} a_ik x_k + a_ij * (x'_j + l_j) {=, <=} beta_i - // or - // sum_{k != j} a_ik x_k + a_ij * x'_j {=, <=} beta_i - a_{ij} l_j - // - // the cost function - // sum_{k != j} c_k x_k + c_j * x_j - // becomes - // sum_{k != j} c_k x_k + c_j (x'_j + l_j) - // - // so we get the constant term c_j * l_j - - std::vector lower_bounds_removed(problem.num_cols, false); - for (i_t j = 0; j < linear_cols; j++) { - if (problem.lower[j] != 0.0 && problem.lower[j] > -inf) { - lower_bounds_removed[j] = true; - presolve_info.removed_lower_bounds[j] = problem.lower[j]; - } - } - - auto old_objective = problem.objective; - if (problem.Q.n > 0) { - for (i_t row = 0; row < linear_cols; row++) { - i_t row_start = problem.Q.row_start[row]; - i_t row_end = problem.Q.row_start[row + 1]; - for (i_t p = row_start; p < row_end; p++) { - i_t col = problem.Q.j[p]; - f_t qij = problem.Q.x[p]; - - if (lower_bounds_removed[row]) { - problem.objective[col] += 0.5 * qij * problem.lower[row]; - } - if (lower_bounds_removed[col]) { - problem.objective[row] += 0.5 * qij * problem.lower[col]; - } - if (lower_bounds_removed[row] && lower_bounds_removed[col]) { - problem.obj_constant += 0.5 * qij * problem.lower[row] * problem.lower[col]; - } - } - } - } - - std::vector kahan_compensation(problem.num_rows, 0.0); - for (i_t j = 0; j < linear_cols; j++) { - if (lower_bounds_removed[j]) { - i_t col_start = problem.A.col_start[j]; - i_t col_end = problem.A.col_start[j + 1]; - for (i_t p = col_start; p < col_end; p++) { - i_t i = problem.A.i[p]; - f_t aij = problem.A.x[p]; - f_t val = -aij * problem.lower[j]; - f_t y = val - kahan_compensation[i]; - f_t t = problem.rhs[i] + y; - kahan_compensation[i] = (t - problem.rhs[i]) - y; - problem.rhs[i] = t; - } - problem.obj_constant += old_objective[j] * problem.lower[j]; - problem.upper[j] -= problem.lower[j]; - problem.lower[j] = 0.0; - } - } + shift_lower_bounds_to_zero(problem, linear_cols, original.num_cols, presolve_info); } // Check for empty cols @@ -1932,6 +2140,24 @@ i_t presolve(const lp_problem_t& original, #endif linear_cols = linear_variable_count(problem); } + + // Fixing a variable can leave further rows with a single variable, so repeat until none do. + i_t num_singleton_rows = 0; + i_t num_fixed_by_rows = 0; + while (settings.barrier_presolve_singleton_rows != 0) { + const i_t rows = remove_singleton_rows(problem, original.num_cols, presolve_info); + if (rows == 0) { break; } + num_singleton_rows += rows; + i_t fixed = 0; + remove_fixed_variables(static_cast(0.0), problem, presolve_info, fixed); + num_fixed_by_rows += fixed; + } + if (num_singleton_rows > 0) { + settings.log.printf("Presolve removed %d singleton rows and %d variables they fixed\n", + num_singleton_rows, + num_fixed_by_rows); + linear_cols = linear_variable_count(problem); + } } // Check for empty rows. This runs after the column reductions above so that rows left empty by @@ -2557,6 +2783,75 @@ void uncrush_solution(const presolve_info_t& presolve_info, } } + csr_matrix_t original_Arow(0, 0, 0); + if (!presolve_info.singleton_rows.empty()) { + original_problem.A.to_compressed_row(original_Arow); + // A slack dropped with its singleton row takes the value that satisfies the row. + for (const auto& record : presolve_info.singleton_rows) { + if (record.slack == -1) { continue; } + f_t activity = 0; + f_t a_is = 0; + for (i_t p = original_Arow.row_start[record.row]; p < original_Arow.row_start[record.row + 1]; + ++p) { + if (original_Arow.j[p] == record.slack) { + a_is = original_Arow.x[p]; + } else { + activity += original_Arow.x[p] * input_x[original_Arow.j[p]]; + } + } + input_x[record.slack] = (original_problem.rhs[record.row] - activity) / a_is; + } + } + + // Variables removed because their bounds were equal were dropped before y was known, so their + // reduced cost could not be recorded at presolve time and removed_reduced_costs holds no valid + // value for them. 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; + } + } + + // A singleton row that supplied the active bound of its variable takes over that bound's + // multiplier: with du = z_j / a_ij, y_i += du and z -= du a_i keep A^T y + z = c + Qx and + // leave z_j = 0. Later rows are undone first, since they may have tightened an earlier bound. + if (!presolve_info.singleton_rows.empty()) { + if (settings.postsolve_info == 1) { + settings.log.printf("Post-solve: Recovering duals for %d singleton rows\n", + static_cast(presolve_info.singleton_rows.size())); + } + for (auto it = presolve_info.singleton_rows.rbegin(); it != presolve_info.singleton_rows.rend(); + ++it) { + const auto& record = *it; + const f_t z_j = input_z[record.variable]; + if (!((z_j > 0 && record.sets_lower) || (z_j < 0 && record.sets_upper))) { continue; } + const i_t row_start = original_Arow.row_start[record.row]; + const i_t row_end = original_Arow.row_start[record.row + 1]; + f_t a_ij = 0; + for (i_t p = row_start; p < row_end; ++p) { + if (original_Arow.j[p] == record.variable) { a_ij = original_Arow.x[p]; } + } + if (a_ij == 0) { continue; } + const f_t du = z_j / a_ij; + input_y[record.row] += du; + for (i_t p = row_start; p < row_end; ++p) { + input_z[original_Arow.j[p]] -= original_Arow.x[p] * du; + } + } + } + // Dual correction for originally free variables that received implied bounds. // Barrier produced (y, z) with z_j != 0 satisfying A^T y + z = c + Qx. // We need corrected (y_bar, z_bar) with z_bar_j = 0 for all j in F_b where @@ -2597,27 +2892,6 @@ 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 and removed_reduced_costs holds no valid - // value for them. 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 db4300b374..c539549b12 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -203,6 +203,20 @@ struct free_variable_elimination_t { std::vector factors; }; +// Row removed because it fixed or bounded a single variable. Indices refer to the problem before +// presolve. +template +struct singleton_row_t { + i_t row; + // The variable fixed or bounded by the row. + i_t variable; + // Zero-cost column singleton dropped with the row, or -1. + i_t slack; + // Whether the row supplied the variable's lower or upper bound. + bool sets_lower; + bool sets_upper; +}; + template struct presolve_info_t { // indices of variables in the original problem that remain in the presolved problem @@ -245,6 +259,9 @@ struct presolve_info_t { // Original indices of variables removed because their bounds were equal. Subset of // removed_variables. std::vector fixed_variables; + + // Rows removed because they fixed or bounded a single variable, in removal order. + std::vector> singleton_rows; }; template diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 963cf8a114..a710b8f2e0 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -91,6 +91,7 @@ struct simplex_solver_settings_t { postsolve_info(-1), barrier_presolve_bound_free_variables(-1), barrier_presolve_fixed_variables(-1), + barrier_presolve_singleton_rows(-1), qcqp_ruiz_equilibration(-1), gpu_ruiz_nnz_threshold(500000), barrier_initial_point_safeguard(10.0), @@ -212,6 +213,7 @@ struct simplex_solver_settings_t { 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 barrier_presolve_singleton_rows; // -1 automatic (enabled), 0 disabled, 1 enabled i_t qcqp_ruiz_equilibration; // -1 automatic (imbalance heuristic), 0 disabled, 1 enabled i_t gpu_ruiz_nnz_threshold; // nnz(A)+nnz(Q) at or above which barrier Ruiz runs on // GPU instead of CPU diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 3d73691885..d3d391ba78 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -243,6 +243,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {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 in barrier presolve: -1 automatic (on), 0 disabled, 1 enabled; never applied when the solve feeds or reuses a barrier cache"}, + {CUOPT_BARRIER_PRESOLVE_SINGLETON_ROWS, &pdlp_settings.barrier_presolve_singleton_rows, -1, 1, -1, "Turn rows with a single variable into bounds in barrier presolve (requires barrier_presolve_fixed_variables): -1 automatic (on), 0 disabled, 1 enabled; never applied when the solve feeds or reuses a barrier cache"}, {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 a92163a381..bb5bc2c6ca 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -560,6 +560,14 @@ i_t effective_fixed_variables(pdlp_solver_settings_t const& settings) return settings.sequence_solve ? 0 : settings.barrier_presolve_fixed_variables; } +// Singleton rows are removed during barrier presolve. Automatic (-1) choice resolves to 0 for a +// sequence solve. +template +i_t effective_singleton_rows(pdlp_solver_settings_t const& settings) +{ + return settings.sequence_solve ? 0 : settings.barrier_presolve_singleton_rows; +} + template std::tuple, simplex::lp_status_t, f_t, f_t, f_t> run_barrier( const simplex::user_problem_t& user_problem, @@ -593,6 +601,7 @@ std::tuple, simplex::lp_status_t, f_t, f_t, f_t barrier_settings.barrier_presolve_bound_free_variables = effective_bound_free_variables(settings); barrier_settings.barrier_eliminate_free_variables = effective_eliminate_free_variables(settings); barrier_settings.barrier_presolve_fixed_variables = effective_fixed_variables(settings); + barrier_settings.barrier_presolve_singleton_rows = effective_singleton_rows(settings); barrier_settings.barrier_initial_point_safeguard = settings.barrier_initial_point_safeguard; barrier_settings.barrier = true; barrier_settings.barrier_presolve = true; From f7cb448b5011ffceeb06e74fb8ef1f380ae65fca Mon Sep 17 00:00:00 2001 From: yuwenchen95 Date: Thu, 22 Oct 2026 05:01:39 -0700 Subject: [PATCH 5/7] fix the uncrush bug: uncrush in reverse order Signed-off-by: yuwenchen95 --- cpp/src/dual_simplex/presolve.cpp | 4 +- cpp/tests/socp/solve_barrier_socp.cu | 99 ++++++++++++++++++++++++++++ 2 files changed, 102 insertions(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 356a3cac7f..aecc5f433a 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -2787,7 +2787,9 @@ void uncrush_solution(const presolve_info_t& presolve_info, if (!presolve_info.singleton_rows.empty()) { original_problem.A.to_compressed_row(original_Arow); // A slack dropped with its singleton row takes the value that satisfies the row. - for (const auto& record : presolve_info.singleton_rows) { + for (auto it = presolve_info.singleton_rows.rbegin(); it != presolve_info.singleton_rows.rend(); + ++it) { + const auto& record = *it; if (record.slack == -1) { continue; } f_t activity = 0; f_t a_is = 0; diff --git a/cpp/tests/socp/solve_barrier_socp.cu b/cpp/tests/socp/solve_barrier_socp.cu index 8a9d42bea5..bebcc8016d 100644 --- a/cpp/tests/socp/solve_barrier_socp.cu +++ b/cpp/tests/socp/solve_barrier_socp.cu @@ -1387,6 +1387,105 @@ TEST(barrier, free_variable_substitution_postsolve_kkt) EXPECT_NEAR((vector_norm_inf(solved_dual)), 0.0, 1e-5); } +TEST(barrier, chained_singleton_row_slacks_are_uncrushed_in_reverse_order) +{ + // minimize x0 s.t. x0 + x1 == 3, x1 <= 5, 0 <= x0, x1 <= 10, with x1 zero-cost. The optimum + // is x0 = 0, x1 = 3, so the <= row's slack is s = 2. + // + // Presolve empties both rows in two passes: pass 1 takes the <= row, bounding x1 and + // recording its slack s. That leaves x1 a zero-cost column singleton, so pass 2 takes the + // equality and records x1 as its own slack. Postsolve must replay the records in reverse, + // because s is computed from x1 and only the later record restores x1. + raft::handle_t handle{}; + init_handler(&handle); + + constexpr int m = 2; + constexpr int n = 2; + constexpr int nz = 3; + + user_problem_t user_problem(&handle); + user_problem.num_rows = m; + user_problem.num_cols = n; + // x1 carries no cost, which is what lets it stand in as a slack in the second pass. + user_problem.objective = {1.0, 0.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + // Columns: x0 (r0 only), x1 (r0 and r1). + user_problem.A.col_start = {0, 1, 3}; + user_problem.A.i = {0, 0, 1}; + user_problem.A.x = {1.0, 1.0, 1.0}; + + user_problem.rhs = {3.0, 5.0}; + user_problem.row_sense = {'E', 'L'}; + user_problem.lower = {0.0, 0.0}; + user_problem.upper = {10.0, 10.0}; + + user_problem.num_range_rows = 0; + user_problem.problem_name = "chained_singleton_rows"; + user_problem.cone_var_start = n; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.barrier_presolve_fixed_variables = 1; + settings.barrier_presolve_singleton_rows = 1; + 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); + // x0, x1 and the slack of r1. + ASSERT_EQ(original_lp.num_cols, 3); + ASSERT_EQ(original_lp.num_rows, 2); + + 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); + + // One record per removed row, and both carry a slack to recover -- that pairing is the chain. + ASSERT_EQ(presolve_info.singleton_rows.size(), 2u); + EXPECT_NE(presolve_info.singleton_rows[0].slack, -1); + EXPECT_NE(presolve_info.singleton_rows[1].slack, -1); + // The later record restores the variable the earlier record's slack is computed from. + EXPECT_EQ(presolve_info.singleton_rows[1].slack, presolve_info.singleton_rows[0].variable); + + // Both rows are gone, so only x0 survives and the crushed optimum is x0 = 0. + ASSERT_EQ(presolved_lp.num_rows, 0); + ASSERT_EQ(presolved_lp.num_cols, 1); + std::vector crushed_x{0.0}; + std::vector crushed_y(presolved_lp.num_rows, 0.0); + std::vector crushed_z{1.0}; + + 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, + crushed_x, + crushed_y, + crushed_z, + uncrushed_x, + uncrushed_y, + uncrushed_z); + + EXPECT_NEAR(uncrushed_x[0], 0.0, 1e-12); + EXPECT_NEAR(uncrushed_x[1], 3.0, 1e-12); + // s = 5 - x1. Recovered forward it comes out 5, because x1 is still its placeholder. + EXPECT_NEAR(uncrushed_x[2], 2.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((vector_norm_inf(primal_residual)), 0.0, 1e-12); +} + TEST(barrier, cached_solve_can_enable_iterative_refinement) { // The first solve has no GMRES workspace. Enabling refinement on a cached solve must From 383af0063b0d8846a4812741365de938e606d96c Mon Sep 17 00:00:00 2001 From: yuwenchen95 Date: Thu, 22 Oct 2026 06:08:58 -0700 Subject: [PATCH 6/7] Update regularization in ADAT path Signed-off-by: yuwenchen95 --- .../mathematical_optimization/constants.h | 1 + .../pdlp/solver_settings.hpp | 2 ++ cpp/src/barrier/barrier.cu | 18 +++++++----------- .../dual_simplex/simplex_solver_settings.hpp | 3 +++ cpp/src/math_optimization/solver_settings.cu | 1 + cpp/src/pdlp/solve.cu | 2 ++ 6 files changed, 16 insertions(+), 11 deletions(-) diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 12009a7df9..8567c92f23 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -56,6 +56,7 @@ #define CUOPT_BARRIER_ADAPTIVE_REGULARIZATION "barrier_adaptive_regularization" #define CUOPT_BARRIER_PRIMAL_REGULARIZATION "barrier_primal_regularization" #define CUOPT_BARRIER_DUAL_REGULARIZATION "barrier_dual_regularization" +#define CUOPT_BARRIER_ADAT_REGULARIZATION "barrier_adat_regularization" #define CUOPT_BARRIER_STEP_SCALE "barrier_step_scale" #define CUOPT_ELIMINATE_DENSE_COLUMNS "eliminate_dense_columns" #define CUOPT_CUDSS_DETERMINISTIC "cudss_deterministic" diff --git a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp index 29bfb1a2e4..9cc163c619 100644 --- a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp @@ -334,6 +334,8 @@ class pdlp_solver_settings_t { // iterations). -1 automatic (uses the built-in heuristic), else the literal starting value. f_t barrier_primal_regularization{-1.0}; f_t barrier_dual_regularization{-1.0}; + // Relative diagonal shift of the QP ADAT factorization (ADAT_ii scaled by 1 + value); 0 disables. + f_t barrier_adat_regularization{1e-14}; i_t barrier_soc_threshold{100}; f_t barrier_step_scale{0.9}; bool save_best_primal_so_far{false}; diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index 80efb79005..0ec800e7d9 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -1252,8 +1252,7 @@ class iteration_data_t { raft::common::nvtx::range fun_scope("Barrier: Form ADAT"); float64_t start_form_adat = tic(); const i_t m = AD.m; - // QPs factor A (D + dual_perturb)^{-1} A^T + primal_perturb I, the ADAT form of the regularized - // augmented system. As there, it only preconditions IR, which solves the unregularized system. + // QPs factor A (D + dual_perturb)^{-1} A^T with a small relative diagonal shift. const bool regularize = Q.n > 0; { @@ -1349,22 +1348,19 @@ class iteration_data_t { multiply_kernels(handle_ptr, device_ADAT, spgemm_info()); handle_ptr->sync_stream(); } - if (regularize && primal_perturb > 0) { + if (regularize && settings_.barrier_adat_regularization > 0) { raft::common::nvtx::range scope("Barrier: Form ADAT: regularize"); - // Scaling with ADAT_ii keeps the shift above pivot roundoff on dependent rows as D spreads. - constexpr f_t relative_perturb = 1e-14; + const f_t relative_perturb = settings_.barrier_adat_regularization; thrust::for_each_n(rmm::exec_policy_nosync(stream_view_), thrust::make_counting_iterator(0), m, - [row_start = cuopt::make_span(device_ADAT.row_start), - col_index = cuopt::make_span(device_ADAT.j), - values = cuopt::make_span(device_ADAT.x), - primal_perturb_value = primal_perturb, + [row_start = cuopt::make_span(device_ADAT.row_start), + col_index = cuopt::make_span(device_ADAT.j), + values = cuopt::make_span(device_ADAT.x), relative_perturb] __device__(i_t i) { for (i_t p = row_start[i]; p < row_start[i + 1]; ++p) { if (col_index[p] == i) { - values[p] += - cuda::std::max(primal_perturb_value, relative_perturb * values[p]); + values[p] += relative_perturb * values[p]; break; } } diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index a710b8f2e0..5125dc38ef 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -77,6 +77,7 @@ struct simplex_solver_settings_t { barrier_adaptive_regularization(-1), barrier_primal_regularization(-1.0), barrier_dual_regularization(-1.0), + barrier_adat_regularization(1e-14), barrier_step_scale(0.9), barrier_soc_threshold(100), num_gpus(1), @@ -197,6 +198,8 @@ struct simplex_solver_settings_t { // initial dual regularization (augmented system's (1,1) block // diagonal) for the first barrier factorization. Same adaptive // caveat as above. + f_t barrier_adat_regularization; // relative diagonal shift of the QP ADAT factorization, 0 + // disables f_t barrier_step_scale; // step scale for barrier method i_t barrier_soc_threshold; // SOC dimension above which rank-2 sparse scaling is used int num_gpus; // Number of GPUs to use (maximum of 2 gpus are supported at the moment) diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index d3d391ba78..a84dced1b1 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -143,6 +143,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_CUT_MIN_ORTHOGONALITY, &mip_settings.cut_min_orthogonality, f_t(0.0), f_t(1.0), f_t(0.5)}, {CUOPT_BARRIER_PRIMAL_REGULARIZATION, &pdlp_settings.barrier_primal_regularization, f_t(-1.0), std::numeric_limits::infinity(), f_t(-1.0), "initial primal regularization for the augmented system; -1 automatic"}, {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_ADAT_REGULARIZATION, &pdlp_settings.barrier_adat_regularization, f_t(0.0), std::numeric_limits::infinity(), f_t(1e-14), "relative diagonal shift of the ADAT factorization for QPs; 0 disables"}, {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"}, // MIP heuristic hyper-parameters (hidden from default --help: name contains "hyper_") diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index bb5bc2c6ca..ac9f2b3952 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -611,6 +611,7 @@ std::tuple, simplex::lp_status_t, f_t, f_t, f_t 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_adat_regularization = settings.barrier_adat_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; @@ -865,6 +866,7 @@ static optimization_problem_solution_t run_pdlp_solver_in_fp32( fs.barrier_adaptive_regularization = settings.barrier_adaptive_regularization; fs.barrier_primal_regularization = settings.barrier_primal_regularization; fs.barrier_dual_regularization = settings.barrier_dual_regularization; + fs.barrier_adat_regularization = settings.barrier_adat_regularization; fs.barrier_step_scale = settings.barrier_step_scale; fs.pdlp_precision = pdlp_precision_t::DefaultPrecision; fs.method = method_t::PDLP; From 6ba91aff762da6df83c24d1600edfe25a72a5a8c Mon Sep 17 00:00:00 2001 From: yuwenchen95 Date: Thu, 22 Oct 2026 08:11:52 -0700 Subject: [PATCH 7/7] Fix CI failure issues Signed-off-by: yuwenchen95 --- .../linear_programming/test_barrier_sequence_solve.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/python/cuopt/cuopt/tests/linear_programming/test_barrier_sequence_solve.py b/python/cuopt/cuopt/tests/linear_programming/test_barrier_sequence_solve.py index bad2f1c4b8..e042a238a6 100644 --- a/python/cuopt/cuopt/tests/linear_programming/test_barrier_sequence_solve.py +++ b/python/cuopt/cuopt/tests/linear_programming/test_barrier_sequence_solve.py @@ -179,7 +179,9 @@ def test_update_rhs_keeps_dropped_empty_row(capfd): first, log = _solve(model, settings, capfd) assert first.get_termination_reason() == "Optimal" - assert "empty rows" in log, "presolve did not drop the empty row" + assert "Presolve eliminated 1 constraints" in log, ( + "presolve did not drop the empty row" + ) rhs = [0.0, 4.0] model.update_rhs(np.asarray(rhs, dtype=np.float64)) @@ -329,7 +331,8 @@ def test_lp_barrier_update_rhs_reuses_cache(capfd): assert first.get_termination_reason() == "Optimal" assert REUSE_LOG not in log assert "Using PSLP presolver" not in log - assert "Handling 1 free variables directly" in log + # The free variable reaches the barrier unsplit. + assert "Free variables : 1\n" in log assert "Folding:" not in log rhs = [6.0]