Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
35 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
51b7dc8
Eliminate zero-cost free linear variables by equality substitution in…
rg20 Sep 16, 2026
cd9e70c
Do not sort row entries before the CSR to CSC conversion
rg20 Sep 16, 2026
02a258e
Merge upstream/main into eliminate_free_variables_in_qp_socp
rg20 Sep 16, 2026
168a1b4
Replace per-row hash maps with flat CSR/CSC arenas for free-variable …
rg20 Sep 16, 2026
711f1f8
Merge remote-tracking branch 'upstream/main' into eliminate_free_vari…
rg20 Sep 21, 2026
6b0f332
Use threshold pivoting when eliminating free variables by substitution.
rg20 Sep 22, 2026
d397bbd
Hardcode the free-variable elimination pivot thresholds.
rg20 Sep 22, 2026
89a6f8d
Merge remote-tracking branch 'upstream/main' into eliminate_free_vari…
rg20 Sep 22, 2026
79a2898
Fix an indexing bug
rg20 Sep 22, 2026
f8668f6
Use TMPDIR instead of explicit tmp
rg20 Sep 22, 2026
3cd2e17
Resize uncrush vectors appropriately
rg20 Sep 23, 2026
1ede03d
Merge remote-tracking branch 'upstream/main' into eliminate_free_vari…
rg20 Sep 23, 2026
8dd5d19
Fix CI failures
rg20 Sep 23, 2026
e3e1e0f
Disable CR scaling for PDLP test
rg20 Sep 23, 2026
47e106d
Address PR review
rg20 Sep 23, 2026
9625c54
disable CR scaling for warm start test
rg20 Sep 23, 2026
97065f7
Address remaining PR review comments.
rg20 Sep 23, 2026
04e8c46
Cache row maxima when choosing a substitution pivot.
rg20 Sep 23, 2026
4a72eb3
Merge remote-tracking branch 'upstream/main' into qp_objective_checks
rg20 Sep 23, 2026
3c17339
fix style checks
rg20 Sep 23, 2026
f899ef4
Merge remote-tracking branch 'upstream/main' into eliminate_free_vari…
rg20 Sep 23, 2026
0837cad
Merge branch 'eliminate_free_variables_in_qp_socp' into qp_objective_…
rg20 Sep 23, 2026
e3d42bb
Merge remote-tracking branch 'upstream/main' into qp_objective_checks
rg20 Sep 24, 2026
74be02a
Merge remote-tracking branch 'upstream/main' into qp_objective_checks
rg20 Sep 28, 2026
1ffe501
Remove default handling of IR tolerance
rg20 Sep 28, 2026
6f2eb49
Implement adaptive IR
rg20 Sep 28, 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
104 changes: 75 additions & 29 deletions cpp/src/barrier/barrier.cu
Original file line number Diff line number Diff line change
Expand Up @@ -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<i_t, f_t> diag;
pinned_dense_vector_t<i_t, f_t> inv_diag;
dense_vector_t<i_t, f_t> inv_sqrt_diag;
Expand Down Expand Up @@ -2581,7 +2583,8 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_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<i_t, f_t, op_t>(op, rhs, soln, ir_tol);
}

Expand Down Expand Up @@ -3175,7 +3178,8 @@ i_t barrier_solver_t<i_t, f_t>::gpu_compute_search_direction(iteration_data_t<i_
} op(data);
if (settings.barrier_iterative_refinement) {
raft::common::nvtx::range fun_scope("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);
const f_t solve_err = iterative_refinement<i_t, f_t, op_t>(
op, data.d_augmented_rhs_, data.d_augmented_soln_, ir_tol);
if (solve_err > 1e-1) {
Expand Down Expand Up @@ -3263,8 +3267,9 @@ i_t barrier_solver_t<i_t, f_t>::gpu_compute_search_direction(iteration_data_t<i_
data_.gpu_solve_adat(b, x);
}
} adat_op(data);
const f_t ir_tol = 1e-8;
const f_t adat_solve_err =
iterative_refinement<i_t, f_t, adat_op_t>(adat_op, data.d_h_, data.d_dy_);
iterative_refinement<i_t, f_t, adat_op_t>(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);
}
Expand Down Expand Up @@ -4266,16 +4271,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) ||

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

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 @@ -4300,24 +4309,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 =

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.

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?

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.

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

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 @@ -4340,14 +4364,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 @@ -4449,18 +4478,18 @@ lp_status_t barrier_solver_t<i_t, f_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;
Expand All @@ -4477,16 +4506,19 @@ lp_status_t barrier_solver_t<i_t, f_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) ||

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.

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.

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 @@ -4535,9 +4567,11 @@ lp_status_t barrier_solver_t<i_t, f_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) {
Expand Down Expand Up @@ -4575,9 +4609,11 @@ lp_status_t barrier_solver_t<i_t, f_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;
Expand Down Expand Up @@ -4605,17 +4641,20 @@ lp_status_t barrier_solver_t<i_t, f_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)));

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 (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 &&
Expand Down Expand Up @@ -4663,9 +4702,11 @@ lp_status_t barrier_solver_t<i_t, f_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);
}

Expand All @@ -4683,7 +4724,8 @@ lp_status_t barrier_solver_t<i_t, f_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) ||

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 @@ -4701,6 +4743,8 @@ lp_status_t barrier_solver_t<i_t, f_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);
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 @@ -4735,9 +4779,11 @@ lp_status_t barrier_solver_t<i_t, f_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);
}
}
Expand Down
2 changes: 2 additions & 0 deletions cpp/src/barrier/barrier.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<i_t, f_t>& solution);

const simplex::lp_problem_t<i_t, f_t>& lp;
Expand Down
6 changes: 3 additions & 3 deletions cpp/src/barrier/iterative_refinement.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ template <typename i_t, typename f_t, typename T>
f_t iterative_refinement_simple(T& op,
const rmm::device_uvector<f_t>& b,
rmm::device_uvector<f_t>& x,
f_t tol = 1e-8)
f_t tol)
{
rmm::device_uvector<f_t> x_sav(x, x.stream());

Expand Down Expand Up @@ -125,7 +125,7 @@ template <typename i_t, typename f_t, typename T>
f_t iterative_refinement_gmres(T& op,
const rmm::device_uvector<f_t>& b,
rmm::device_uvector<f_t>& 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
Expand Down Expand Up @@ -383,7 +383,7 @@ template <typename i_t, typename f_t, typename T>
f_t iterative_refinement(T& op,
const rmm::device_uvector<f_t>& b,
rmm::device_uvector<f_t>& x,
f_t tol = 1e-8)
f_t tol)
{
return iterative_refinement_gmres<i_t, f_t, T>(op, b, x, tol);
}
Expand Down
5 changes: 5 additions & 0 deletions cpp/src/dual_simplex/simplex_solver_settings.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<f_t>::infinity()),
steepest_edge_ratio(0.5),
steepest_edge_primal_tol(1e-9),
Expand Down Expand Up @@ -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
Expand Down
Loading
Loading