Skip to content

Include objective gap in QP termination criteria #1733

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rg20 wants to merge 9 commits into
base: main
Choose a base branch
from
Open

Include objective gap in QP termination criteria #1733

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
e1b24f2
Add objective gap termination criteria for QP
rg20 Aug 14, 2026
269cbd9
Add objective gap reporting and checks in suboptimal cases
rg20 Aug 15, 2026
29bce4f
Rename the relaxed objective tolerance
rg20 Aug 17, 2026
ed622ed
Adjust the tolerance to 1e-8
rg20 Aug 17, 2026
1fedd5a
Revert the tolerance
rg20 Aug 18, 2026
3936874
Merge remote-tracking branch 'upstream/main' into qp_objective_checks
rg20 Aug 20, 2026
d9513a7
Merge remote-tracking branch 'upstream/main' into qp_objective_checks
rg20 Aug 24, 2026
de0bd77
Merge remote-tracking branch 'upstream/main' into qp_objective_checks
rg20 Aug 26, 2026
6742aae
Address PR comments
rg20 Aug 26, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
90 changes: 63 additions & 27 deletions cpp/src/barrier/barrier.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4015,16 +4015,20 @@ lp_status_t barrier_solver_t<i_t, f_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<i_t, f_t>& solution)
{
raft::common::nvtx::range fun_scope("Barrier: check_for_suboptimal_solution");
bool small_gap = (!data.has_cones() && data.Q.n == 0) ||

@yuwenchen95 yuwenchen95 Aug 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do we need (!data.has_cones() && data.Q.n == 0)? I think we should have duality gap check for all kinds of problems solved by barrier.

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 &&

@chris-maes chris-maes Aug 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The purpose of primal_objective == primal_objective is not record solutions that lead to NaN in the objective. Maybe this is no longer necessary with small_gap.

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_);
Expand All @@ -4049,24 +4053,39 @@ lp_status_t barrier_solver_t<i_t, f_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<i_t, f_t> Qx_save(data.Q.n);
dense_vector_t<i_t, f_t> x_save_host(data.Q.n);
std::copy(data.x_save.begin(), data.x_save.begin() + data.Q.n, x_save_host.begin());
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 =

@yuwenchen95 yuwenchen95 Aug 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Only primal info is used in computingrelative_objective_gap_save. We should also use dual info for the denominator.

@chris-maes chris-maes Aug 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How would you use that info in the denominator? Take the min over the primal and dual objectives?

@yuwenchen95 yuwenchen95 Aug 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we can take the max over the absolute values of primal and dual, and then min operation over solver objective and user objective.

@rg20 rg20 Aug 20, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why max though? min is much more conservative

@yuwenchen95 yuwenchen95 Aug 20, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, min is stricter. The point is to include both primal and dual info for computation at this line.

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,
Expand All @@ -4089,14 +4108,19 @@ lp_status_t barrier_solver_t<i_t, f_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;
Expand Down Expand Up @@ -4236,14 +4260,14 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
matrix_vector_multiply(data.Q, 1.0, data.x, 0.0, Qx);
quad_objective = 0.5 * data.x.inner_product(Qx);
}
f_t primal_objective = data.c.inner_product(data.x) + quad_objective;
f_t primal_objective = data.c.inner_product(data.x) + quad_objective;
f_t user_primal_objective = compute_user_objective(lp, primal_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)));

@yuwenchen95 yuwenchen95 Aug 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should also use dual info for the denominator here.


dense_vector_t<i_t, f_t> upper(lp.upper);
data.gather_upper_bounds(upper, data.restrict_u_);
Expand All @@ -4252,11 +4276,11 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
data.d_restrict_u_.data(), data.restrict_u_.data(), data.restrict_u_.size(), stream_view_);
f_t dual_objective =
data.b.inner_product(data.y) - data.restrict_u_.inner_product(data.v) - quad_objective;
f_t user_dual_objective = compute_user_objective(lp, dual_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;
Expand All @@ -4273,16 +4297,19 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
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) ||

@yuwenchen95 yuwenchen95 Aug 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

!data.has_cones() && data.Q.n == 0): we may want duality check for all problems.

@chris-maes chris-maes Aug 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's test if we can do this for all problems. Hopefully, we can.

@rg20 rg20 Aug 20, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not yet, we have regressions on LP barrier

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;

Expand Down Expand Up @@ -4327,9 +4354,11 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
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) {
Expand Down Expand Up @@ -4367,9 +4396,11 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
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;
Expand All @@ -4396,17 +4427,15 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,

compute_primal_dual_objective(data, 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)));

@yuwenchen95 yuwenchen95 Aug 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same issue for missing dual info.


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 (relative_primal_residual < settings.barrier_relaxed_feasibility_tol &&
relative_dual_residual < settings.barrier_relaxed_optimality_tol &&
Expand Down Expand Up @@ -4454,9 +4483,11 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
primal_residual_norm,
dual_residual_norm,
complementarity_residual_norm,
objective_gap,
relative_primal_residual,
relative_dual_residual,
relative_complementarity_residual,
relative_objective_gap,
solution);
}

Expand All @@ -4474,7 +4505,8 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
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) ||

@yuwenchen95 yuwenchen95 Aug 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same concern for (!data.has_cones() && data.Q.n == 0) above.

relative_objective_gap < settings.barrier_relative_objective_gap_tol;

converged = primal_feasible && dual_feasible && small_gap && small_objective_gap;

Expand All @@ -4492,6 +4524,8 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
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);
Comment thread
coderabbitai[bot] marked this conversation as resolved.
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_);
Expand Down Expand Up @@ -4526,9 +4560,11 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
primal_residual_norm,
dual_residual_norm,
complementarity_residual_norm,
objective_gap,
relative_primal_residual,
relative_dual_residual,
relative_complementarity_residual,
relative_objective_gap,
solution);
}
}
Expand Down
2 changes: 2 additions & 0 deletions cpp/src/barrier/barrier.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,9 +109,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<i_t, f_t>& solution);

const simplex::lp_problem_t<i_t, f_t>& lp;
Expand Down
5 changes: 5 additions & 0 deletions cpp/src/dual_simplex/simplex_solver_settings.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,9 +48,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<f_t>::infinity()),
steepest_edge_ratio(0.5),
steepest_edge_primal_tol(1e-9),
Expand Down Expand Up @@ -141,9 +143,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
Expand Down
24 changes: 24 additions & 0 deletions cpp/src/dual_simplex/solve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,24 @@ f_t compute_user_objective(const lp_problem_t<i_t, f_t>& lp, f_t obj)
return user_obj;
}

template <typename i_t, typename f_t>
void compute_objective_gap(const lp_problem_t<i_t, f_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 <typename i_t, typename f_t>
f_t compute_presolved_objective(const lp_problem_t<i_t, f_t>& lp, f_t user_obj)
{
Expand Down Expand Up @@ -819,6 +837,12 @@ template double compute_user_objective<int, double>(const lp_problem_t<int, doub

template double compute_user_objective(const lp_problem_t<int, double>& lp, double obj);

template void compute_objective_gap<int, double>(const lp_problem_t<int, double>& problem,
double primal_obj,
double dual_obj,
double& objective_gap,
double& relative_objective_gap);

template double compute_presolved_objective(const lp_problem_t<int, double>& lp, double user_obj);

template lp_status_t solve_linear_program_advanced(
Expand Down
7 changes: 7 additions & 0 deletions cpp/src/dual_simplex/solve.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,13 @@ f_t compute_user_objective(const lp_problem_t<i_t, f_t>& lp, const std::vector<f
template <typename i_t, typename f_t>
f_t compute_user_objective(const lp_problem_t<i_t, f_t>& lp, f_t obj);

template <typename i_t, typename f_t>
void compute_objective_gap(const lp_problem_t<i_t, f_t>& problem,
f_t primal_obj,
f_t dual_obj,
f_t& objective_gap,
f_t& relative_objective_gap);

template <typename i_t, typename f_t>
f_t compute_presolved_objective(const lp_problem_t<i_t, f_t>& lp, f_t user_obj);

Expand Down
2 changes: 2 additions & 0 deletions cpp/src/pdlp/solve.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -523,6 +523,8 @@ std::tuple<simplex::lp_solution_t<i_t, f_t>, 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;
Expand Down