diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index f17e53ab98..93de531290 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -2203,6 +2203,8 @@ class iteration_data_t { f_t dual_residual_norm_save; f_t complementarity_residual_norm_save; + bool use_high_accuracy_ir = false; + dense_vector_t diag; pinned_dense_vector_t inv_diag; dense_vector_t inv_sqrt_diag; @@ -2581,7 +2583,8 @@ int barrier_solver_t::initial_point(iteration_data_t& data) } op(data); if (settings.barrier_iterative_refinement) { - const f_t ir_tol = data.has_sparse_cones() ? f_t(1e-12) : f_t(1e-8); + const f_t ir_tol = + (data.has_sparse_cones() || data.use_high_accuracy_ir) ? f_t(1e-12) : f_t(1e-8); iterative_refinement(op, rhs, soln, ir_tol); } @@ -3175,7 +3178,8 @@ i_t barrier_solver_t::gpu_compute_search_direction(iteration_data_t( op, data.d_augmented_rhs_, data.d_augmented_soln_, ir_tol); if (solve_err > 1e-1) { @@ -3263,8 +3267,9 @@ i_t barrier_solver_t::gpu_compute_search_direction(iteration_data_t(adat_op, data.d_h_, data.d_dy_); + iterative_refinement(adat_op, data.d_h_, data.d_dy_, ir_tol); if (adat_solve_err > 1e-1) { settings.log.debug("||ADAT*dy - h|| %e after IR\n", adat_solve_err); } @@ -4266,16 +4271,20 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( f_t& primal_residual_norm, f_t& dual_residual_norm, f_t& complementarity_residual_norm, + f_t& objective_gap, f_t& relative_primal_residual, f_t& relative_dual_residual, f_t& relative_complementarity_residual, + f_t& relative_objective_gap, lp_solution_t& solution) { raft::common::nvtx::range fun_scope("Barrier: check_for_suboptimal_solution"); + bool small_gap = (!data.has_cones() && data.Q.n == 0) || + relative_objective_gap < settings.barrier_relaxed_relative_objective_gap_tol; if (relative_primal_residual < settings.barrier_relaxed_feasibility_tol && relative_dual_residual < settings.barrier_relaxed_optimality_tol && relative_complementarity_residual < settings.barrier_relaxed_complementarity_tol && - primal_objective == primal_objective) { + small_gap && primal_objective == primal_objective) { raft::copy(data.x.data(), data.d_x_.data(), data.d_x_.size(), stream_view_); raft::copy(data.y.data(), data.d_y_.data(), data.d_y_.size(), stream_view_); raft::copy(data.z.data(), data.d_z_.data(), data.d_z_.size(), stream_view_); @@ -4300,12 +4309,16 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( settings.log.printf("Complementarity gap (abs/rel): %8.2e/%8.2e\n", complementarity_residual_norm, relative_complementarity_residual); + settings.log.printf( + "Objective gap (abs/rel): %8.2e/%8.2e\n", objective_gap, relative_objective_gap); settings.log.printf("\n"); return lp_status_t::OPTIMAL; // TODO: Barrier should probably have a separate suboptimal // status } f_t primal_objective_save = data.c.inner_product(data.x_save); + f_t dual_objective_save = + data.b.inner_product(data.y_save) - data.restrict_u_.inner_product(data.v_save); if (data.Q.n > 0) { dense_vector_t Qx_save(data.Q.n); dense_vector_t x_save_host(data.Q.n); @@ -4313,11 +4326,22 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( matrix_vector_multiply(data.Q, 1.0, x_save_host, 0.0, Qx_save); f_t quad_objective = 0.5 * x_save_host.inner_product(Qx_save); primal_objective_save += quad_objective; + dual_objective_save -= quad_objective; } + f_t objective_gap_save = std::abs(primal_objective_save - dual_objective_save); + f_t user_primal_objective_save = compute_user_objective(lp, primal_objective_save); + f_t relative_objective_gap_save = + objective_gap_save / + (1.0 + std::min(std::abs(user_primal_objective_save), std::abs(primal_objective_save))); + bool small_gap_save = + (!data.has_cones() && data.Q.n == 0) || + relative_objective_gap_save < settings.barrier_relaxed_relative_objective_gap_tol; + if (data.relative_primal_residual_save < settings.barrier_relaxed_feasibility_tol && data.relative_dual_residual_save < settings.barrier_relaxed_optimality_tol && - data.relative_complementarity_residual_save < settings.barrier_relaxed_complementarity_tol) { + data.relative_complementarity_residual_save < settings.barrier_relaxed_complementarity_tol && + small_gap_save) { settings.log.printf("Restoring previous solution\n"); data.restore_saved_iterate(); data.to_solution(lp, @@ -4340,14 +4364,19 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( settings.log.printf("Complementarity gap (abs/rel): %8.2e/%8.2e\n", data.complementarity_residual_norm_save, data.relative_complementarity_residual_save); + settings.log.printf("Objective gap (abs/rel): %8.2e/%8.2e\n", + objective_gap_save, + relative_objective_gap_save); settings.log.printf("\n"); return lp_status_t::OPTIMAL; // TODO: Barrier should probably have a separate suboptimal // status } else { - settings.log.printf("Primal residual %.2e dual residual %.2e complementarity residual %.2e\n", - relative_primal_residual, - relative_dual_residual, - relative_complementarity_residual); + settings.log.printf( + "Primal residual %.2e dual residual %.2e complementarity residual %.2e objective gap %.2e\n", + relative_primal_residual, + relative_dual_residual, + relative_complementarity_residual, + relative_objective_gap); } settings.log.printf("Search direction computation failed\n"); return lp_status_t::NUMERICAL_ISSUES; @@ -4449,18 +4478,18 @@ lp_status_t barrier_solver_t::barrier_advanced_solve(f_t start_time, mu, primal_objective, dual_objective); + f_t user_primal_objective = compute_user_objective(lp, primal_objective); + f_t user_dual_objective = compute_user_objective(lp, dual_objective); f_t relative_primal_residual = primal_residual_norm / (1.0 + norm_b); f_t relative_dual_residual = dual_residual_norm / (1.0 + norm_c); f_t relative_complementarity_residual = complementarity_residual_norm / - (1.0 + std::min(std::abs(compute_user_objective(lp, primal_objective)), - std::abs(primal_objective))); + (1.0 + std::min(std::abs(user_primal_objective), std::abs(primal_objective))); - f_t objective_gap_abs = std::abs(primal_objective - dual_objective); - f_t objective_gap_rel = - objective_gap_abs / - std::max(f_t(1), std::min(std::abs(primal_objective), std::abs(dual_objective))); + f_t objective_gap, relative_objective_gap; + compute_objective_gap( + lp, primal_objective, dual_objective, objective_gap, relative_objective_gap); data.w_save = data.w; data.x_save = data.x; @@ -4477,16 +4506,19 @@ lp_status_t barrier_solver_t::barrier_advanced_solve(f_t start_time, float64_t elapsed_time = toc(start_time); settings.log.printf("%3d %+.12e %+.12e %.2e %.2e %.2e %.1f\n", iter, - compute_user_objective(lp, primal_objective), - compute_user_objective(lp, dual_objective), + user_primal_objective, + user_dual_objective, relative_primal_residual, relative_dual_residual, relative_complementarity_residual, elapsed_time); - bool converged = primal_residual_norm < settings.barrier_relative_feasibility_tol && - dual_residual_norm < settings.barrier_relative_optimality_tol && - complementarity_residual_norm < settings.barrier_relative_complementarity_tol; + bool small_gap = (!data.has_cones() && data.Q.n == 0) || + relative_objective_gap < settings.barrier_relaxed_relative_objective_gap_tol; + bool converged = + primal_residual_norm < settings.barrier_relative_feasibility_tol && + dual_residual_norm < settings.barrier_relative_optimality_tol && + complementarity_residual_norm < settings.barrier_relative_complementarity_tol && small_gap; const i_t iteration_limit = settings.iteration_limit; @@ -4535,9 +4567,11 @@ lp_status_t barrier_solver_t::barrier_advanced_solve(f_t start_time, primal_residual_norm, dual_residual_norm, complementarity_residual_norm, + objective_gap, relative_primal_residual, relative_dual_residual, relative_complementarity_residual, + relative_objective_gap, solution); } if (toc(start_time) > settings.time_limit) { @@ -4575,9 +4609,11 @@ lp_status_t barrier_solver_t::barrier_advanced_solve(f_t start_time, primal_residual_norm, dual_residual_norm, complementarity_residual_norm, + objective_gap, relative_primal_residual, relative_dual_residual, relative_complementarity_residual, + relative_objective_gap, solution); } data.has_factorization = false; @@ -4605,17 +4641,20 @@ lp_status_t barrier_solver_t::barrier_advanced_solve(f_t start_time, primal_objective, dual_objective); - relative_primal_residual = primal_residual_norm / (1.0 + norm_b); - relative_dual_residual = dual_residual_norm / (1.0 + norm_c); + f_t user_primal_objective = compute_user_objective(lp, primal_objective); + relative_primal_residual = primal_residual_norm / (1.0 + norm_b); + relative_dual_residual = dual_residual_norm / (1.0 + norm_c); relative_complementarity_residual = complementarity_residual_norm / - (1.0 + std::min(std::abs(compute_user_objective(lp, primal_objective)), - std::abs(primal_objective))); + (1.0 + std::min(std::abs(user_primal_objective), std::abs(primal_objective))); - objective_gap_abs = std::abs(primal_objective - dual_objective); - objective_gap_rel = - objective_gap_abs / - std::max(f_t(1), std::min(std::abs(primal_objective), std::abs(dual_objective))); + compute_objective_gap( + lp, primal_objective, dual_objective, objective_gap, relative_objective_gap); + + if (data.has_cones() && !data.use_high_accuracy_ir && relative_primal_residual < 1e-6 && + relative_dual_residual < 1e-6 && relative_complementarity_residual < 1e-6) { + data.use_high_accuracy_ir = true; + } if (relative_primal_residual < settings.barrier_relaxed_feasibility_tol && relative_dual_residual < settings.barrier_relaxed_optimality_tol && @@ -4663,9 +4702,11 @@ lp_status_t barrier_solver_t::barrier_advanced_solve(f_t start_time, primal_residual_norm, dual_residual_norm, complementarity_residual_norm, + objective_gap, relative_primal_residual, relative_dual_residual, relative_complementarity_residual, + relative_objective_gap, solution); } @@ -4683,7 +4724,8 @@ lp_status_t barrier_solver_t::barrier_advanced_solve(f_t start_time, bool small_gap = relative_complementarity_residual < settings.barrier_relative_complementarity_tol; bool small_objective_gap = - !data.has_cones() || objective_gap_rel < settings.barrier_relaxed_complementarity_tol; + (!data.has_cones() && data.Q.n == 0) || + relative_objective_gap < settings.barrier_relative_objective_gap_tol; converged = primal_feasible && dual_feasible && small_gap && small_objective_gap; @@ -4701,6 +4743,8 @@ lp_status_t barrier_solver_t::barrier_advanced_solve(f_t start_time, settings.log.printf("Complementarity gap (abs/rel): %8.2e/%8.2e\n", complementarity_residual_norm, relative_complementarity_residual); + settings.log.printf( + "Objective gap (abs/rel): %8.2e/%8.2e\n", objective_gap, relative_objective_gap); settings.log.printf("\n"); raft::copy(data.x.data(), data.d_x_.data(), data.d_x_.size(), stream_view_); raft::copy(data.y.data(), data.d_y_.data(), data.d_y_.size(), stream_view_); @@ -4735,9 +4779,11 @@ lp_status_t barrier_solver_t::barrier_advanced_solve(f_t start_time, primal_residual_norm, dual_residual_norm, complementarity_residual_norm, + objective_gap, relative_primal_residual, relative_dual_residual, relative_complementarity_residual, + relative_objective_gap, solution); } } diff --git a/cpp/src/barrier/barrier.hpp b/cpp/src/barrier/barrier.hpp index 58fb6d63f5..b6ce252012 100644 --- a/cpp/src/barrier/barrier.hpp +++ b/cpp/src/barrier/barrier.hpp @@ -108,9 +108,11 @@ class barrier_solver_t { f_t& primal_residual_norm, f_t& dual_residual_norm, f_t& complementarity_residual_norm, + f_t& objective_gap, f_t& relative_primal_residual, f_t& relative_dual_residual, f_t& relative_complementarity_residual, + f_t& relative_objective_gap, simplex::lp_solution_t& solution); const simplex::lp_problem_t& lp; diff --git a/cpp/src/barrier/iterative_refinement.hpp b/cpp/src/barrier/iterative_refinement.hpp index 9391dd0b6d..7680dfe610 100644 --- a/cpp/src/barrier/iterative_refinement.hpp +++ b/cpp/src/barrier/iterative_refinement.hpp @@ -58,7 +58,7 @@ template f_t iterative_refinement_simple(T& op, const rmm::device_uvector& b, rmm::device_uvector& x, - f_t tol = 1e-8) + f_t tol) { rmm::device_uvector x_sav(x, x.stream()); @@ -125,7 +125,7 @@ template f_t iterative_refinement_gmres(T& op, const rmm::device_uvector& b, rmm::device_uvector& x, - f_t tol = 1e-8) + f_t tol) { // Parameters // Ideally, we do not need to restart here. But having restarts helps as a checkpoint to get @@ -383,7 +383,7 @@ template f_t iterative_refinement(T& op, const rmm::device_uvector& b, rmm::device_uvector& x, - f_t tol = 1e-8) + f_t tol) { return iterative_refinement_gmres(op, b, x, tol); } diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 647b11f48e..33d8f56371 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -49,9 +49,11 @@ struct simplex_solver_settings_t { barrier_relative_feasibility_tol(1e-8), barrier_relative_optimality_tol(1e-8), barrier_relative_complementarity_tol(1e-8), + barrier_relative_objective_gap_tol(1e-6), barrier_relaxed_feasibility_tol(1e-4), barrier_relaxed_optimality_tol(1e-4), barrier_relaxed_complementarity_tol(1e-4), + barrier_relaxed_relative_objective_gap_tol(1e-4), cut_off(std::numeric_limits::infinity()), steepest_edge_ratio(0.5), steepest_edge_primal_tol(1e-9), @@ -149,9 +151,12 @@ struct simplex_solver_settings_t { f_t barrier_relative_optimality_tol; // Relative optimality tolerance for barrier method f_t barrier_relative_complementarity_tol; // Relative complementarity tolerance for barrier method + f_t barrier_relative_objective_gap_tol; // Relative objective gap tolerance for barrier method f_t barrier_relaxed_feasibility_tol; // Relative feasibility tolerance for barrier method f_t barrier_relaxed_optimality_tol; // Relative optimality tolerance for barrier method f_t barrier_relaxed_complementarity_tol; // Relative complementarity tolerance for barrier method + f_t barrier_relaxed_relative_objective_gap_tol; // Relative objective gap tolerance for barrier + // method f_t cut_off; // If the dual objective is greater than the cutoff we stop f_t steepest_edge_ratio; // the ratio of computed steepest edge mismatch from updated steepest edge diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 0fcb973ea5..12f042a931 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -221,6 +221,24 @@ f_t compute_user_objective(const lp_problem_t& lp, f_t obj) return user_obj; } +template +void compute_objective_gap(const lp_problem_t& problem, + f_t primal_obj, + f_t dual_obj, + f_t& objective_gap, + f_t& relative_objective_gap) +{ + objective_gap = std::abs(primal_obj - dual_obj); + f_t user_primal_obj = compute_user_objective(problem, primal_obj); + f_t user_dual_obj = compute_user_objective(problem, dual_obj); + + f_t denom_1 = std::min(std::abs(user_primal_obj), std::abs(primal_obj)); + f_t denom_2 = std::min(std::abs(user_dual_obj), std::abs(dual_obj)); + f_t denom = 1.0 + std::max(denom_1, denom_2); + + relative_objective_gap = objective_gap / denom; +} + template f_t compute_presolved_objective(const lp_problem_t& lp, f_t user_obj) { @@ -1133,6 +1151,12 @@ template double compute_user_objective(const lp_problem_t& lp, double obj); +template void compute_objective_gap(const lp_problem_t& problem, + double primal_obj, + double dual_obj, + double& objective_gap, + double& relative_objective_gap); + template double compute_presolved_objective(const lp_problem_t& lp, double user_obj); template lp_status_t solve_linear_program_advanced( diff --git a/cpp/src/dual_simplex/solve.hpp b/cpp/src/dual_simplex/solve.hpp index a45913c144..08ae11a9dc 100644 --- a/cpp/src/dual_simplex/solve.hpp +++ b/cpp/src/dual_simplex/solve.hpp @@ -67,6 +67,13 @@ f_t compute_user_objective(const lp_problem_t& lp, const std::vector f_t compute_user_objective(const lp_problem_t& lp, f_t obj); +template +void compute_objective_gap(const lp_problem_t& problem, + f_t primal_obj, + f_t dual_obj, + f_t& objective_gap, + f_t& relative_objective_gap); + template f_t compute_presolved_objective(const lp_problem_t& lp, f_t user_obj); diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index a6d7b61011..8c40913f06 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -562,6 +562,8 @@ std::tuple, simplex::lp_status_t, f_t, f_t, f_t barrier_settings.barrier_relaxed_feasibility_tol = settings.tolerances.relative_primal_tolerance; barrier_settings.barrier_relaxed_optimality_tol = settings.tolerances.relative_dual_tolerance; barrier_settings.barrier_relaxed_complementarity_tol = settings.tolerances.relative_gap_tolerance; + barrier_settings.barrier_relaxed_relative_objective_gap_tol = + settings.tolerances.relative_gap_tolerance; if (barrier_settings.concurrent_halt != nullptr) { // Don't show the barrier log in concurrent mode. Show the PDLP log instead barrier_settings.log.log = false;