diff --git a/cpp/cmake/thirdparty/FindCUDSS.cmake b/cpp/cmake/thirdparty/FindCUDSS.cmake index 29021142c6..cf9fbdfcde 100644 --- a/cpp/cmake/thirdparty/FindCUDSS.cmake +++ b/cpp/cmake/thirdparty/FindCUDSS.cmake @@ -1,6 +1,6 @@ # ============================================================================= # cmake-format: off -# SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION. +# SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION. # SPDX-License-Identifier: Apache-2.0 # cmake-format: on # ============================================================================= diff --git a/cpp/src/barrier/translate_soc.hpp b/cpp/src/barrier/translate_soc.hpp index 3209dee3c6..bd2dc26a88 100644 --- a/cpp/src/barrier/translate_soc.hpp +++ b/cpp/src/barrier/translate_soc.hpp @@ -16,6 +16,8 @@ #include #include +#include + #include #include #include @@ -63,6 +65,10 @@ void convert_quadratic_constraints_to_second_order_cones( error_type_t::ValidationError, "Quadratic-constraint flag is set, but no constraints were provided"); + user_problem.original_num_rows = csr_A.m; + user_problem.qc_dual_recovery.clear(); + user_problem.qc_dual_recovery.resize(qcs.size()); + // Use a practical tolerance for text-parsed MPS numeric values. const f_t tol = std::numeric_limits::epsilon() * 2; @@ -564,14 +570,25 @@ void convert_quadratic_constraints_to_second_order_cones( cone_vars.push_back(std::move(cone)); cone_is_rotated.push_back(is_rotated); + auto& dual_meta = user_problem.qc_dual_recovery[qc_i]; + dual_meta.uniform_s = uniform_s; + dual_meta.cone_index = static_cast(cone_dims.size()) - 1; + if (has_linear_part) { + dual_meta.path = dual_simplex::user_problem_t::qc_soc_recognition_path_t::AFFINE; + } else if (is_rotated) { + dual_meta.path = dual_simplex::user_problem_t::qc_soc_recognition_path_t::ROTATED; + } else { + dual_meta.path = dual_simplex::user_problem_t::qc_soc_recognition_path_t::LORENTZ; + } + } else { // ========================================================================= - // General convex quadratic constraint path: - // x^T Q x + c^T x <= alpha + // General convex quadratic constraint path (QCMATRIX storage, no extra 1/2 on quad part): + // sum_k v_k x_{r_k} x_{c_k} + c^T x <= alpha // Invalidate affine head for this QC — the general path handles the linear part directly qc_affine_heads[qc_i] = -1; - // where Q is (possibly unsymmetric) and H = Q + Q^T must be PSD. + // H below satisfies (1/2) x^T H x = that QCMATRIX sum; H must be PSD. // ========================================================================= const f_t alpha = qc.rhs_value; @@ -798,6 +815,13 @@ void convert_quadratic_constraints_to_second_order_cones( cone_dims.push_back(cone_dim); cone_vars.push_back(std::move(cone)); cone_is_rotated.push_back(is_rotated); + + auto& dual_meta = user_problem.qc_dual_recovery[qc_i]; + dual_meta.path = dual_simplex::user_problem_t::qc_soc_recognition_path_t::GENERAL; + dual_meta.uniform_s = 1; + dual_meta.cone_index = static_cast(cone_dims.size()) - 1; + dual_meta.s0_link_row = m_before + r; + dual_meta.sr1_link_row = m_before + r + 1; } } @@ -865,6 +889,7 @@ void convert_quadratic_constraints_to_second_order_cones( user_problem.row_names[row_write_cursor] = "_CUOPT_qc_linear_link_" + qc.constraint_row_name; } + user_problem.qc_dual_recovery[qc_i].affine_link_row = row_write_cursor; ++row_write_cursor; } @@ -965,10 +990,13 @@ void convert_quadratic_constraints_to_second_order_cones( is_cone_var[slack_base] = 1; is_cone_var[slack_base + 1] = 1; - eq_row.n = n_prob; + eq_row.n = n_prob; + auto& rsoc_meta = user_problem.qc_dual_recovery[static_cast(ci)]; // If the second head is not constant half, we need to lift it. if (!rc.head1_is_constant_half) { - const f_t h = inv_sqrt_2 * rc.head_lift_sqrt_ratio; + const f_t h = inv_sqrt_2 * rc.head_lift_sqrt_ratio; + rsoc_meta.rsoc_head_lift_h = h; + rsoc_meta.rsoc_head1_is_constant_half = false; // s_0 - h * x_h0 - h * x_h1 = 0 (h = inv_sqrt_2 * sqrt(d/s)) eq_row.i = {rc.head0, rc.head1, slack_base}; eq_row.x = {-h, -h, f_t(1)}; @@ -976,6 +1004,7 @@ void convert_quadratic_constraints_to_second_order_cones( csr_A.append_row(eq_row); user_problem.row_sense[row_idx] = 'E'; user_problem.rhs[row_idx] = 0; + rsoc_meta.rsoc_s0_lift_row = row_idx; ++row_idx; // s_1 - h * x_h0 + h * x_h1 = 0 @@ -985,12 +1014,15 @@ void convert_quadratic_constraints_to_second_order_cones( csr_A.append_row(eq_row); user_problem.row_sense[row_idx] = 'E'; user_problem.rhs[row_idx] = 0; + rsoc_meta.rsoc_s1_lift_row = row_idx; ++row_idx; is_cone_var[rc.head0] = 0; is_cone_var[rc.head1] = 0; } else { // One head is constant half, so we can lift it directly. + rsoc_meta.rsoc_head_lift_h = inv_sqrt_2; + rsoc_meta.rsoc_head1_is_constant_half = true; // s_0 - inv_sqrt_2 * x_h0 = inv_sqrt_2 * (1/2) eq_row.i = {rc.head0, slack_base}; eq_row.x = {-inv_sqrt_2, f_t(1)}; @@ -998,6 +1030,7 @@ void convert_quadratic_constraints_to_second_order_cones( csr_A.append_row(eq_row); user_problem.row_sense[row_idx] = 'E'; user_problem.rhs[row_idx] = inv_sqrt_2 * half; + rsoc_meta.rsoc_s0_lift_row = row_idx; ++row_idx; // s_1 - inv_sqrt_2 * x_h0 = -inv_sqrt_2 * (1/2) @@ -1007,6 +1040,7 @@ void convert_quadratic_constraints_to_second_order_cones( csr_A.append_row(eq_row); user_problem.row_sense[row_idx] = 'E'; user_problem.rhs[row_idx] = -inv_sqrt_2 * half; + rsoc_meta.rsoc_s1_lift_row = row_idx; ++row_idx; is_cone_var[rc.head0] = 0; @@ -1352,4 +1386,152 @@ void project_barrier_solution_to_model_variables( solution.z = std::move(model_z); } +/** + * True when the Lorentz cone block is at (or numerically indistinguishable from) the apex. + * Feasibility with tiny head implies tiny tails; grad g ≈ 0 so the QC multiplier is not unique. + */ +template +bool lorentz_cone_block_at_apex(const std::vector& x, + i_t cone_offset, + i_t cone_dim, + f_t apex_tol = 1e-6) +{ + f_t cone_inf = 0; + for (i_t k = 0; k < cone_dim; ++k) { + cone_inf = std::max(cone_inf, std::abs(x[cone_offset + k])); + } + return cone_inf <= apex_tol; +} + +/** + * Recover the SOC KKT multiplier from barrier primal/reduced-cost pair (x_head, z_head). + * + * Barrier Lorentz block (uniform scale s > 0): + * phi = -s * x_head^2 + s * sum_{i in tail} x_i^2 <= 0. + * + * At a non-apex point, cone stationarity in the barrier's (x, z) coordinates gives + * nu_bar = z_head / (2 * s * x_head), + * where nu_bar is the cone dual in barrier export coordinates (not necessarily the textbook + * multiplier nu on phi <= 0 with the same sign). LORENTZ / ROTATED / AFFINE paths validate + * that nu_bar equals the exported user QC dual when the user row is already in phi-form. + * + * At the apex (||x||_inf <= apex_tol), grad phi = 0 so the multiplier is not unique; return 0. + * + * --- GENERAL (LDLT) path: extra factor of 2 (see path_t::GENERAL in dual projection) --- + * + * User QC (QCMATRIX / Python API) is stored without an extra 1/2 on the quadratic part: + * g(x) = sum_k v_k x_{r_k} x_{c_k} + c^T x - alpha <= 0 + * (MPS QCMATRIX coeffs are literal; unlike QUADOBJ/QMATRIX objectives where MPS uses 1/2 x^T Q x + * and the parser applies a 0.5 scale — see mps_parser.cpp value_scale notes.) + * LDLT builds symmetric H so (1/2) x^T H x equals that QCMATRIX sum; equivalently + * g(x) = (1/2) x^T H x + c^T x - alpha <= 0. + * LDLT lift (rank r) adds y, s_0, s_{r+1} and equalities + * y_k = sqrt(D_k) (L^T P x)_k, + * s_0 + c^T x = alpha + 1/2, s_{r+1} + c^T x = alpha - 1/2. + * On that manifold: (1/2)||y||^2 = (1/2) x^T H x and s_0 - s_{r+1} = 1, hence + * g = (1/2)||y||^2 + 1/2 - s_0 and phi = -s_0^2 + ||y||^2 + s_{r+1}^2 <= 0 + * describe the same QC inequality (g <= 0 <=> phi <= 0). + * + * Sensitivity / scaling on the LDLT coords (not MPS 1/2 objective vs constraint, and not a + * joint (lambda + 2*nu) = 0 identity with textbook signs): d g / d y_k = y_k because + * g carries (1/2)||y||^2, while d phi / d y_k = 2 y_k on tails. With lambda_user >= 0 + * the exported KKT multiplier on g and nu_bar from the head formula satisfy + * lambda_user = 2 * nu_bar (sign included: both are >= 0 at an active optimum in tests). + * + * nu_bar is recovered from (x, z) on the packed cone head via the formula above; export + * lambda_user = 2 * nu_bar = z_head / (s * x_head). + */ +template +f_t qc_multiplier_from_lorentz_soc_kkt(const std::vector& x, + const std::vector& z, + i_t cone_offset, + i_t cone_dim, + f_t uniform_s, + f_t apex_tol = 1e-6) +{ + if (cone_dim <= 0 || uniform_s <= 0) { return f_t(0); } + const i_t head_col = cone_offset; + cuopt_expects(head_col >= 0 && head_col < static_cast(x.size()), + error_type_t::RuntimeError, + "SOC head column index out of range"); + cuopt_expects(static_cast(z.size()) == static_cast(x.size()), + error_type_t::RuntimeError, + "SOC KKT recovery requires x and z of equal length"); + // primal solution is 0, so mu is not unique and we return 0. + if (lorentz_cone_block_at_apex(x, cone_offset, cone_dim, apex_tol)) { return f_t(0); } + + // otherwise, we can compute mu = z_head / (2 s x_head) + const f_t denom_scale = f_t(2) * uniform_s; + return z[head_col] / (denom_scale * x[head_col]); +} + +/** + * Map barrier duals from the expanded SOC model to the original QCQP layout. + * + * Must run before @ref project_barrier_solution_to_model_variables while @p solution.z still + * includes the trailing cone block (cone duals live in z[cone_var_start:)). + * + * Output layout for @p solution.y: + * [0, original_num_rows) — dual multipliers for user linear rows + * [original_num_rows, original_num_rows + n_qc) — KKT multiplier mu >= 0 on each g(x) <= 0 + */ +template +void project_barrier_qcqp_duals_to_model(const dual_simplex::user_problem_t& user_problem, + dual_simplex::lp_solution_t& solution) +{ + const i_t m_linear = user_problem.original_num_rows; + const i_t n_qc = static_cast(user_problem.qc_dual_recovery.size()); + if (m_linear <= 0 && n_qc == 0) { return; } + // Duals are in expanded-model row order (same layout as uncrushed barrier y). User linear + // rows stay at indices [0, m_linear); QC multipliers are appended after projection. + if (static_cast(solution.y.size()) < m_linear) { return; } + if (static_cast(solution.z.size()) < user_problem.num_cols) { return; } + if (user_problem.second_order_cone_dims.size() != static_cast(n_qc)) { return; } + + const i_t n_out = m_linear + n_qc; + std::vector model_y(n_out, f_t(0)); + for (i_t i = 0; i < m_linear && i < static_cast(solution.y.size()); ++i) { + model_y[i] = solution.y[i]; + } + + i_t cone_offset = user_problem.cone_var_start; + using path_t = typename dual_simplex::user_problem_t::qc_soc_recognition_path_t; + + for (i_t qi = 0; qi < n_qc; ++qi) { + const auto& entry = user_problem.qc_dual_recovery[static_cast(qi)]; + cuopt_expects(entry.cone_index >= 0 && + entry.cone_index < static_cast(user_problem.second_order_cone_dims.size()), + error_type_t::RuntimeError, + "Invalid cone index %d for quadratic constraint dual recovery", + static_cast(entry.cone_index)); + const i_t dim = user_problem.second_order_cone_dims[static_cast(entry.cone_index)]; + cuopt_expects(cone_offset + dim <= user_problem.num_cols, + error_type_t::RuntimeError, + "Cone dual block exceeds expanded column count"); + const f_t s = entry.uniform_s > 0 ? entry.uniform_s : f_t(1); + f_t lambda = 0; + + switch (entry.path) { + case path_t::LORENTZ: + case path_t::ROTATED: + case path_t::AFFINE: + // Lorentz-head recovery on [s0, s1, tails…]. Link-row duals are equality multipliers. + lambda = + qc_multiplier_from_lorentz_soc_kkt(solution.x, solution.z, cone_offset, dim, s); + break; + case path_t::GENERAL: + // User g uses (1/2)||y||^2; expanded phi uses ||y||^2 on tails => lambda_user = 2*nu_bar. + // nu_bar from cone head (z, x); see @ref qc_multiplier_from_lorentz_soc_kkt comment block. + lambda = f_t(2) * qc_multiplier_from_lorentz_soc_kkt( + solution.x, solution.z, cone_offset, dim, s); + break; + } + + model_y[m_linear + qi] = lambda; + cone_offset += dim; + } + + solution.y = std::move(model_y); +} + } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/dual_simplex/user_problem.hpp b/cpp/src/dual_simplex/user_problem.hpp index 548ec9e449..e908dcc32e 100644 --- a/cpp/src/dual_simplex/user_problem.hpp +++ b/cpp/src/dual_simplex/user_problem.hpp @@ -69,6 +69,35 @@ struct user_problem_t { // expanded layout (num_cols) and must be projected back via original_col_to_expanded_col. i_t original_num_cols{0}; std::vector original_col_to_expanded_col; + // Linear constraint count before QCMATRIX->SOC expansion (user-facing rows). + i_t original_num_rows{0}; + /** How each quadratic constraint was recognized during SOC conversion. */ + enum class qc_soc_recognition_path_t : int8_t { + LORENTZ = 0, + AFFINE = 1, + ROTATED = 2, + GENERAL = 3, + }; + /** Per-QC metadata for @ref project_barrier_qcqp_duals_to_model (barrier SOCP -> user QCQP). */ + struct qc_dual_recovery_entry_t { + qc_soc_recognition_path_t path{qc_soc_recognition_path_t::LORENTZ}; + /** Uniform diagonal scale s in -s x_head^2 + s||tail||^2 <= 0 (Lorentz/rotated/affine). */ + f_t uniform_s{1}; + i_t cone_index{-1}; + /** Expanded-model equality row for t = -(1/s)a^T x (affine path); -1 if none. */ + i_t affine_link_row{-1}; + /** Expanded-model rows for s0 + c^T x = alpha + 1/2 and s_{r+1} + c^T x = alpha - 1/2. */ + i_t s0_link_row{-1}; + i_t sr1_link_row{-1}; + /** Rotated lift: equality rows s0 = h(x_h0+x_h1), s1 = h(x_h0-x_h1) (-1 if not rotated). */ + i_t rsoc_s0_lift_row{-1}; + i_t rsoc_s1_lift_row{-1}; + /** h = (1/sqrt(2))*sqrt(d/s) in the two-head lift; inv_sqrt(2) for constant-half affine lift. + */ + f_t rsoc_head_lift_h{0}; + bool rsoc_head1_is_constant_half{false}; + }; + std::vector qc_dual_recovery; }; } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 5d29e97fa0..fb990e43e4 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -525,6 +525,9 @@ run_barrier(dual_simplex::user_problem_t& user_problem, auto status = dual_simplex::solve_linear_program_with_barrier( user_problem, barrier_settings, timer.get_tic_start(), solution); + if (!user_problem.qc_dual_recovery.empty()) { + detail::project_barrier_qcqp_duals_to_model(user_problem, solution); + } detail::project_barrier_solution_to_model_variables(user_problem, solution); CUOPT_LOG_CONDITIONAL_INFO( @@ -1844,17 +1847,9 @@ optimization_problem_solution_t solve_qcqp( method_t::Barrier); if (has_qc) { - CUOPT_LOG_INFO("Dual variables for problems with quadratic constraints not returned."); - const f_t nan_val = std::numeric_limits::quiet_NaN(); - auto stream = op_problem.get_handle_ptr()->get_stream(); - thrust::fill(rmm::exec_policy(stream), - solution.get_dual_solution().begin(), - solution.get_dual_solution().end(), - nan_val); - thrust::fill(rmm::exec_policy(stream), - solution.get_reduced_cost().begin(), - solution.get_reduced_cost().end(), - nan_val); + CUOPT_LOG_INFO( + "Quadratic-constraint dual multipliers recovered from barrier SOC solution " + "(linear rows + one entry per QCMATRIX row)."); } if (settings.sol_file != "") { diff --git a/cpp/tests/socp/general_quadratic_test.cu b/cpp/tests/socp/general_quadratic_test.cu index 2918d80a7d..ad88c42d8e 100644 --- a/cpp/tests/socp/general_quadratic_test.cu +++ b/cpp/tests/socp/general_quadratic_test.cu @@ -17,6 +17,7 @@ #include #include +#include #include namespace cuopt::linear_programming::detail::test { @@ -24,6 +25,8 @@ namespace cuopt::linear_programming::detail::test { using i_t = int; using f_t = double; using qc_t = optimization_problem_interface_t::quadratic_constraint_t; +using qc_soc_path_t = + cuopt::linear_programming::dual_simplex::user_problem_t::qc_soc_recognition_path_t; static void init_handler(const raft::handle_t* handle_ptr) { @@ -33,6 +36,56 @@ static void init_handler(const raft::handle_t* handle_ptr) handle_ptr->get_cusparse_handle(), CUSPARSE_POINTER_MODE_DEVICE, handle_ptr->get_stream())); } +/** Barrier solve on expanded SOC model, then project duals/primals back to the user QCQP. */ +static void solve_qc_soc_barrier_with_dual_recovery( + cuopt::linear_programming::dual_simplex::user_problem_t& user_problem, + cuopt::linear_programming::dual_simplex::lp_solution_t& solution) +{ + using namespace cuopt::linear_programming::dual_simplex; + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + + const auto status = solve_linear_program_with_barrier(user_problem, settings, solution); + ASSERT_EQ(status, lp_status_t::OPTIMAL); + + project_barrier_qcqp_duals_to_model(user_problem, solution); + project_barrier_solution_to_model_variables(user_problem, solution); +} + +static f_t qc_dual_from_projected_solution( + const cuopt::linear_programming::dual_simplex::user_problem_t& user_problem, + const cuopt::linear_programming::dual_simplex::lp_solution_t& solution, + i_t qc_index = 0) +{ + const i_t m_linear = user_problem.original_num_rows; + return solution.y[m_linear + qc_index]; +} + +static void print_lp_solution_debug( + const char* label, + const cuopt::linear_programming::dual_simplex::lp_solution_t& solution) +{ + std::cout << "=== " << label << " ===\n"; + std::cout << " objective = " << solution.objective << '\n'; + std::cout << " x[" << solution.x.size() << "]:"; + for (size_t i = 0; i < solution.x.size(); ++i) { + std::cout << ' ' << i << ':' << solution.x[i]; + } + std::cout << '\n'; + std::cout << " y[" << solution.y.size() << "]:"; + for (size_t i = 0; i < solution.y.size(); ++i) { + std::cout << ' ' << i << ':' << solution.y[i]; + } + std::cout << '\n'; + std::cout << " z[" << solution.z.size() << "]:"; + for (size_t i = 0; i < solution.z.size(); ++i) { + std::cout << ' ' << i << ':' << solution.z[i]; + } + std::cout << '\n'; +} + // Test: general convex quadratic constraint with dense PD Q matrix. // minimize x0 + x1 // subject to x^T Q x <= 1 where Q = [2 1; 1 2] (PD, eigenvalues 1 and 3) @@ -766,6 +819,14 @@ TEST(general_quadratic, rotated_soc_heads_nonneg_accepted) (convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem))); EXPECT_GT(user_problem.second_order_cone_dims.size(), 0u); + ASSERT_EQ(user_problem.qc_dual_recovery.size(), 1u); + const auto& meta = user_problem.qc_dual_recovery[0]; + using path_t = user_problem_t::qc_soc_recognition_path_t; + EXPECT_EQ(meta.path, path_t::ROTATED); + EXPECT_GE(meta.rsoc_s0_lift_row, 0); + EXPECT_GE(meta.rsoc_s1_lift_row, 0); + EXPECT_GT(meta.rsoc_head_lift_h, 0.0); + EXPECT_FALSE(meta.rsoc_head1_is_constant_half); } // Test: x0^2 + x1^2 - 2*y*z <= 0 with y free, z free should be rejected (non-convex). @@ -824,4 +885,640 @@ TEST(general_quadratic, rotated_soc_heads_free_rejected) cuopt::logic_error); } +// End-to-end barrier dual recovery (mirrors python/cuopt/cuopt/tests/socp/test_socp.py path tests). +// min t s.t. x1=1, x2=0, x1^2+x2^2 <= t^2, t>=0 => t=1, mu_QC=1/2. +TEST(qc_dual_recovery, lorentz_path) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 2; + constexpr int n = 3; + constexpr int nz = 2; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {0.0, 0.0, 1.0}; // min t + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 1, 2}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 1; + user_problem.A.x[1] = 1.0; + + user_problem.rhs = {1.0, 0.0}; + user_problem.row_sense = {'E', 'E'}; + user_problem.lower = {-inf, -inf, 0.0}; + user_problem.upper = {inf, inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "lorentz_soc"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 0.0; + qc.rows = {0, 1, 2}; + qc.cols = {0, 1, 2}; + qc.vals = {1.0, 1.0, -1.0}; + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 1, 2}; + csr_A.j = {0, 1}; + csr_A.x = {1.0, 1.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + csr_A.to_compressed_col(user_problem.A); + + ASSERT_EQ(user_problem.qc_dual_recovery.size(), 1u); + EXPECT_EQ(user_problem.qc_dual_recovery[0].path, qc_soc_path_t::LORENTZ); + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + solve_qc_soc_barrier_with_dual_recovery(user_problem, solution); + + EXPECT_NEAR(solution.objective, 1.0, 1e-4); + EXPECT_NEAR(solution.x[2], 1.0, 1e-4); // t + EXPECT_NEAR(qc_dual_from_projected_solution(user_problem, solution), 0.5, 1e-4); +} + +// min t s.t. x1=0, x2=0, x1^2+x2^2 <= t^2, t>=0 => apex; mu_QC exported as 0. +TEST(qc_dual_recovery, lorentz_degenerate) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 2; + constexpr int n = 3; + constexpr int nz = 2; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {0.0, 0.0, 1.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 1, 2}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 1; + user_problem.A.x[1] = 1.0; + + user_problem.rhs = {0.0, 0.0}; + user_problem.row_sense = {'E', 'E'}; + user_problem.lower = {-inf, -inf, 0.0}; + user_problem.upper = {inf, inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "lorentz_soc"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 0.0; + qc.rows = {0, 1, 2}; + qc.cols = {0, 1, 2}; + qc.vals = {1.0, 1.0, -1.0}; + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 1, 2}; + csr_A.j = {0, 1}; + csr_A.x = {1.0, 1.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + csr_A.to_compressed_col(user_problem.A); + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + solve_qc_soc_barrier_with_dual_recovery(user_problem, solution); + + EXPECT_NEAR(solution.objective, 0.0, 1e-4); + EXPECT_NEAR(solution.x[0], 0.0, 1e-4); + EXPECT_NEAR(solution.x[1], 0.0, 1e-4); + EXPECT_NEAR(solution.x[2], 0.0, 1e-4); + EXPECT_NEAR(qc_dual_from_projected_solution(user_problem, solution), 0.0, 1e-4); +} + +// min x1 s.t. x0=1, x1^2+x2^2 <= x0 => x1=-1, mu_QC=1/2. +TEST(qc_dual_recovery, affine_path) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 1; + constexpr int n = 3; // x0, x1, x2 + constexpr int nz = 1; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {0.0, 1.0, 0.0}; // min x1 + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 1, 1, 1}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + + user_problem.rhs = {1.0}; + user_problem.row_sense = {'E'}; + user_problem.lower = {-inf, -inf, -inf}; + user_problem.upper = {inf, inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "affine_soc"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 0.0; + qc.rows = {1, 2}; // x1^2 + x2^2 + qc.cols = {1, 2}; + qc.vals = {1.0, 1.0}; + qc.linear_indices = {0}; // -x0 + qc.linear_values = {-1.0}; + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 1}; + csr_A.j = {0}; + csr_A.x = {1.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + csr_A.to_compressed_col(user_problem.A); + + ASSERT_EQ(user_problem.qc_dual_recovery.size(), 1u); + EXPECT_EQ(user_problem.qc_dual_recovery[0].path, qc_soc_path_t::AFFINE); + EXPECT_GE(user_problem.qc_dual_recovery[0].affine_link_row, 0); + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + solve_qc_soc_barrier_with_dual_recovery(user_problem, solution); + + EXPECT_NEAR(solution.objective, -1.0, 1e-4); + EXPECT_NEAR(solution.x[0], 1.0, 1e-4); // x0 + EXPECT_NEAR(solution.x[1], -1.0, 1e-4); // x1 + EXPECT_NEAR(solution.x[2], 0.0, 1e-4); // x2 + EXPECT_NEAR(solution.y[0], -0.5, 1e-4); // EQ x0 = 1 + EXPECT_NEAR(qc_dual_from_projected_solution(user_problem, solution), 0.5, 1e-4); +} + +// min x0+x1 s.t. x0=0, x1=0, x2^2+x3^2 <= x0*x1, x0,x1>=0 => lifted s0=s1=0; mu_QC=0. +TEST(qc_dual_recovery, rotated_degenerate) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 2; + constexpr int n = 4; + constexpr int nz = 2; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {1.0, 1.0, 0.0, 0.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 1, 2, 2, 2}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 1; + user_problem.A.x[1] = 1.0; + + user_problem.rhs = {0.0, 0.0}; + user_problem.row_sense = {'E', 'E'}; + user_problem.lower = {0.0, 0.0, -inf, -inf}; + user_problem.upper = {inf, inf, inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "rotated_soc"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 0.0; + qc.rows = {2, 3, 0, 1}; + qc.cols = {2, 3, 1, 0}; + qc.vals = {1.0, 1.0, -0.5, -0.5}; + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 1, 2}; + csr_A.j = {0, 1}; + csr_A.x = {1.0, 1.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + csr_A.to_compressed_col(user_problem.A); + + ASSERT_EQ(user_problem.qc_dual_recovery.size(), 1u); + EXPECT_EQ(user_problem.qc_dual_recovery[0].path, qc_soc_path_t::ROTATED); + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + solve_qc_soc_barrier_with_dual_recovery(user_problem, solution); + + EXPECT_NEAR(solution.objective, 0.0, 1e-4); + EXPECT_NEAR(solution.x[0], 0.0, 1e-4); + EXPECT_NEAR(solution.x[1], 0.0, 1e-4); + EXPECT_NEAR(solution.x[2], 0.0, 1e-4); + EXPECT_NEAR(solution.x[3], 0.0, 1e-4); + EXPECT_NEAR(qc_dual_from_projected_solution(user_problem, solution), 0.0, 1e-4); +} + +// min x0+3*x2 s.t. x1=2, x3=0, x2^2+x3^2 <= x0*x1, x0,x1>=0 => obj=-9/2, mu_QC=1/2. +TEST(qc_dual_recovery, rotated_path) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 2; + constexpr int n = 4; + constexpr int nz = 2; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {1.0, 0.0, 3.0, 0.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 0, 1, 1, 1}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 1; + user_problem.A.x[1] = 1.0; + + user_problem.rhs = {2.0, 0.0}; + user_problem.row_sense = {'E', 'E'}; + user_problem.lower = {0.0, 0.0, -inf, -inf}; + user_problem.upper = {inf, inf, inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "rotated_soc"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 0.0; + qc.rows = {2, 3, 0, 1}; + qc.cols = {2, 3, 1, 0}; + qc.vals = {1.0, 1.0, -0.5, -0.5}; + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 1, 2}; + csr_A.j = {1, 3}; + csr_A.x = {1.0, 1.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + csr_A.to_compressed_col(user_problem.A); + + ASSERT_EQ(user_problem.qc_dual_recovery.size(), 1u); + EXPECT_EQ(user_problem.qc_dual_recovery[0].path, qc_soc_path_t::ROTATED); + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + solve_qc_soc_barrier_with_dual_recovery(user_problem, solution); + + EXPECT_NEAR(solution.objective, -4.5, 1e-4); + EXPECT_NEAR(solution.x[0], 4.5, 1e-4); + EXPECT_NEAR(solution.x[1], 2.0, 1e-4); + EXPECT_NEAR(solution.x[2], -3.0, 1e-4); + EXPECT_NEAR(solution.x[3], 0.0, 1e-4); + EXPECT_NEAR(qc_dual_from_projected_solution(user_problem, solution), 0.5, 1e-4); +} + +// min 2*x0+x1 s.t. 2*x0^2+2*x0*x1+2*x1^2<=1, x0-3*x1=0 +// => x0=-3/sqrt(26), x1=-1/sqrt(26), mu_QC=7/(2*sqrt(26)), y_KKT=-3/26. +TEST(qc_dual_recovery, general_path) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 1; + constexpr int n = 2; + constexpr int nz = 2; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {2.0, 1.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 1, 2}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 0; + user_problem.A.x[1] = -3.0; + + user_problem.rhs = {0.0}; + user_problem.row_sense = {'E'}; + user_problem.lower = {-inf, -inf}; + user_problem.upper = {inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "general_qc"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 1.0; + qc.rows = {0, 1, 1}; + qc.cols = {0, 0, 1}; + qc.vals = {2.0, 2.0, 2.0}; // 2*x0^2 + 2*x0*x1 + 2*x1^2 + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 2}; + csr_A.j = {0, 1}; + csr_A.x = {1.0, -3.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + csr_A.to_compressed_col(user_problem.A); + + ASSERT_EQ(user_problem.qc_dual_recovery.size(), 1u); + const auto& meta = user_problem.qc_dual_recovery[0]; + EXPECT_EQ(meta.path, qc_soc_path_t::GENERAL); + EXPECT_GE(meta.s0_link_row, 0); + EXPECT_GE(meta.sr1_link_row, 0); + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + ASSERT_EQ(solve_linear_program_with_barrier(user_problem, settings, solution), + lp_status_t::OPTIMAL); + + std::cout << "qc_dual_recovery.general_path debug:\n"; + std::cout << " expanded num_rows=" << user_problem.num_rows + << " num_cols=" << user_problem.num_cols + << " original_num_rows=" << user_problem.original_num_rows + << " original_num_cols=" << user_problem.original_num_cols + << " cone_var_start=" << user_problem.cone_var_start << '\n'; + std::cout << " s0_link_row=" << meta.s0_link_row << " sr1_link_row=" << meta.sr1_link_row + << " cone_index=" << meta.cone_index << '\n'; + if (!user_problem.second_order_cone_dims.empty()) { + std::cout << " cone_dim=" << user_problem.second_order_cone_dims[0] << '\n'; + } + print_lp_solution_debug("expanded barrier (pre-projection)", solution); + + project_barrier_qcqp_duals_to_model(user_problem, solution); + print_lp_solution_debug("after project_barrier_qcqp_duals_to_model", solution); + + project_barrier_solution_to_model_variables(user_problem, solution); + print_lp_solution_debug("user model (final)", solution); + + const f_t inv_sqrt26 = 1.0 / std::sqrt(26.0); + EXPECT_NEAR(solution.objective, -7.0 * inv_sqrt26, 1e-4); + EXPECT_NEAR(solution.x[0], -3.0 * inv_sqrt26, 1e-4); + EXPECT_NEAR(solution.x[1], -inv_sqrt26, 1e-4); + EXPECT_GT(std::abs(solution.x[0] - solution.x[1]), 1e-3); + EXPECT_NEAR(solution.y[0], 3.0 / 26.0, 1e-4); // EQ x0 - 3*x1 = 0 (y_KKT = -3/26) + EXPECT_NEAR( + qc_dual_from_projected_solution(user_problem, solution), 7.0 * inv_sqrt26 / 2.0, 1e-4); +} + +// min x0+2*x1 s.t. (x0+x1)^2<=1, x0-2*x1=0 (rank-1 PSD QC, LDLT rank r=1) +// => x0=-2/3, x1=-1/3, obj=-4/3, mu_QC=2/3, y_KKT=1/3. +TEST(qc_dual_recovery, general_path_rank_one) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 1; + constexpr int n = 2; + constexpr int nz = 2; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {1.0, 2.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 1, 2}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 0; + user_problem.A.x[1] = -2.0; + + user_problem.rhs = {0.0}; + user_problem.row_sense = {'E'}; + user_problem.lower = {-inf, -inf}; + user_problem.upper = {inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + // (x0+x1)^2 <= 1 via Q COO (0,0,1), (1,0,2), (1,1,1) -> H = [2,2;2,2], rank 1. + // Cross term needs coefficient 2 (not 1): vals {1,1,1} would give x0^2+x0*x1+x1^2, full rank. + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "rank1_general_qc"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 1.0; + qc.rows = {0, 1, 1}; + qc.cols = {0, 0, 1}; + qc.vals = {1.0, 2.0, 1.0}; + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 2}; + csr_A.j = {0, 1}; + csr_A.x = {1.0, -2.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + csr_A.to_compressed_col(user_problem.A); + + ASSERT_EQ(user_problem.qc_dual_recovery.size(), 1u); + const auto& meta = user_problem.qc_dual_recovery[0]; + EXPECT_EQ(meta.path, qc_soc_path_t::GENERAL); + EXPECT_GE(meta.s0_link_row, 0); + EXPECT_GE(meta.sr1_link_row, 0); + ASSERT_EQ(user_problem.second_order_cone_dims.size(), 1u); + EXPECT_EQ(user_problem.second_order_cone_dims[0], 3); // rank 1 + 2 slacks + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + ASSERT_EQ(solve_linear_program_with_barrier(user_problem, settings, solution), + lp_status_t::OPTIMAL); + + std::cout << "qc_dual_recovery.general_path_rank_one debug:\n"; + std::cout << " expanded num_rows=" << user_problem.num_rows + << " num_cols=" << user_problem.num_cols + << " original_num_rows=" << user_problem.original_num_rows + << " original_num_cols=" << user_problem.original_num_cols + << " cone_var_start=" << user_problem.cone_var_start << '\n'; + std::cout << " s0_link_row=" << meta.s0_link_row << " sr1_link_row=" << meta.sr1_link_row + << " cone_index=" << meta.cone_index << '\n'; + if (!user_problem.second_order_cone_dims.empty()) { + std::cout << " cone_dim=" << user_problem.second_order_cone_dims[0] << '\n'; + } + print_lp_solution_debug("expanded barrier (pre-projection)", solution); + + project_barrier_qcqp_duals_to_model(user_problem, solution); + print_lp_solution_debug("after project_barrier_qcqp_duals_to_model", solution); + + project_barrier_solution_to_model_variables(user_problem, solution); + print_lp_solution_debug("user model (final)", solution); + + EXPECT_NEAR(solution.objective, -4.0 / 3.0, 1e-4); + EXPECT_NEAR(solution.x[0], -2.0 / 3.0, 1e-4); + EXPECT_NEAR(solution.x[1], -1.0 / 3.0, 1e-4); + EXPECT_NEAR(solution.y[0], -1.0 / 3.0, 1e-4); // EQ x0 - 2*x1 = 0 (y_KKT = 1/3) + EXPECT_NEAR(qc_dual_from_projected_solution(user_problem, solution), 2.0 / 3.0, 1e-4); +} + +// min x0+x1 s.t. 2*x0^2+2*x0*x1+2*x1^2+x0+x1<=3, x0-x1=0 (QC linear part c!=0) +// => x0=x1=(-1-sqrt(19))/6, obj=(-1-sqrt(19))/3, mu_QC=1/sqrt(19), y_KKT=0. +TEST(qc_dual_recovery, general_path_linear_qc) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 1; + constexpr int n = 2; + constexpr int nz = 2; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {1.0, 1.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 1, 2}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 0; + user_problem.A.x[1] = -1.0; + + user_problem.rhs = {0.0}; + user_problem.row_sense = {'E'}; + user_problem.lower = {-inf, -inf}; + user_problem.upper = {inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "linear_general_qc"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 3.0; + qc.rows = {0, 1, 1}; + qc.cols = {0, 0, 1}; + qc.vals = {2.0, 2.0, 2.0}; + qc.linear_indices = {0, 1}; + qc.linear_values = {1.0, 1.0}; // + x0 + x1 + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 2}; + csr_A.j = {0, 1}; + csr_A.x = {1.0, -1.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + csr_A.to_compressed_col(user_problem.A); + + ASSERT_EQ(user_problem.qc_dual_recovery.size(), 1u); + const auto& meta = user_problem.qc_dual_recovery[0]; + EXPECT_EQ(meta.path, qc_soc_path_t::GENERAL); + EXPECT_GE(meta.s0_link_row, 0); + EXPECT_GE(meta.sr1_link_row, 0); + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + ASSERT_EQ(solve_linear_program_with_barrier(user_problem, settings, solution), + lp_status_t::OPTIMAL); + + std::cout << "qc_dual_recovery.general_path_linear_qc debug:\n"; + std::cout << " expanded num_rows=" << user_problem.num_rows + << " num_cols=" << user_problem.num_cols + << " original_num_rows=" << user_problem.original_num_rows + << " original_num_cols=" << user_problem.original_num_cols + << " cone_var_start=" << user_problem.cone_var_start << '\n'; + std::cout << " s0_link_row=" << meta.s0_link_row << " sr1_link_row=" << meta.sr1_link_row + << " cone_index=" << meta.cone_index << '\n'; + if (!user_problem.second_order_cone_dims.empty()) { + std::cout << " cone_dim=" << user_problem.second_order_cone_dims[0] << '\n'; + } + print_lp_solution_debug("expanded barrier (pre-projection)", solution); + + project_barrier_qcqp_duals_to_model(user_problem, solution); + print_lp_solution_debug("after project_barrier_qcqp_duals_to_model", solution); + + project_barrier_solution_to_model_variables(user_problem, solution); + print_lp_solution_debug("user model (final)", solution); + + const f_t sqrt19 = std::sqrt(19.0); + const f_t t_star = (-1.0 - sqrt19) / 6.0; + const f_t mu_qc = 1.0 / sqrt19; + + EXPECT_NEAR(solution.objective, 2.0 * t_star, 1e-4); + EXPECT_NEAR(solution.x[0], t_star, 1e-4); + EXPECT_NEAR(solution.x[1], t_star, 1e-4); + EXPECT_NEAR(solution.y[0], 0.0, 1e-4); // EQ x0 - x1 = 0 + EXPECT_NEAR(qc_dual_from_projected_solution(user_problem, solution), mu_qc, 1e-4); +} + } // namespace cuopt::linear_programming::detail::test diff --git a/python/cuopt/cuopt/linear_programming/problem.py b/python/cuopt/cuopt/linear_programming/problem.py index 10600a543b..e153aeb569 100644 --- a/python/cuopt/cuopt/linear_programming/problem.py +++ b/python/cuopt/cuopt/linear_programming/problem.py @@ -2186,14 +2186,22 @@ def populate_solution(self, solution): dual_sol = None if not IsMIP: dual_sol = solution.get_dual_solution() + n_linear = sum(1 for c in self.constrs if not c.is_quadratic) + qc_dual_idx = 0 linear_row = 0 for constr in self.constrs: if constr.is_quadratic: - continue - if dual_sol is not None and len(dual_sol) > linear_row: + if ( + dual_sol is not None + and len(dual_sol) > n_linear + qc_dual_idx + ): + constr.DualValue = dual_sol[n_linear + qc_dual_idx] + qc_dual_idx += 1 + elif dual_sol is not None and len(dual_sol) > linear_row: constr.DualValue = dual_sol[linear_row] constr.Slack = constr.compute_slack() - linear_row += 1 + if not constr.is_quadratic: + linear_row += 1 self.solved = True def solve(self, settings=solver_settings.SolverSettings()): diff --git a/python/cuopt/cuopt/tests/socp/test_socp.py b/python/cuopt/cuopt/tests/socp/test_socp.py index 87d8f3b5c9..c3306847e8 100644 --- a/python/cuopt/cuopt/tests/socp/test_socp.py +++ b/python/cuopt/cuopt/tests/socp/test_socp.py @@ -2,11 +2,11 @@ # SPDX-License-Identifier: Apache-2.0 """ -Barrier SOCP tests via the Problem Python API. +Barrier SOCP / QCQP tests via the Problem Python API. -Checks that the barrier solution is mapped back to the original model variables -after SOC conversion (see ``project_barrier_solution_to_model_variables`` in -``cpp/src/barrier/translate_soc.hpp``). +Python ``dual_api_`` tests exercise each SOC translation path through the Problem +API and verify dual extraction only: ``get_dual_solution()`` length/layout and +``Constraint.DualValue`` agree for linear and quadratic rows. """ from __future__ import annotations @@ -31,6 +31,7 @@ OBJ_TOL = 1e-6 PRIMAL_TOL = 1e-6 FEAS_TOL = 1e-6 +DUAL_TOL = 1e-5 def _barrier_settings() -> SolverSettings: @@ -39,8 +40,278 @@ def _barrier_settings() -> SolverSettings: return settings +def _quadratic_constraint_violation(constr, variables) -> float: + """QCMATRIX row value minus rhs (feasible when <= 0 for L rows).""" + vals = [var.Value for var in variables] + quad = 0.0 + for k in range(len(constr.vals)): + i = int(constr.rows[k]) + j = int(constr.cols[k]) + quad += float(constr.vals[k]) * vals[i] * vals[j] + lin = 0.0 + for k in range(len(constr.linear_values)): + lin += ( + float(constr.linear_values[k]) + * vals[int(constr.linear_indices[k])] + ) + return quad + lin - float(constr.rhs_value) + + +def _assert_feasible(problem: Problem) -> None: + variables = problem.getVariables() + for constr in problem.getConstraints(): + if constr.is_quadratic: + assert ( + _quadratic_constraint_violation(constr, variables) <= FEAS_TOL + ) + continue + slack = constr.compute_slack() + if constr.Sense == LE: + assert slack >= -FEAS_TOL + elif constr.Sense == GE: + assert slack <= FEAS_TOL + else: + assert constr.Sense == EQ + assert slack == pytest.approx(0.0, abs=FEAS_TOL) + + +def _assert_solution_on_original_model(problem: Problem, solution) -> None: + primal = solution.get_primal_solution() + assert len(primal) == problem.NumVariables + assert problem.ObjValue == pytest.approx( + solution.get_primal_objective(), rel=0, abs=OBJ_TOL + ) + assert problem.ObjValue == pytest.approx( + problem.getObjective().getValue(), rel=0, abs=OBJ_TOL + ) + + +def _assert_dual_layout(problem: Problem, solution) -> None: + """``get_dual_solution()`` matches per-constraint ``DualValue`` on the user model.""" + dual = solution.get_dual_solution() + n_linear = sum(1 for c in problem.getConstraints() if not c.is_quadratic) + n_qc = sum(1 for c in problem.getConstraints() if c.is_quadratic) + assert len(dual) == n_linear + n_qc + + lin_idx = 0 + qc_idx = 0 + for constr in problem.getConstraints(): + if constr.is_quadratic: + assert constr.DualValue == pytest.approx( + dual[n_linear + qc_idx], rel=0, abs=DUAL_TOL + ) + qc_idx += 1 + else: + assert constr.DualValue == pytest.approx( + dual[lin_idx], rel=0, abs=DUAL_TOL + ) + lin_idx += 1 + + +def _assert_duals_finite(problem: Problem) -> None: + for constr in problem.getConstraints(): + assert np.isfinite(constr.DualValue) + + +def _solve_barrier(problem: Problem): + solution = problem.solve(_barrier_settings()) + assert problem.Status.name == "Optimal" + _assert_solution_on_original_model(problem, solution) + _assert_feasible(problem) + return solution + + +def _solve_and_assert_dual_api(problem: Problem) -> None: + solution = _solve_barrier(problem) + _assert_dual_layout(problem, solution) + _assert_duals_finite(problem) + assert any(c.is_quadratic for c in problem.getConstraints()) + + +# --------------------------------------------------------------------------- +# Path-specific builders +# --------------------------------------------------------------------------- + + +def build_lorentz_path_problem() -> tuple[Problem, tuple]: + """ + LORENTZ translation path (mirrors C++ ``qc_dual_recovery/lorentz_path``). + + min t + s.t. x1 = 1, x2 = 0, x1^2 + x2^2 <= t^2, t >= 0. + + Optimum: t = 1, obj = 1, mu_QC = 1/2. + """ + problem = Problem("path_lorentz") + x1 = problem.addVariable(lb=-np.inf, name="x1") + x2 = problem.addVariable(lb=-np.inf, name="x2") + t = problem.addVariable(lb=0, name="t") + problem.setObjective(t) + problem.addConstraint(x1 == 1) + problem.addConstraint(x2 == 0) + problem.addConstraint(x1 * x1 + x2 * x2 - t * t <= 0, name="lorentz_soc") + return problem, (x1, x2, t) + + +def build_lorentz_degenerate_problem() -> tuple[Problem, tuple]: + """ + LORENTZ apex case (mirrors C++ ``qc_dual_recovery/lorentz_degenerate``). + + min t + s.t. x1 = 0, x2 = 0, x1^2 + x2^2 <= t^2, t >= 0. + + Optimum: x1 = x2 = t = 0, obj = 0 at the cone apex. grad g = 0 and the QC + multiplier is not unique; recovery exports mu_QC = 0 by convention. + """ + problem = Problem("path_lorentz_degenerate") + x1 = problem.addVariable(lb=-np.inf, name="x1") + x2 = problem.addVariable(lb=-np.inf, name="x2") + t = problem.addVariable(lb=0, name="t") + problem.setObjective(t) + problem.addConstraint(x1 == 0) + problem.addConstraint(x2 == 0) + problem.addConstraint(x1 * x1 + x2 * x2 - t * t <= 0, name="lorentz_soc") + return problem, (x1, x2, t) + + +def build_affine_path_problem() -> tuple[Problem, tuple]: + """ + AFFINE translation path (mirrors C++ ``qc_dual_recovery/affine_path``). + + min x2 + s.t. x0 = 1, x2^2 + x3^2 <= x0. + + Optimum: x0 = 1, x2 = -1, x3 = 0, obj = -1, mu_QC = 1/2. + (C++ names the objective tail x1 and uses x1^2 + x2^2 <= x0; same model.) + """ + problem = Problem("path_affine") + x0 = problem.addVariable(lb=-np.inf, name="x0") + x2 = problem.addVariable(lb=-np.inf, name="x2") + x3 = problem.addVariable(lb=-np.inf, name="x3") + problem.setObjective(x2) + problem.addConstraint(x0 == 1) + problem.addConstraint(x2 * x2 + x3 * x3 - x0 <= 0, name="affine_soc") + return problem, (x0, x2, x3) + + +def build_rotated_path_problem() -> tuple[Problem, tuple]: + """ + ROTATED translation path with free heads in the objective + (mirrors C++ ``qc_dual_recovery/rotated_path``). + + min x0 + 3*x2 + s.t. x1 = 2, x3 = 0, x2^2 + x3^2 <= x0*x1, x0, x1 >= 0. + + Active QC: x2^2 = 2*x0 => x2 = -sqrt(2*x0). Minimizing + f(x0) = x0 - 3*sqrt(2*x0) gives x0 = 9/2, x2 = -3, obj = -9/2. + + KKT (lambda on g = x2^2 + x3^2 - x0*x1 <= 0): + x0: 1 - lambda*x1 = 0 => lambda = 1/2 + x2: 3 + lambda*2*x2 = 0 + EQ x1: DualValue = -lambda*x0 = -9/4; EQ x3: 0. + + At the optimum the lifted slacks are s0 = h*(x0+x1) > 0 and s1 = h*(x0-x1) > 0 + with h = 1/2, so recovery uses the rotated path with nonzero s1. + """ + problem = Problem("path_rotated") + x0 = problem.addVariable(lb=0, name="x0") + x1 = problem.addVariable(lb=0, name="x1") + x2 = problem.addVariable(lb=-np.inf, name="x2") + x3 = problem.addVariable(lb=-np.inf, name="x3") + problem.setObjective(x0 + 3 * x2) + problem.addConstraint(x1 == 2) + problem.addConstraint(x3 == 0) + problem.addConstraint( + x2 * x2 + x3 * x3 - 0.5 * x0 * x1 - 0.5 * x1 * x0 <= 0, + name="rotated_soc", + ) + return problem, (x0, x1, x2, x3) + + +def build_rotated_degenerate_problem() -> tuple[Problem, tuple]: + """ + ROTATED apex case (mirrors C++ ``qc_dual_recovery/rotated_degenerate``). + + min x0 + x1 + s.t. x0 = 0, x1 = 0, x2^2 + x3^2 <= x0*x1, x0, x1 >= 0. + + Optimum: all variables zero, obj = 0 at the lifted cone apex. Both rotated + heads are pinned to zero (s0 = s1 = 0), grad g = 0, and the QC multiplier + is not unique; recovery exports mu_QC = 0 by convention. + """ + problem = Problem("path_rotated_degenerate") + x0 = problem.addVariable(lb=0, name="x0") + x1 = problem.addVariable(lb=0, name="x1") + x2 = problem.addVariable(lb=-np.inf, name="x2") + x3 = problem.addVariable(lb=-np.inf, name="x3") + problem.setObjective(x0 + x1) + problem.addConstraint(x0 == 0) + problem.addConstraint(x1 == 0) + problem.addConstraint( + x2 * x2 + x3 * x3 - 0.5 * x0 * x1 - 0.5 * x1 * x0 <= 0, + name="rotated_soc", + ) + return problem, (x0, x1, x2, x3) + + +def build_general_path_problem() -> tuple[Problem, tuple]: + """ + GENERAL (LDLT) translation path, full-rank PSD QC + (mirrors C++ ``qc_dual_recovery/general_path``). + + min 2*x0 + x1 + s.t. 2*x0^2 + 2*x0*x1 + 2*x1^2 <= 1, x0 - 3*x1 = 0. + + Optimum: x0 = -3/sqrt(26), x1 = -1/sqrt(26), obj = -7/sqrt(26), + mu_QC = 7/(2*sqrt(26)), y_KKT = -3/26 on the equality row. + """ + problem = Problem("path_general") + x0 = problem.addVariable(lb=-np.inf, name="x0") + x1 = problem.addVariable(lb=-np.inf, name="x1") + problem.setObjective(2 * x0 + x1) + problem.addConstraint( + 2 * x0 * x0 + 2 * x0 * x1 + 2 * x1 * x1 - 1 <= 0, name="general_qc" + ) + problem.addConstraint(x0 - 3 * x1 == 0) + return problem, (x0, x1) + + +def test_socp_dual_api_lorentz_path(): + problem, _ = build_lorentz_path_problem() + _solve_and_assert_dual_api(problem) + + +def test_socp_dual_api_lorentz_degenerate(): + problem, _ = build_lorentz_degenerate_problem() + _solve_and_assert_dual_api(problem) + + +def test_socp_dual_api_affine_path(): + problem, _ = build_affine_path_problem() + _solve_and_assert_dual_api(problem) + + +def test_socp_dual_api_rotated_degenerate(): + problem, _ = build_rotated_degenerate_problem() + _solve_and_assert_dual_api(problem) + + +def test_socp_dual_api_rotated_path(): + problem, _ = build_rotated_path_problem() + _solve_and_assert_dual_api(problem) + + +def test_socp_dual_api_general_path(): + problem, _ = build_general_path_problem() + _solve_and_assert_dual_api(problem) + + +# --------------------------------------------------------------------------- +# Larger Lorentz-style regression models (primal only) +# --------------------------------------------------------------------------- + + def _soc_two_dim_constraint(problem, x0, x1, mat, head) -> None: - """Encode ||mat @ [x0, x1]||_2 <= head as a standard Lorentz cone in (head, z0, z1).""" z0 = problem.addVariable(lb=-np.inf) z1 = problem.addVariable(lb=-np.inf) problem.addConstraint(z0 == mat[0, 0] * x0 + mat[0, 1] * x1) @@ -49,7 +320,6 @@ def _soc_two_dim_constraint(problem, x0, x1, mat, head) -> None: def build_socp_1() -> tuple[Problem, tuple]: - """Min 3*x0+2*x1+x2 s.t. ||x||_2 <= y, x0+x1+3*x2 >= 1, 0 <= y <= 5.""" problem = Problem("socp_1") x0 = problem.addVariable(lb=-np.inf, name="x0") x1 = problem.addVariable(lb=-np.inf, name="x1") @@ -64,7 +334,6 @@ def build_socp_1() -> tuple[Problem, tuple]: def build_socp_3() -> tuple[Problem, tuple]: - """Min -x0+2*x1 s.t. ||M_i x||_2 <= 1 for three fixed 2x2 maps M_i.""" root2 = np.sqrt(2.0) u = np.array([[1 / root2, -1 / root2], [1 / root2, 1 / root2]]) mat1 = np.diag([root2, 1 / root2]) @ u.T @@ -84,64 +353,9 @@ def build_socp_3() -> tuple[Problem, tuple]: return problem, (x0, x1, h1, h2, h3) -def _quadratic_constraint_violation(constr, variables) -> float: - """QCMATRIX row value minus rhs (should be <= 0 for L rows).""" - vals = [var.Value for var in variables] - quad = 0.0 - for k in range(len(constr.vals)): - i = int(constr.rows[k]) - j = int(constr.cols[k]) - quad += float(constr.vals[k]) * vals[i] * vals[j] - lin = 0.0 - for k in range(len(constr.linear_values)): - lin += ( - float(constr.linear_values[k]) - * vals[int(constr.linear_indices[k])] - ) - return quad + lin - float(constr.rhs_value) - - -def _assert_solution_on_original_model(problem: Problem, solution) -> None: - primal = solution.get_primal_solution() - assert len(primal) == problem.NumVariables - assert problem.ObjValue == pytest.approx( - solution.get_primal_objective(), rel=0, abs=OBJ_TOL - ) - assert problem.ObjValue == pytest.approx( - problem.getObjective().getValue(), rel=0, abs=OBJ_TOL - ) - - -def _assert_feasible(problem: Problem) -> None: - variables = problem.getVariables() - for constr in problem.getConstraints(): - if constr.is_quadratic: - assert ( - _quadratic_constraint_violation(constr, variables) <= FEAS_TOL - ) - continue - slack = constr.compute_slack() - if constr.Sense == LE: - assert slack >= -FEAS_TOL - elif constr.Sense == GE: - assert slack <= FEAS_TOL - else: - assert constr.Sense == EQ - assert slack == pytest.approx(0.0, abs=FEAS_TOL) - - -def _solve(problem: Problem): - solution = problem.solve(_barrier_settings()) - assert problem.Status.name == "Optimal" - return solution - - def test_socp_1_barrier_solution(): problem, (x0, x1, x2, y) = build_socp_1() - solution = _solve(problem) - _assert_solution_on_original_model(problem, solution) - _assert_feasible(problem) - + _solve_barrier(problem) assert problem.ObjValue == pytest.approx( EXPECTED_SOCP_1_OBJECTIVE, abs=OBJ_TOL ) @@ -153,10 +367,7 @@ def test_socp_1_barrier_solution(): def test_socp_3_barrier_solution(): problem, (x0, x1, h1, h2, h3) = build_socp_3() - solution = _solve(problem) - _assert_solution_on_original_model(problem, solution) - _assert_feasible(problem) - + _solve_barrier(problem) assert problem.ObjValue == pytest.approx( EXPECTED_SOCP_3_OBJECTIVE, abs=OBJ_TOL ) @@ -167,79 +378,27 @@ def test_socp_3_barrier_solution(): assert h3.Value == pytest.approx(1.0, abs=PRIMAL_TOL) -def test_general_quadratic_unsymmetric(): - """ - Min x0 + x1 - s.t. 2*x0^2 + 3*x0*x1 + 2*x1^2 <= 1 (unsymmetric Q: cross term only as x0*x1) - x0 - x1 = 0 - - Q is given unsymmetrically: the 3*x0*x1 term is stored as a single - entry (row=0, col=1, val=3) rather than symmetric (0,1,1.5)+(1,0,1.5). - After symmetrization H = [4 3; 3 4], eigenvalues 1 and 7 (PD). - - With x0 = x1 = t: 2t^2 + 3t^2 + 2t^2 = 7t^2 <= 1 - min 2t at t = -1/sqrt(7), obj = -2/sqrt(7) ≈ -0.755929 - """ - problem = Problem("general_qc_unsymmetric") - x0 = problem.addVariable(lb=-np.inf, name="x0") - x1 = problem.addVariable(lb=-np.inf, name="x1") - problem.setObjective(x0 + x1) - problem.addConstraint(2 * x0 * x0 + 3 * x0 * x1 + 2 * x1 * x1 <= 1) - problem.addConstraint(x0 - x1 == 0) - - solution = _solve(problem) - _assert_solution_on_original_model(problem, solution) - _assert_feasible(problem) - - expected_obj = -2.0 / np.sqrt(7.0) - expected_x = -1.0 / np.sqrt(7.0) - assert problem.ObjValue == pytest.approx(expected_obj, abs=OBJ_TOL) - assert x0.Value == pytest.approx(expected_x, abs=PRIMAL_TOL) - assert x1.Value == pytest.approx(expected_x, abs=PRIMAL_TOL) - - def test_maximize_with_quadratic_constraint(): - """ - Maximize x + y - s.t. x + y <= 10 - 2*x^2 + 2*x*y + 2*y^2 <= 6 - - The quadratic constraint is the binding one. - With x = y = t: 2t^2 + 2t^2 + 2t^2 = 6t^2 <= 6 => t in [-1, 1]. - Maximizing 2t gives t = 1, obj = 2. - - Minimizing gives t = -1, obj = -2. - - This test verifies that MAXIMIZE is respected when quadratic constraints - are present (regression for a bug where the QCQP path ignored the - objective sense). - """ from cuopt.linear_programming.problem import MAXIMIZE, MINIMIZE - # Solve as MINIMIZE first to establish baseline prob_min = Problem("qc_maximize_min") x = prob_min.addVariable(lb=-np.inf, name="x") y = prob_min.addVariable(lb=-np.inf, name="y") prob_min.addConstraint(x + y <= 10) prob_min.addConstraint(2 * x * x + 2 * x * y + 2 * y * y <= 6) prob_min.setObjective(x + y, sense=MINIMIZE) - _solve(prob_min) - _assert_feasible(prob_min) - + _solve_barrier(prob_min) assert prob_min.ObjValue == pytest.approx(-2.0, abs=OBJ_TOL) assert x.Value == pytest.approx(-1.0, abs=PRIMAL_TOL) assert y.Value == pytest.approx(-1.0, abs=PRIMAL_TOL) - # Solve as MAXIMIZE - should give the opposite optimum prob_max = Problem("qc_maximize_max") x = prob_max.addVariable(lb=-np.inf, name="x") y = prob_max.addVariable(lb=-np.inf, name="y") prob_max.addConstraint(x + y <= 10) prob_max.addConstraint(2 * x * x + 2 * x * y + 2 * y * y <= 6) prob_max.setObjective(x + y, sense=MAXIMIZE) - _solve(prob_max) - _assert_feasible(prob_max) - + _solve_barrier(prob_max) assert prob_max.ObjValue == pytest.approx(2.0, abs=OBJ_TOL) assert x.Value == pytest.approx(1.0, abs=PRIMAL_TOL) assert y.Value == pytest.approx(1.0, abs=PRIMAL_TOL)