diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 47d879ff42..d2731999ac 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -851,6 +851,67 @@ i_t remove_fixed_variables(f_t fixed_tolerance, return 0; } +// Adds one slack column s_k per entry of `rows`; slack k has a single nonzero in row i = rows[k]: +// "<=" row (is_range = false): a_i^T x + s_k = rhs_i, s_k >= 0 +// range row (is_range = true): a_i^T x - s_k = 0, lower[k] <= s_k <= upper[k] +// - rows: indices of the constraint rows of A that receive a slack. +// - is_range: selects the form above; the slack's entry in A is +1 or -1 respectively. +// - lower/upper: bounds of each new slack, indexed like rows. +// - new_slacks: receives the column index of each new slack. +// The slacks are inserted after the linear variables and before the cone variables, giving the +// layout [linear | slacks | cone]; slacks have zero objective cost and empty rows/columns in Q. +template +static void insert_slack_columns(lp_problem_t& problem, + const std::vector& rows, + bool is_range, + const std::vector& lower, + const std::vector& upper, + std::vector& new_slacks) +{ + const i_t num_new = static_cast(rows.size()); + if (num_new == 0) { return; } + const i_t old_num_cols = problem.num_cols; + const i_t insert_at = linear_variable_count(problem); + const i_t num_cols = old_num_cols + num_new; + + csc_matrix_t& A = problem.A; + const i_t nz_before = A.col_start[insert_at]; + // Insert the pre-slack columns into the matrix A + A.i.insert(A.i.begin() + nz_before, rows.begin(), rows.end()); + A.x.insert(A.x.begin() + nz_before, num_new, is_range ? f_t(-1) : f_t(1)); + A.col_start.insert(A.col_start.begin() + insert_at, num_new, 0); + for (i_t k = 0; k < num_new; ++k) { + A.col_start[insert_at + k] = nz_before + k; + new_slacks.push_back(insert_at + k); + } + for (i_t j = insert_at + num_new; j <= num_cols; ++j) { + A.col_start[j] += num_new; + } + A.n = num_cols; + A.nz_max += num_new; + // Insert the slack columns into the objective, lower, and upper bound vectors + problem.objective.insert(problem.objective.begin() + insert_at, num_new, f_t(0)); + problem.lower.insert(problem.lower.begin() + insert_at, lower.begin(), lower.end()); + problem.upper.insert(problem.upper.begin() + insert_at, upper.begin(), upper.end()); + + if (problem.Q.n > 0) { + // Slacks have empty Q rows and columns: add the empty rows, then shift the column indices + // that point past them. + csr_matrix_t& Q = problem.Q; + const i_t q_start = Q.row_start[insert_at]; + Q.row_start.insert(Q.row_start.begin() + insert_at + 1, num_new, q_start); + for (i_t& col : Q.j) { + if (col >= insert_at) { col += num_new; } + } + Q.m = num_cols; + Q.n = num_cols; + } + + problem.num_cols = num_cols; + // Update the cone variable start index for second-order cone constraints + if (!problem.second_order_cone_dims.empty()) { problem.cone_var_start += num_new; } +} + template i_t convert_less_than_to_equal(const user_problem_t& user_problem, std::vector& row_sense, @@ -864,144 +925,21 @@ i_t convert_less_than_to_equal(const user_problem_t& user_problem, } // We must convert rows in the form: a_i^T x <= beta // into: a_i^T x + s_i = beta, s_i >= 0 - - if (!problem.second_order_cone_dims.empty()) { - const i_t old_num_cols = problem.num_cols; - const i_t linear_cols = linear_variable_count(problem); - const i_t num_slacks = less_rows; - const i_t num_cols = old_num_cols + num_slacks; - const i_t old_nnz = problem.A.col_start[old_num_cols]; - const i_t nnz = old_nnz + num_slacks; - const i_t new_cone_start = linear_cols + num_slacks; - - auto old_A = problem.A; - csc_matrix_t expanded_A(problem.A.m, num_cols, nnz); - - std::vector objective(num_cols, 0.0); - std::vector lower(num_cols, 0.0); - std::vector upper(num_cols, INFINITY); - std::vector old_to_new(old_num_cols, -1); - - for (i_t j = 0; j < linear_cols; ++j) { - old_to_new[j] = j; - objective[j] = problem.objective[j]; - lower[j] = problem.lower[j]; - upper[j] = problem.upper[j]; - } - for (i_t j = linear_cols; j < old_num_cols; ++j) { - old_to_new[j] = j + num_slacks; - objective[old_to_new[j]] = problem.objective[j]; - lower[old_to_new[j]] = problem.lower[j]; - upper[old_to_new[j]] = problem.upper[j]; - } - - i_t nz = 0; - for (i_t j = 0; j < linear_cols; ++j) { - expanded_A.col_start[j] = nz; - for (i_t p = old_A.col_start[j]; p < old_A.col_start[j + 1]; ++p) { - expanded_A.i[nz] = old_A.i[p]; - expanded_A.x[nz] = old_A.x[p]; - ++nz; - } - } - - i_t slack_col = linear_cols; - for (i_t i = 0; i < problem.num_rows; i++) { - if (row_sense[i] == 'L') { - expanded_A.col_start[slack_col] = nz; - expanded_A.i[nz] = i; - expanded_A.x[nz] = 1.0; - new_slacks.push_back(slack_col); - row_sense[i] = 'E'; - ++slack_col; - ++nz; - --less_rows; - } - } - - for (i_t j = linear_cols; j < old_num_cols; ++j) { - i_t new_j = old_to_new[j]; - expanded_A.col_start[new_j] = nz; - for (i_t p = old_A.col_start[j]; p < old_A.col_start[j + 1]; ++p) { - expanded_A.i[nz] = old_A.i[p]; - expanded_A.x[nz] = old_A.x[p]; - ++nz; - } - } - expanded_A.col_start[num_cols] = nz; - assert(less_rows == 0); - assert(slack_col == new_cone_start); - assert(nz == nnz); - - if (problem.Q.n > 0) { - const auto old_Q = problem.Q; - const i_t q_nnz = old_Q.row_start[old_num_cols]; - - problem.Q.row_start.assign(num_cols + 1, 0); - for (i_t row = 0; row < old_num_cols; ++row) { - i_t new_row = old_to_new[row]; - problem.Q.row_start[new_row + 1] = old_Q.row_start[row + 1] - old_Q.row_start[row]; - } - for (i_t row = 0; row < num_cols; ++row) { - problem.Q.row_start[row + 1] += problem.Q.row_start[row]; - } - - problem.Q.j.resize(q_nnz); - problem.Q.x.resize(q_nnz); - auto row_starts = problem.Q.row_start; - for (i_t row = 0; row < old_num_cols; ++row) { - i_t new_row = old_to_new[row]; - for (i_t p = old_Q.row_start[row]; p < old_Q.row_start[row + 1]; ++p) { - problem.Q.j[row_starts[new_row]] = old_to_new[old_Q.j[p]]; - problem.Q.x[row_starts[new_row]] = old_Q.x[p]; - ++row_starts[new_row]; - } - } - problem.Q.m = num_cols; - problem.Q.n = num_cols; - problem.Q.nz_max = q_nnz; - } - - problem.A = expanded_A; - problem.A.n = num_cols; - problem.objective = objective; - problem.lower = lower; - problem.upper = upper; - problem.num_cols = num_cols; - problem.cone_var_start = new_cone_start; - return 0; - } - - i_t num_cols = problem.num_cols + less_rows; - i_t nnz = problem.A.col_start[problem.num_cols] + less_rows; - problem.A.col_start.resize(num_cols + 1); - problem.A.i.resize(nnz); - problem.A.x.resize(nnz); - problem.lower.resize(num_cols); - problem.upper.resize(num_cols); - problem.objective.resize(num_cols); - - i_t p = problem.A.col_start[problem.num_cols]; - i_t j = problem.num_cols; + std::vector rows; + rows.reserve(less_rows); for (i_t i = 0; i < problem.num_rows; i++) { if (row_sense[i] == 'L') { - problem.lower[j] = 0.0; - problem.upper[j] = INFINITY; - problem.objective[j] = 0.0; - problem.A.i[p] = i; - problem.A.x[p] = 1.0; - new_slacks.push_back(j); - problem.A.col_start[j++] = p++; - row_sense[i] = 'E'; - less_rows--; + rows.push_back(i); + row_sense[i] = 'E'; } } - problem.A.col_start[num_cols] = p; - assert(less_rows == 0); - assert(p == nnz); - problem.A.n = num_cols; - problem.num_cols = num_cols; + insert_slack_columns(problem, + rows, + false, + std::vector(rows.size(), 0.0), + std::vector(rows.size(), INFINITY), + new_slacks); return 0; } @@ -1074,23 +1012,17 @@ i_t convert_range_rows(const user_problem_t& user_problem, // // The values of h_i and u_i are determined by the b_i (RHS) and r_i (RANGES) // associated with the ith constraint as well as the row sense - i_t num_cols = problem.num_cols + user_problem.num_range_rows; - i_t num_range_rows = user_problem.num_range_rows; - i_t nnz = problem.A.col_start[problem.num_cols] + num_range_rows; - problem.A.col_start.resize(num_cols + 1); - problem.A.i.resize(nnz); - problem.A.x.resize(nnz); - problem.lower.resize(num_cols); - problem.upper.resize(num_cols); - problem.objective.resize(num_cols); - - i_t p = problem.A.col_start[problem.num_cols]; - i_t j = problem.num_cols; + const i_t num_range_rows = user_problem.num_range_rows; + std::vector rows(num_range_rows); + std::vector lower(num_range_rows); + std::vector upper(num_range_rows); for (i_t k = 0; k < num_range_rows; k++) { const i_t i = user_problem.range_rows[k]; const f_t r = user_problem.range_value[k]; const f_t b = problem.rhs[i]; - auto [lower, upper] = get_range_bounds_from_sense(row_sense[i], b, r); + const auto [lo, hi] = get_range_bounds_from_sense(row_sense[i], b, r); + lower[k] = lo; + upper[k] = hi; if (row_sense[i] == 'L') { less_rows--; equal_rows++; @@ -1098,21 +1030,11 @@ i_t convert_range_rows(const user_problem_t& user_problem, greater_rows--; equal_rows++; } - problem.lower[j] = lower; - problem.upper[j] = upper; - problem.objective[j] = 0.0; - problem.A.i[p] = i; - problem.A.x[p] = -1.0; - new_slacks.push_back(j); - problem.A.col_start[j++] = p++; - problem.rhs[i] = 0.0; - row_sense[i] = 'E'; + rows[k] = i; + problem.rhs[i] = 0.0; + row_sense[i] = 'E'; } - problem.A.col_start[num_cols] = p; - assert(p == nnz); - problem.A.n = num_cols; - problem.num_cols = num_cols; - + insert_slack_columns(problem, rows, true, lower, upper, new_slacks); return 0; } @@ -1327,6 +1249,18 @@ void convert_user_problem(const user_problem_t& user_problem, problem.cone_var_start = user_problem.cone_var_start; problem.second_order_cone_dims = user_problem.second_order_cone_dims; + if (user_problem.Q_values.size() > 0) { + settings.log.debug("Converting problem with %d quadratic nonzeros\n", + user_problem.Q_values.size()); + problem.Q.m = user_problem.num_cols; + problem.Q.n = user_problem.num_cols; + problem.Q.nz_max = user_problem.Q_values.size(); + problem.Q.row_start.assign(user_problem.Q_offsets.begin(), + user_problem.Q_offsets.begin() + user_problem.num_cols + 1); + problem.Q.j = user_problem.Q_indices; + problem.Q.x = user_problem.Q_values; + } + // Make a copy of row_sense so we can modify it std::vector row_sense = user_problem.row_sense; @@ -1534,26 +1468,6 @@ void convert_user_problem(const user_problem_t& user_problem, convert_less_than_to_equal(user_problem, row_sense, problem, less_rows, new_slacks); } - if (user_problem.Q_values.size() > 0) { - settings.log.debug("Converting problem with %d quadratic nonzeros\n", - user_problem.Q_values.size()); - settings.log.debug( - "problem.num_cols: %d user_problem.num_cols: %d\n", problem.num_cols, user_problem.num_cols); - problem.Q.m = problem.num_cols; - problem.Q.n = problem.num_cols; - problem.Q.nz_max = user_problem.Q_values.size(); - problem.Q.row_start.resize(problem.num_cols + 1); - for (i_t j = 0; j < user_problem.num_cols; j++) { - problem.Q.row_start[j] = user_problem.Q_offsets[j]; - } - i_t nz = user_problem.Q_offsets[user_problem.num_cols]; - for (i_t j = user_problem.num_cols; j <= problem.num_cols; j++) { - problem.Q.row_start[j] = nz; - } - problem.Q.j = user_problem.Q_indices; - problem.Q.x = user_problem.Q_values; - } - // Add artifical variables if (!settings.barrier_presolve) { add_artifical_variables(problem, user_problem.range_rows, equality_rows, new_slacks); diff --git a/cpp/tests/socp/solve_barrier_socp.cu b/cpp/tests/socp/solve_barrier_socp.cu index 550a7e3021..9a9bbf908e 100644 --- a/cpp/tests/socp/solve_barrier_socp.cu +++ b/cpp/tests/socp/solve_barrier_socp.cu @@ -159,6 +159,145 @@ TEST(barrier, cone_metadata_reindexed_when_slack_is_inserted_before_cones) EXPECT_EQ(barrier_lp.cone_var_start, 2); } +TEST(barrier, range_and_less_rows_with_soc_block) +{ + // Variables ordered as [a, b | t, u, v], where (t, u, v) \in Q^3. + // + // minimize -a + t + // subject to a + b <= 3 (<= row) + // b - u = 0 + // 1 <= b <= 2 (range row) + // (t, u, v) in Q^3 + // + // The <= row and the range row each add a slack column, and both must land in front of the + // cone block. With v = 0 and t = u = b the objective is -a + b, and a <= 3 - b gives + // -3 + 2b, minimized at the range row's lower end b = 1. + // + // Optimal: a* = 2, b* = 1, t* = 1, u* = 1, v* = 0, obj* = -1. + raft::handle_t handle{}; + init_handler(&handle); + + user_problem_t user_problem(&handle); + + constexpr int m = 3; + constexpr int n = 5; + constexpr int nz = 5; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {-1.0, 0.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); + // Columns: a, b, t, u, v + user_problem.A.col_start = {0, 1, 4, 4, 5, 5}; + user_problem.A.i = {0, 0, 1, 2, 1}; + user_problem.A.x = {1.0, 1.0, 1.0, 1.0, -1.0}; + + // Range rows arrive as 'E' rows with rhs = lower bound and range = upper - lower. + user_problem.rhs = {3.0, 0.0, 1.0}; + user_problem.row_sense = {'L', 'E', 'E'}; + user_problem.lower.assign(n, 0.0); + user_problem.upper.assign(n, inf); + + user_problem.num_range_rows = 1; + user_problem.range_rows = {2}; + user_problem.range_value = {1.0}; + user_problem.problem_name = "range_and_less_rows_with_soc_block"; + user_problem.cone_var_start = 2; + user_problem.second_order_cone_dims = {3}; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + + lp_solution_t solution(m, n); + auto status = solve_linear_program_with_barrier(user_problem, settings, solution); + + EXPECT_EQ(status, lp_status_t::OPTIMAL); + EXPECT_NEAR(solution.objective, -1.0, 1e-4); + EXPECT_NEAR(solution.x[0], 2.0, 1e-4); + EXPECT_NEAR(solution.x[1], 1.0, 1e-4); + EXPECT_NEAR(solution.x[2], 1.0, 1e-4); + EXPECT_NEAR(solution.x[3], 1.0, 1e-4); + EXPECT_NEAR(std::abs(solution.x[4]), 0.0, 1e-4); +} + +TEST(barrier, qp_with_range_and_less_rows_and_soc_block) +{ + // Variables ordered as [l | t, u, v], where (t, u, v) \in Q^3. + // + // minimize 0.5 l^2 + 0.5 t^2 + // subject to l - u <= -1 (<= row) + // 2 <= l + u <= 4 (range row) + // (t, u, v) in Q^3 + // + // Q has an entry on the cone head t, so it has to move along with the cone block when the two + // slack columns are inserted in front of it. Since t >= |u|, v = 0 and t = u at the optimum, + // leaving 0.5 l^2 + 0.5 u^2. The unconstrained point l = u = 1 violates l - u <= -1, so the + // optimum sits where u - l = 1 meets the range row's lower end l + u = 2 (multipliers 1 and + // 0.5, both strictly positive). + // + // Optimal: l* = 0.5, t* = 1.5, u* = 1.5, v* = 0, obj* = 1.25. + raft::handle_t handle{}; + init_handler(&handle); + + user_problem_t user_problem(&handle); + + constexpr int m = 2; + constexpr int n = 4; + constexpr int nz = 4; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective.assign(n, 0.0); + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + // Columns: l, t, u, v + user_problem.A.col_start = {0, 2, 2, 4, 4}; + user_problem.A.i = {0, 1, 0, 1}; + user_problem.A.x = {1.0, 1.0, -1.0, 1.0}; + + user_problem.rhs = {-1.0, 2.0}; + user_problem.row_sense = {'L', 'E'}; + user_problem.lower.assign(n, 0.0); + user_problem.upper.assign(n, inf); + + user_problem.Q_offsets = {0, 1, 2, 2, 2}; + user_problem.Q_indices = {0, 1}; + user_problem.Q_values = {1.0, 1.0}; + + user_problem.num_range_rows = 1; + user_problem.range_rows = {1}; + user_problem.range_value = {2.0}; + user_problem.problem_name = "qp_with_range_and_less_rows_and_soc_block"; + user_problem.cone_var_start = 1; + user_problem.second_order_cone_dims = {3}; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + + lp_solution_t solution(m, n); + auto status = solve_linear_program_with_barrier(user_problem, settings, solution); + + EXPECT_EQ(status, lp_status_t::OPTIMAL); + EXPECT_NEAR(solution.objective, 1.25, 1e-4); + EXPECT_NEAR(solution.x[0], 0.5, 1e-4); + EXPECT_NEAR(solution.x[1], 1.5, 1e-4); + EXPECT_NEAR(solution.x[2], 1.5, 1e-4); + EXPECT_NEAR(std::abs(solution.x[3]), 0.0, 1e-4); +} + TEST(barrier, presolve_reindexes_cone_start_after_empty_column_removal) { raft::handle_t handle{};