From 278d02a191c4466bae3c09efcc708f3e4365c4eb Mon Sep 17 00:00:00 2001 From: yuwenchen95 Date: Tue, 14 Jul 2026 07:42:07 -0700 Subject: [PATCH 1/4] Extend Sedumi initialization to LP/QP Signed-off-by: yuwenchen95 --- cpp/src/barrier/barrier.cu | 56 +++++++++++-------- .../dual_simplex/simplex_solver_settings.hpp | 3 +- cpp/src/math_optimization/solver_settings.cu | 2 +- .../cuopt_server/tests/test_lp.py | 2 +- .../linear_programming/data_definition.py | 3 +- 5 files changed, 38 insertions(+), 28 deletions(-) diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index 784f6c0901..e55abddce3 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -2155,41 +2155,49 @@ int barrier_solver_t::initial_point(iteration_data_t& data) const bool use_augmented = data.use_augmented; const bool has_direct_free_linear = data.n_direct_free_linear > 0; - // SOCP: data-dependent initial point following SeDuMi (Sturm, 1999). - // mu = sqrt((1 + ||b||_inf) * (1 + ||c||_inf)) - // primal and dual: x = mu * e_K, z = mu * e_K + const i_t init_strategy = data.has_cones() ? 2 : settings.barrier_dual_initial_point; + + // Option 2: Sturm/SeDuMi-style mu-based primal+dual initial point. + // mu = sqrt((1 + ||b||_inf) * (1 + ||c||_inf)); x = z = mu * e_K. // where e_K is the identity of the symmetric cone: // LP block: e = 1, SOC block: e = (sqrt(2), 0, ..., 0) - if (data.has_cones()) { - const i_t cs = data.cone_start(); - const f_t norm_b = vector_norm_inf(lp.rhs); - const f_t norm_c = vector_norm_inf(lp.objective); - const f_t mu = std::sqrt((1.0 + norm_b) * (1.0 + norm_c)); - const f_t sqrt2 = std::sqrt(2.0); - const f_t x_soc = mu * sqrt2; - const f_t z_soc = mu * sqrt2; - // Linear orthant - for (i_t j = 0; j < cs; ++j) { + // Full primal+dual point; no factorization/solve (main loop factorizes later). + if (init_strategy == 2) { + const f_t norm_b = vector_norm_inf(lp.rhs); + const f_t norm_c = vector_norm_inf(lp.objective); + const f_t mu = std::sqrt((1.0 + norm_b) * (1.0 + norm_c)); + const f_t sqrt2 = std::sqrt(2.0); + const i_t linear_end = data.linear_xz_size(lp.num_cols); + + // Linear orthant: x = z = mu * e, with e = 1 + for (i_t j = 0; j < linear_end; ++j) { data.x[j] = mu; data.z[j] = mu; } if (has_direct_free_linear) { for (i_t j : presolve_info.direct_free_variables) { - if (j < cs) { data.z[j] = 0.0; } + if (j < linear_end) { data.z[j] = 0.0; } } } - // SOC blocks - i_t off = 0; - for (size_t k = 0; k < lp.second_order_cone_dims.size(); k++) { - i_t q_k = lp.second_order_cone_dims[k]; - data.x[cs + off] = x_soc; - data.z[cs + off] = z_soc; - for (i_t j = 1; j < q_k; ++j) { - data.x[cs + off + j] = 0.0; - data.z[cs + off + j] = 0.0; + + // SOC blocks: x = z = mu * e, with e = (sqrt(2), 0, ..., 0) + if (data.has_cones()) { + const i_t cs = data.cone_start(); + const f_t x_soc = mu * sqrt2; + const f_t z_soc = mu * sqrt2; + i_t off = 0; + for (size_t k = 0; k < lp.second_order_cone_dims.size(); k++) { + i_t q_k = lp.second_order_cone_dims[k]; + data.x[cs + off] = x_soc; + data.z[cs + off] = z_soc; + for (i_t j = 1; j < q_k; ++j) { + data.x[cs + off + j] = 0.0; + data.z[cs + off + j] = 0.0; + } + off += q_k; } - off += q_k; } + data.y.set_scalar(0.0); if (data.n_upper_bounds > 0) { data.w.set_scalar(mu); diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index d5834f323f..fe8c5474e3 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -164,7 +164,8 @@ struct simplex_solver_settings_t { i_t dualize; // -1 automatic, 0 to not dualize, 1 to dualize i_t ordering; // -1 automatic, 0 to use nested dissection, 1 to use AMD i_t barrier_dual_initial_point; // -1 automatic, 0 to use Lustig, Marsten, and Shanno initial - // point, 1 to use initial point form dual least squares problem + // point, 1 to use initial point form dual least squares problem, + // 2 to use Sturm/SeDuMi mu-based primal+dual point bool check_Q; // true to check if Q is positive semidefinite bool crossover; // true to do crossover, false to not i_t refactor_frequency; // number of basis updates before refactorization diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 8bfd7cd8ba..816414c893 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -131,7 +131,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_FOLDING, &pdlp_settings.folding, -1, 1, -1}, {CUOPT_DUALIZE, &pdlp_settings.dualize, -1, 1, -1}, {CUOPT_ORDERING, &pdlp_settings.ordering, -1, 1, -1}, - {CUOPT_BARRIER_DUAL_INITIAL_POINT, &pdlp_settings.barrier_dual_initial_point, -1, 1, -1}, + {CUOPT_BARRIER_DUAL_INITIAL_POINT, &pdlp_settings.barrier_dual_initial_point, -1, 2, -1}, {CUOPT_MIP_CUT_PASSES, &mip_settings.max_cut_passes, -1, std::numeric_limits::max(), 10}, {CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS, &mip_settings.mir_cuts, -1, 1, -1}, {CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS, &mip_settings.mixed_integer_gomory_cuts, -1, 1, -1}, diff --git a/python/cuopt_server/cuopt_server/tests/test_lp.py b/python/cuopt_server/cuopt_server/tests/test_lp.py index e3a683f8de..8ea85b60ab 100644 --- a/python/cuopt_server/cuopt_server/tests/test_lp.py +++ b/python/cuopt_server/cuopt_server/tests/test_lp.py @@ -180,7 +180,7 @@ def test_barrier_solver_options( - cudss_deterministic: True for deterministic, False for nondeterministic - barrier_dual_initial_point: (-1) automatic, (0) Lustig-Marsten-Shanno, - (1) dual least squares + (1) dual least squares, (2) Sturm/SeDuMi mu-based primal+dual """ data = get_std_data_for_lp() diff --git a/python/cuopt_server/cuopt_server/utils/linear_programming/data_definition.py b/python/cuopt_server/cuopt_server/utils/linear_programming/data_definition.py index 6cd8f7828a..75998a347a 100644 --- a/python/cuopt_server/cuopt_server/utils/linear_programming/data_definition.py +++ b/python/cuopt_server/cuopt_server/utils/linear_programming/data_definition.py @@ -501,7 +501,8 @@ class SolverConfig(BaseModel): description="Set the type of dual initial point to use for the barrier" "solver. -1 for automatic, 0 to use Lustig, Marsten, and Shanno" "initial point, 1 to use initial point from a dual least squares" - "problem", + "problem, 2 to use Sturm/SeDuMi mu-based primal+dual" + "point", ) eliminate_dense_columns: Optional[bool] = Field( default=True, From 8e29f927a73e64046bfc38693e6b4a3f78365eaf Mon Sep 17 00:00:00 2001 From: yuwenchen95 Date: Wed, 15 Jul 2026 09:58:14 -0700 Subject: [PATCH 2/4] Extend dual initial point strategies for socp Signed-off-by: yuwenchen95 --- cpp/src/barrier/barrier.cu | 144 ++++++++++++++++++++++++------------- 1 file changed, 96 insertions(+), 48 deletions(-) diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index e55abddce3..bae7e258ca 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -94,6 +94,39 @@ bool validate_barrier_cone_layout(const lp_problem_t& problem, return true; } +// Push entries into interior of nonnegative orthant and SOC. +template +static void ensure_initial_point_interior(dense_vector_t& values, + f_t epsilon_adjust, + const std::vector& linear_mask, + i_t linear_end, + const std::vector& cone_dims) +{ + f_t min_linear = inf; + for (i_t j = 0; j < linear_end; ++j) { + if (linear_mask[j]) { min_linear = std::min(min_linear, values[j]); } + } + if (min_linear <= epsilon_adjust) { + const f_t delta = -min_linear + epsilon_adjust; + for (i_t j = 0; j < linear_end; ++j) { + if (linear_mask[j]) { values[j] += delta; } + } + } + + i_t off = 0; + for (i_t q_k : cone_dims) { + const i_t base = linear_end + off; + f_t tail_sq = 0.0; + for (i_t j = 1; j < q_k; ++j) { + const f_t t = values[base + j]; + tail_sq += t * t; + } + const f_t tail_norm = std::sqrt(tail_sq); + if (values[base] <= tail_norm + epsilon_adjust) { values[base] = tail_norm + epsilon_adjust; } + off += q_k; + } +} + template [[maybe_unused]] static void pairwise_multiply( f_t* a, f_t* b, f_t* out, int size, rmm::cuda_stream_view stream) @@ -239,7 +272,6 @@ class iteration_data_t { complementarity_residual_norm_save(inf), diag(lp.num_cols), inv_diag(lp.num_cols), - inv_sqrt_diag(lp.num_cols), AD(lp.num_cols, lp.num_rows, 0), AT(lp.num_rows, lp.num_cols, 0), ADAT(lp.num_rows, lp.num_rows, 0), @@ -540,12 +572,9 @@ class iteration_data_t { } inv_diag.set_scalar(1.0); - if (use_augmented) { diag.multiply_scalar(-1.0); } if (n_upper_bounds > 0 || (has_Q && !use_augmented)) { diag.inverse(inv_diag); } // TMP diag and inv_diag should directly created and filled on the GPU raft::copy(d_inv_diag.data(), inv_diag.data(), inv_diag.size(), stream_view_); - inv_sqrt_diag.set_scalar(1.0); - if (n_upper_bounds > 0 || (has_Q && !use_augmented)) { inv_diag.sqrt(inv_sqrt_diag); } } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; } @@ -1912,7 +1941,6 @@ class iteration_data_t { dense_vector_t diag; pinned_dense_vector_t inv_diag; - dense_vector_t inv_sqrt_diag; rmm::device_uvector d_original_A_values; @@ -2237,15 +2265,18 @@ int barrier_solver_t::initial_point(iteration_data_t& data) dense_vector_t DinvFu(lp.num_cols); // DinvFu = Dinv * Fu data.inv_diag.pairwise_product(Fu, DinvFu); dense_vector_t q(lp.num_rows); + const i_t aug_base = lp.num_cols + lp.num_rows; + const i_t aug_size = use_augmented ? std::max(aug_base, data.device_augmented.n) : aug_base; if (use_augmented) { - dense_vector_t rhs(lp.num_cols + lp.num_rows); + dense_vector_t rhs(aug_size); + rhs.set_scalar(0.0); for (i_t k = 0; k < lp.num_cols; k++) { rhs[k] = -Fu[k]; } for (i_t k = 0; k < lp.num_rows; k++) { rhs[lp.num_cols + k] = rhs_x[k]; } - dense_vector_t soln(lp.num_cols + lp.num_rows); + dense_vector_t soln(aug_size); i_t solve_status = data.chol->solve(rhs, soln); struct op_t { op_t(iteration_data_t& data) : data_(data) {} @@ -2263,7 +2294,11 @@ int barrier_solver_t::initial_point(iteration_data_t& data) } } op(data); - if (settings.barrier_iterative_refinement) { + // Initial-point IR is LP/QP-only. Cone problems and direct-free linear vars are excluded. + if (settings.barrier_iterative_refinement && !data.has_cones() && + data.n_direct_free_linear == 0) { + raft::copy( + data.d_diag_.data(), data.diag.data(), data.diag.size(), data.handle_ptr->get_stream()); iterative_refinement(op, rhs, soln); } @@ -2317,27 +2352,19 @@ int barrier_solver_t::initial_point(iteration_data_t& data) } } - // Verify A*x = b - dense_vector_t init_primal_residual(lp.num_rows); - init_primal_residual = lp.rhs; - data.cusparse_view_.spmv(1.0, data.x, -1.0, init_primal_residual); - data.handle_ptr->get_stream().synchronize(); -#ifdef PRINT_INFO - settings.log.printf("||b - A * x||: %.16e\n", vector_norm2(init_primal_residual)); -#endif - - if (data.n_upper_bounds > 0) { - dense_vector_t init_bound_residual(data.n_upper_bounds); - for (i_t k = 0; k < data.n_upper_bounds; k++) { - i_t j = data.upper_bounds[k]; - init_bound_residual[k] = lp.upper[j] - data.w[k] - data.x[j]; - } -#ifdef PRINT_INFO - settings.log.printf("|| u - w - x||: %e\n", vector_norm2(init_bound_residual)); -#endif - } - float64_t epsilon_adjust = 10.0; + // Push entries into interior of nonnegative orthant and SOC. + const bool has_soc = data.has_cones(); + const i_t linear_end = has_soc ? data.cone_start() : lp.num_cols; + auto ensure_interior = [&](dense_vector_t& values, + const std::vector& linear_mask) { + if (has_soc) { + ensure_initial_point_interior( + values, epsilon_adjust, linear_mask, linear_end, lp.second_order_cone_dims); + } else { + values.ensure_positive(epsilon_adjust, linear_mask); + } + }; if (settings.barrier_dual_initial_point == -1 || settings.barrier_dual_initial_point == 0) { // Use the dual starting point suggested by the paper @@ -2388,12 +2415,12 @@ int barrier_solver_t::initial_point(iteration_data_t& data) } } } else if (use_augmented) { - dense_vector_t dual_rhs(lp.num_cols + lp.num_rows); + dense_vector_t dual_rhs(aug_size); dual_rhs.set_scalar(0.0); for (i_t k = 0; k < lp.num_cols; k++) { dual_rhs[k] = data.c[k]; } - dense_vector_t py(lp.num_cols + lp.num_rows); + dense_vector_t py(aug_size); data.chol->solve(dual_rhs, py); for (i_t k = 0; k < lp.num_cols; k++) { data.z[k] = py[k]; @@ -2407,7 +2434,6 @@ int barrier_solver_t::initial_point(iteration_data_t& data) data.v.multiply_scalar(-1.0); data.v.ensure_positive(epsilon_adjust); - data.z.ensure_positive(epsilon_adjust, nonnegative_z); } else { // First compute rhs = A*Dinv*c dense_vector_t rhs(lp.num_rows); @@ -2431,24 +2457,8 @@ int barrier_solver_t::initial_point(iteration_data_t& data) data.gather_upper_bounds(data.z, data.v); data.v.multiply_scalar(-1.0); data.v.ensure_positive(epsilon_adjust); - data.z.ensure_positive(epsilon_adjust, nonnegative_z); } - // Verify A'*y + z - E*v - Q*x = c - dense_vector_t init_dual_residual(lp.num_cols); - data.z.pairwise_subtract(data.c, init_dual_residual); - if (data.Q.n > 0) { matrix_vector_multiply(data.Q, -1.0, data.x, 1.0, init_dual_residual); } - data.cusparse_view_.transpose_spmv(1.0, data.y, 1.0, init_dual_residual); - if (data.n_upper_bounds > 0) { - for (i_t k = 0; k < data.n_upper_bounds; k++) { - i_t j = data.upper_bounds[k]; - init_dual_residual[j] -= data.v[k]; - } - } -#ifdef PRINT_INFO - settings.log.printf("||A^T y + z - E*v - Q*x - c ||: %e\n", - vector_norm2(init_dual_residual)); -#endif // Make sure (w, x, v, z) > 0. Skip free variables being handled directly. data.w.ensure_positive(epsilon_adjust); std::vector nonnegative_variables(data.x.size(), 1); @@ -2457,7 +2467,8 @@ int barrier_solver_t::initial_point(iteration_data_t& data) nonnegative_variables[j] = 0; } } - data.x.ensure_positive(epsilon_adjust, nonnegative_variables); + ensure_interior(data.z, nonnegative_z); + ensure_interior(data.x, nonnegative_variables); // Direct free variables: reduced cost z = 0 (no complementarity condition). if (has_direct_free_linear) { for (i_t j : presolve_info.direct_free_variables) { @@ -2468,6 +2479,43 @@ int barrier_solver_t::initial_point(iteration_data_t& data) settings.log.printf("min v %e min z %e\n", data.v.minimum(), data.z.minimum()); #endif + // Residual checks below reflect the final initial point, after positivity shifts. + // Verify A*x = b + dense_vector_t init_primal_residual(lp.num_rows); + init_primal_residual = lp.rhs; + data.cusparse_view_.spmv(1.0, data.x, -1.0, init_primal_residual); + data.handle_ptr->get_stream().synchronize(); +#ifdef PRINT_INFO + settings.log.printf("||b - A * x||: %.16e\n", vector_norm2(init_primal_residual)); +#endif + + if (data.n_upper_bounds > 0) { + dense_vector_t init_bound_residual(data.n_upper_bounds); + for (i_t k = 0; k < data.n_upper_bounds; k++) { + i_t j = data.upper_bounds[k]; + init_bound_residual[k] = lp.upper[j] - data.w[k] - data.x[j]; + } +#ifdef PRINT_INFO + settings.log.printf("|| u - w - x||: %e\n", vector_norm2(init_bound_residual)); +#endif + } + + // Verify A'*y + z - E*v - Q*x = c + dense_vector_t init_dual_residual(lp.num_cols); + data.z.pairwise_subtract(data.c, init_dual_residual); + if (data.Q.n > 0) { matrix_vector_multiply(data.Q, -1.0, data.x, 1.0, init_dual_residual); } + data.cusparse_view_.transpose_spmv(1.0, data.y, 1.0, init_dual_residual); + if (data.n_upper_bounds > 0) { + for (i_t k = 0; k < data.n_upper_bounds; k++) { + i_t j = data.upper_bounds[k]; + init_dual_residual[j] -= data.v[k]; + } + } +#ifdef PRINT_INFO + settings.log.printf("||A^T y + z - E*v - Q*x - c ||: %e\n", + vector_norm2(init_dual_residual)); +#endif + return 0; } From d39ddefa3986a4e3361665fbb048e3d4151c218d Mon Sep 17 00:00:00 2001 From: yuwenchen95 Date: Thu, 16 Jul 2026 11:37:56 -0700 Subject: [PATCH 3/4] Use ENUM for initial point strategies Signed-off-by: yuwenchen95 --- .../cuopt/mathematical_optimization/constants.h | 5 +++++ .../utilities/internals.hpp | 15 +++++++++++++++ cpp/src/barrier/barrier.cu | 16 ++++++++++++---- cpp/src/dual_simplex/simplex_solver_settings.hpp | 4 +--- 4 files changed, 33 insertions(+), 7 deletions(-) diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 5e720a645e..8d896d38e9 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -183,6 +183,11 @@ #define CUOPT_METHOD_BARRIER 3 #define CUOPT_METHOD_UNSET 4 +#define CUOPT_BARRIER_DUAL_INITIAL_POINT_AUTOMATIC -1 +#define CUOPT_BARRIER_DUAL_INITIAL_POINT_LUSTIG_MARSTEN_SHANNO 0 +#define CUOPT_BARRIER_DUAL_INITIAL_POINT_DUAL_LEAST_SQUARES 1 +#define CUOPT_BARRIER_DUAL_INITIAL_POINT_SEDUMI_MU 2 + /* @brief PDLP precision mode constants */ #define CUOPT_PDLP_DEFAULT_PRECISION -1 #define CUOPT_PDLP_SINGLE_PRECISION 0 diff --git a/cpp/include/cuopt/mathematical_optimization/utilities/internals.hpp b/cpp/include/cuopt/mathematical_optimization/utilities/internals.hpp index aaec8ef842..f5602b0ff0 100644 --- a/cpp/include/cuopt/mathematical_optimization/utilities/internals.hpp +++ b/cpp/include/cuopt/mathematical_optimization/utilities/internals.hpp @@ -142,5 +142,20 @@ enum presolver_t : int { PSLP = CUOPT_PRESOLVE_PSLP }; +/** + * @brief Barrier primal-dual initial-point strategy. + * + * Automatic: use Lustig-Marsten-Shanno for LP/QP; Sturm/SeDuMi mu-based point for conic problems. + * LustigMarstenShanno: Mehrotra-style dual start (Lustig, Marsten, Shanno, SIAM J. Optim. 1992). + * DualLeastSquares: solve augmented or ADAT dual least-squares system. + * SedumiMu: Sturm/SeDuMi mu-based primal+dual point (no factorization). + */ +enum barrier_dual_initial_point_t : int { + Automatic = CUOPT_BARRIER_DUAL_INITIAL_POINT_AUTOMATIC, + LustigMarstenShanno = CUOPT_BARRIER_DUAL_INITIAL_POINT_LUSTIG_MARSTEN_SHANNO, + DualLeastSquares = CUOPT_BARRIER_DUAL_INITIAL_POINT_DUAL_LEAST_SQUARES, + SedumiMu = CUOPT_BARRIER_DUAL_INITIAL_POINT_SEDUMI_MU +}; + } // namespace mathematical_optimization } // namespace cuopt diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index bae7e258ca..563292206c 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -22,6 +22,7 @@ #include #include +#include #include #include #include @@ -2183,14 +2184,20 @@ int barrier_solver_t::initial_point(iteration_data_t& data) const bool use_augmented = data.use_augmented; const bool has_direct_free_linear = data.n_direct_free_linear > 0; - const i_t init_strategy = data.has_cones() ? 2 : settings.barrier_dual_initial_point; + const barrier_dual_initial_point_t input_strategy = + static_cast(settings.barrier_dual_initial_point); - // Option 2: Sturm/SeDuMi-style mu-based primal+dual initial point. + const barrier_dual_initial_point_t init_strategy = + (data.has_cones() && input_strategy == barrier_dual_initial_point_t::Automatic) + ? barrier_dual_initial_point_t::SedumiMu + : input_strategy; + + // SedumiMu: Sturm/SeDuMi-style mu-based primal+dual initial point. // mu = sqrt((1 + ||b||_inf) * (1 + ||c||_inf)); x = z = mu * e_K. // where e_K is the identity of the symmetric cone: // LP block: e = 1, SOC block: e = (sqrt(2), 0, ..., 0) // Full primal+dual point; no factorization/solve (main loop factorizes later). - if (init_strategy == 2) { + if (init_strategy == barrier_dual_initial_point_t::SedumiMu) { const f_t norm_b = vector_norm_inf(lp.rhs); const f_t norm_c = vector_norm_inf(lp.objective); const f_t mu = std::sqrt((1.0 + norm_b) * (1.0 + norm_c)); @@ -2366,7 +2373,8 @@ int barrier_solver_t::initial_point(iteration_data_t& data) } }; - if (settings.barrier_dual_initial_point == -1 || settings.barrier_dual_initial_point == 0) { + if (init_strategy == barrier_dual_initial_point_t::Automatic || + init_strategy == barrier_dual_initial_point_t::LustigMarstenShanno) { // Use the dual starting point suggested by the paper // On Implementing Mehrotra’s Predictor–Corrector Interior-Point Method for Linear Programming // Irvin J. Lustig, Roy E. Marsten, and David F. Shanno diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index fe8c5474e3..821ba462ef 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -163,9 +163,7 @@ struct simplex_solver_settings_t { i_t augmented; // -1 automatic, 0 to solve with ADAT, 1 to solve with augmented system i_t dualize; // -1 automatic, 0 to not dualize, 1 to dualize i_t ordering; // -1 automatic, 0 to use nested dissection, 1 to use AMD - i_t barrier_dual_initial_point; // -1 automatic, 0 to use Lustig, Marsten, and Shanno initial - // point, 1 to use initial point form dual least squares problem, - // 2 to use Sturm/SeDuMi mu-based primal+dual point + i_t barrier_dual_initial_point; // barrier_dual_initial_point_t; see internals.hpp bool check_Q; // true to check if Q is positive semidefinite bool crossover; // true to do crossover, false to not i_t refactor_frequency; // number of basis updates before refactorization From 4f5e8fde594797eb40c5719624af412b9fb28665 Mon Sep 17 00:00:00 2001 From: YUWEN Chen Date: Thu, 23 Jul 2026 07:42:59 -0700 Subject: [PATCH 4/4] restore residual check in initialization and clean up code related to initial point Signed-off-by: YUWEN Chen --- .../utilities/internals.hpp | 15 ---- cpp/src/barrier/barrier.cu | 74 +++++++++---------- cpp/src/barrier/barrier.hpp | 16 ++++ .../dual_simplex/simplex_solver_settings.hpp | 3 +- 4 files changed, 54 insertions(+), 54 deletions(-) diff --git a/cpp/include/cuopt/mathematical_optimization/utilities/internals.hpp b/cpp/include/cuopt/mathematical_optimization/utilities/internals.hpp index f5602b0ff0..aaec8ef842 100644 --- a/cpp/include/cuopt/mathematical_optimization/utilities/internals.hpp +++ b/cpp/include/cuopt/mathematical_optimization/utilities/internals.hpp @@ -142,20 +142,5 @@ enum presolver_t : int { PSLP = CUOPT_PRESOLVE_PSLP }; -/** - * @brief Barrier primal-dual initial-point strategy. - * - * Automatic: use Lustig-Marsten-Shanno for LP/QP; Sturm/SeDuMi mu-based point for conic problems. - * LustigMarstenShanno: Mehrotra-style dual start (Lustig, Marsten, Shanno, SIAM J. Optim. 1992). - * DualLeastSquares: solve augmented or ADAT dual least-squares system. - * SedumiMu: Sturm/SeDuMi mu-based primal+dual point (no factorization). - */ -enum barrier_dual_initial_point_t : int { - Automatic = CUOPT_BARRIER_DUAL_INITIAL_POINT_AUTOMATIC, - LustigMarstenShanno = CUOPT_BARRIER_DUAL_INITIAL_POINT_LUSTIG_MARSTEN_SHANNO, - DualLeastSquares = CUOPT_BARRIER_DUAL_INITIAL_POINT_DUAL_LEAST_SQUARES, - SedumiMu = CUOPT_BARRIER_DUAL_INITIAL_POINT_SEDUMI_MU -}; - } // namespace mathematical_optimization } // namespace cuopt diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index 563292206c..e7af2a294b 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -22,7 +22,6 @@ #include #include -#include #include #include #include @@ -2359,6 +2358,26 @@ int barrier_solver_t::initial_point(iteration_data_t& data) } } + // Verify A*x = b + dense_vector_t init_primal_residual(lp.num_rows); + init_primal_residual = lp.rhs; + data.cusparse_view_.spmv(1.0, data.x, -1.0, init_primal_residual); + data.handle_ptr->get_stream().synchronize(); +#ifdef PRINT_INFO + settings.log.printf("||b - A * x||: %.16e\n", vector_norm2(init_primal_residual)); +#endif + + if (data.n_upper_bounds > 0) { + dense_vector_t init_bound_residual(data.n_upper_bounds); + for (i_t k = 0; k < data.n_upper_bounds; k++) { + i_t j = data.upper_bounds[k]; + init_bound_residual[k] = lp.upper[j] - data.w[k] - data.x[j]; + } +#ifdef PRINT_INFO + settings.log.printf("|| u - w - x||: %e\n", vector_norm2(init_bound_residual)); +#endif + } + float64_t epsilon_adjust = 10.0; // Push entries into interior of nonnegative orthant and SOC. const bool has_soc = data.has_cones(); @@ -2467,6 +2486,22 @@ int barrier_solver_t::initial_point(iteration_data_t& data) data.v.ensure_positive(epsilon_adjust); } + // Verify A'*y + z - E*v - Q*x = c + dense_vector_t init_dual_residual(lp.num_cols); + data.z.pairwise_subtract(data.c, init_dual_residual); + if (data.Q.n > 0) { matrix_vector_multiply(data.Q, -1.0, data.x, 1.0, init_dual_residual); } + data.cusparse_view_.transpose_spmv(1.0, data.y, 1.0, init_dual_residual); + if (data.n_upper_bounds > 0) { + for (i_t k = 0; k < data.n_upper_bounds; k++) { + i_t j = data.upper_bounds[k]; + init_dual_residual[j] -= data.v[k]; + } + } +#ifdef PRINT_INFO + settings.log.printf("||A^T y + z - E*v - Q*x - c ||: %e\n", + vector_norm2(init_dual_residual)); +#endif + // Make sure (w, x, v, z) > 0. Skip free variables being handled directly. data.w.ensure_positive(epsilon_adjust); std::vector nonnegative_variables(data.x.size(), 1); @@ -2487,43 +2522,6 @@ int barrier_solver_t::initial_point(iteration_data_t& data) settings.log.printf("min v %e min z %e\n", data.v.minimum(), data.z.minimum()); #endif - // Residual checks below reflect the final initial point, after positivity shifts. - // Verify A*x = b - dense_vector_t init_primal_residual(lp.num_rows); - init_primal_residual = lp.rhs; - data.cusparse_view_.spmv(1.0, data.x, -1.0, init_primal_residual); - data.handle_ptr->get_stream().synchronize(); -#ifdef PRINT_INFO - settings.log.printf("||b - A * x||: %.16e\n", vector_norm2(init_primal_residual)); -#endif - - if (data.n_upper_bounds > 0) { - dense_vector_t init_bound_residual(data.n_upper_bounds); - for (i_t k = 0; k < data.n_upper_bounds; k++) { - i_t j = data.upper_bounds[k]; - init_bound_residual[k] = lp.upper[j] - data.w[k] - data.x[j]; - } -#ifdef PRINT_INFO - settings.log.printf("|| u - w - x||: %e\n", vector_norm2(init_bound_residual)); -#endif - } - - // Verify A'*y + z - E*v - Q*x = c - dense_vector_t init_dual_residual(lp.num_cols); - data.z.pairwise_subtract(data.c, init_dual_residual); - if (data.Q.n > 0) { matrix_vector_multiply(data.Q, -1.0, data.x, 1.0, init_dual_residual); } - data.cusparse_view_.transpose_spmv(1.0, data.y, 1.0, init_dual_residual); - if (data.n_upper_bounds > 0) { - for (i_t k = 0; k < data.n_upper_bounds; k++) { - i_t j = data.upper_bounds[k]; - init_dual_residual[j] -= data.v[k]; - } - } -#ifdef PRINT_INFO - settings.log.printf("||A^T y + z - E*v - Q*x - c ||: %e\n", - vector_norm2(init_dual_residual)); -#endif - return 0; } diff --git a/cpp/src/barrier/barrier.hpp b/cpp/src/barrier/barrier.hpp index 07094a54a5..57e691be71 100644 --- a/cpp/src/barrier/barrier.hpp +++ b/cpp/src/barrier/barrier.hpp @@ -8,6 +8,7 @@ #include +#include #include #include #include @@ -18,6 +19,21 @@ #include namespace cuopt::mathematical_optimization::barrier { +/** + * @brief Barrier primal-dual initial-point strategy. + * + * Automatic: use Lustig-Marsten-Shanno for LP/QP; Sturm/SeDuMi mu-based point for conic problems. + * LustigMarstenShanno: Mehrotra-style dual start (Lustig, Marsten, Shanno, SIAM J. Optim. 1992). + * DualLeastSquares: solve augmented or ADAT dual least-squares system. + * SedumiMu: Sturm/SeDuMi mu-based primal+dual point (no factorization). + */ +enum barrier_dual_initial_point_t : int { + Automatic = CUOPT_BARRIER_DUAL_INITIAL_POINT_AUTOMATIC, + LustigMarstenShanno = CUOPT_BARRIER_DUAL_INITIAL_POINT_LUSTIG_MARSTEN_SHANNO, + DualLeastSquares = CUOPT_BARRIER_DUAL_INITIAL_POINT_DUAL_LEAST_SQUARES, + SedumiMu = CUOPT_BARRIER_DUAL_INITIAL_POINT_SEDUMI_MU +}; + /** Validates SOC layout on an simplex::lp_problem_t before barrier presolve/solve. */ template bool validate_barrier_cone_layout(const simplex::lp_problem_t& problem, diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index ca2fc53e0a..646397513e 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -167,7 +167,8 @@ struct simplex_solver_settings_t { i_t augmented; // -1 automatic, 0 to solve with ADAT, 1 to solve with augmented system i_t dualize; // -1 automatic, 0 to not dualize, 1 to dualize i_t ordering; // -1 automatic, 0 to use nested dissection, 1 to use AMD - i_t barrier_dual_initial_point; // barrier_dual_initial_point_t; see internals.hpp + i_t barrier_dual_initial_point; // -1 automatic, 0 Lustig-Marsten-Shanno, + // 1 dual least squares, 2 SeDuMi mu-based bool check_Q; // true to check if Q is positive semidefinite bool crossover; // true to do crossover, false to not i_t refactor_frequency; // number of basis updates before refactorization