diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index f6be07aaa9..9e57e1feed 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -194,6 +194,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/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index 784f6c0901..e7af2a294b 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; @@ -2155,41 +2183,55 @@ 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 barrier_dual_initial_point_t input_strategy = + static_cast(settings.barrier_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) - 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 == 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)); + 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); @@ -2229,15 +2271,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) {} @@ -2255,7 +2300,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); } @@ -2330,8 +2379,21 @@ int barrier_solver_t::initial_point(iteration_data_t& data) } 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) { + 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 @@ -2380,12 +2442,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]; @@ -2399,7 +2461,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); @@ -2423,7 +2484,6 @@ 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 @@ -2441,6 +2501,7 @@ int barrier_solver_t::initial_point(iteration_data_t& data) 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); @@ -2449,7 +2510,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) { 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 6a69cdfcd2..646397513e 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -167,8 +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; // -1 automatic, 0 to use Lustig, Marsten, and Shanno initial - // point, 1 to use initial point form dual least squares problem + 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 diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 9193112d71..1010b1f8c6 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -137,7 +137,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,