duality gap as optional stopping criteria

- use settings to decide whether or not to check the duality gap

Author: fabinsch

include/proxsuite/proxqp/dense/solver.hpp

@@ -1074,8 +1074,13 @@ qp_solve( //
       ruiz.scale_dual_in_place_in(VectorViewMut<T>{ from_eigen, qpresults.z });
     }
     if (is_primal_feasible && is_dual_feasible) {
-      if (std::abs(qpresults.info.duality_gap) <=
-          qpsettings.eps_abs + qpsettings.eps_rel * rhs_duality_gap) {
+      if (qpsettings.check_duality_gap) {
+        if (std::abs(qpresults.info.duality_gap) <=
+            qpsettings.eps_abs + qpsettings.eps_rel * rhs_duality_gap) {
+          qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
+          break;
+        }
+      } else {
         qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
         break;
       }
@@ -1155,8 +1160,13 @@ qp_solve( //
              std::max(dual_feasibility_rhs_1, qpwork.dual_feasibility_rhs_2)));
 
       if (is_dual_feasible) {
-        if (std::abs(qpresults.info.duality_gap) <=
-            qpsettings.eps_abs + qpsettings.eps_rel * rhs_duality_gap) {
+        if (qpsettings.check_duality_gap) {
+          if (std::abs(qpresults.info.duality_gap) <=
+              qpsettings.eps_abs + qpsettings.eps_rel * rhs_duality_gap) {
+            qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
+            break;
+          }
+        } else {
           qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
           break;
         }

include/proxsuite/proxqp/settings.hpp

@@ -84,6 +84,7 @@ struct Settings
   bool update_preconditioner;
   bool compute_preconditioner;
   bool compute_timings;
+  bool check_duality_gap;
 
   isize preconditioner_max_iter;
   T preconditioner_accuracy;
@@ -136,6 +137,8 @@ struct Settings
    * @param compute_timings If set to true, timings will be computed by the
    * solver (setup time, solving time, and run time = setup time + solving
    * time).
+   * @param check_duality_gap If set to true, duality gap will be calculated and
+   * included in the stopping criterion.
    * @param preconditioner_max_iter maximal number of authorized iterations for
    * the preconditioner.
    * @param preconditioner_accuracy accuracy level of the preconditioner.
@@ -183,6 +186,7 @@ struct Settings
     bool update_preconditioner = true,
     bool compute_preconditioner = true,
     bool compute_timings = false,
+    bool check_duality_gap = false,
     isize preconditioner_max_iter = 10,
     T preconditioner_accuracy = 1.e-3,
     T eps_primal_inf = 1.E-4,
@@ -218,6 +222,7 @@ struct Settings
     , update_preconditioner(update_preconditioner)
     , compute_preconditioner(compute_preconditioner)
     , compute_timings(compute_timings)
+    , check_duality_gap(check_duality_gap)
     , preconditioner_max_iter(preconditioner_max_iter)
     , preconditioner_accuracy(preconditioner_accuracy)
     , eps_primal_inf(eps_primal_inf)
@@ -263,6 +268,7 @@ operator==(const Settings<T>& settings1, const Settings<T>& settings2)
     settings1.update_preconditioner == settings2.update_preconditioner &&
     settings1.compute_preconditioner == settings2.compute_preconditioner &&
     settings1.compute_timings == settings2.compute_timings &&
+    settings1.check_duality_gap == settings2.check_duality_gap &&
     settings1.preconditioner_max_iter == settings2.preconditioner_max_iter &&
     settings1.preconditioner_accuracy == settings2.preconditioner_accuracy &&
     settings1.eps_primal_inf == settings2.eps_primal_inf &&

include/proxsuite/proxqp/sparse/solver.hpp

@@ -842,8 +842,15 @@ qp_solve(Results<T>& results,
       }
       if (is_primal_feasible(primal_feasibility_lhs) &&
           is_dual_feasible(dual_feasibility_lhs)) {
-        if (std::abs(results.info.duality_gap) <=
-            settings.eps_abs + settings.eps_rel * rhs_duality_gap) {
+        if (settings.check_duality_gap) {
+          if (std::abs(results.info.duality_gap) <=
+              settings.eps_abs + settings.eps_rel * rhs_duality_gap) {
+            results.info.pri_res = primal_feasibility_lhs;
+            results.info.dua_res = dual_feasibility_lhs;
+            results.info.status = QPSolverOutput::PROXQP_SOLVED;
+            break;
+          }
+        } else {
           results.info.pri_res = primal_feasibility_lhs;
           results.info.dua_res = dual_feasibility_lhs;
           results.info.status = QPSolverOutput::PROXQP_SOLVED;
@@ -1235,8 +1242,15 @@ qp_solve(Results<T>& results,
 
       if (is_primal_feasible(primal_feasibility_lhs_new) &&
           is_dual_feasible(dual_feasibility_lhs_new)) {
-        if (std::fabs(results.info.duality_gap) <=
-            settings.eps_abs + settings.eps_rel * rhs_duality_gap) {
+        if (settings.check_duality_gap) {
+          if (std::fabs(results.info.duality_gap) <=
+              settings.eps_abs + settings.eps_rel * rhs_duality_gap) {
+            results.info.pri_res = primal_feasibility_lhs_new;
+            results.info.dua_res = dual_feasibility_lhs_new;
+            results.info.status = QPSolverOutput::PROXQP_SOLVED;
+            break;
+          }
+        } else {
           results.info.pri_res = primal_feasibility_lhs_new;
           results.info.dua_res = dual_feasibility_lhs_new;
           results.info.status = QPSolverOutput::PROXQP_SOLVED;

include/proxsuite/serialization/settings.hpp

@@ -45,6 +45,7 @@ serialize(Archive& archive, proxsuite::proxqp::Settings<T>& settings)
           CEREAL_NVP(settings.update_preconditioner),
           CEREAL_NVP(settings.compute_preconditioner),
           CEREAL_NVP(settings.compute_timings),
+          CEREAL_NVP(settings.check_duality_gap),
           CEREAL_NVP(settings.preconditioner_max_iter),
           CEREAL_NVP(settings.preconditioner_accuracy),
           CEREAL_NVP(settings.eps_primal_inf),