Merge pull request #163 from jcarpent/devel

Fix computation of duality gap quantity

Author: jcarpent

LICENSE

@@ -1,6 +1,6 @@
 BSD 2-Clause License
 
-Copyright (c) 2022, Inria
+Copyright (c) 2022-2023, Inria
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without

README.md

@@ -9,7 +9,7 @@
   <a href="https://anaconda.org/conda-forge/proxsuite"><img src="https://anaconda.org/conda-forge/proxsuite/badges/version.svg" alt="Conda version" height="20"></a>
 </p>
 
-**ProxSuite** is a collection of open-source, numerically robust, precise and efficient numerical solvers (e.g., LPs, QPs, etc.) rooted in revisited primal-dual proximal algorithms. 
+**ProxSuite** is a collection of open-source, numerically robust, precise and efficient numerical solvers (e.g., LPs, QPs, etc.) rooted in revisited primal-dual proximal algorithms.
 Through **ProxSuite**, we aim to offer the community scalable optimizers that can deal with dense, sparse or matrix-free problems. While the first targeted application is Robotics, **ProxSuite** can be used in other contexts without limits.
 
 **ProxSuite** is actively developped and supported by the [Willow](https://www.di.ens.fr/willow/) and [Sierra](https://www.di.ens.fr/sierra/) research groups, joint research teams between [Inria](https://www.inria.fr/en), [École Normale Supérieure de Paris](https://www.ens.fr) and [Centre National de la Recherche Scientifique](https://www.cnrs.fr) localized in France.

bindings/python/src/expose-results.hpp

@@ -1,5 +1,5 @@
 //
-// Copyright (c) 2022 INRIA
+// Copyright (c) 2022-2023 INRIA
 //
 #include <proxsuite/proxqp/results.hpp>
 #include <pybind11/pybind11.h>
@@ -39,6 +39,7 @@ exposeResults(pybind11::module_ m)
     .def_readwrite("run_time", &Info<T>::run_time)
     .def_readwrite("setup_time", &Info<T>::setup_time)
     .def_readwrite("solve_time", &Info<T>::solve_time)
+    .def_readwrite("duality_gap", &Info<T>::duality_gap)
     .def_readwrite("pri_res", &Info<T>::pri_res)
     .def_readwrite("dua_res", &Info<T>::dua_res)
     .def_readwrite("objValue", &Info<T>::objValue)

include/proxsuite/helpers/common.hpp

@@ -56,6 +56,14 @@ auto
 negative_part(T const& expr)
   PROXSUITE_DEDUCE_RET((expr.array() < 0).select(expr, T::Zero(expr.rows())));
 
+/// @brief \brief Select the components of the expression if the condition is
+/// fullfiled. Otherwise, set the component to value
+template<typename Condition, typename T, typename Scalar>
+auto
+select(Condition const& condition, T const& expr, const Scalar value)
+  PROXSUITE_DEDUCE_RET((condition).select(expr,
+                                          T::Constant(expr.rows(), value)));
+
 } // helpers
 } // proxsuite
 

include/proxsuite/proxqp/dense/solver.hpp

@@ -1065,34 +1065,24 @@ qp_solve( //
                 << std::endl;
       std::cout << std::scientific << std::setw(2) << std::setprecision(2)
                 << "| primal residual=" << qpresults.info.pri_res
-                << "| dual residual=" << qpresults.info.dua_res
-                << "| duality gap=" << qpresults.info.duality_gap
-                << "| mu_in=" << qpresults.info.mu_in
-                << "| rho=" << qpresults.info.rho << std::endl;
+                << " | dual residual=" << qpresults.info.dua_res
+                << " | duality gap=" << qpresults.info.duality_gap
+                << " | mu_in=" << qpresults.info.mu_in
+                << " | rho=" << qpresults.info.rho << std::endl;
       ruiz.scale_primal_in_place(VectorViewMut<T>{ from_eigen, qpresults.x });
       ruiz.scale_dual_in_place_eq(VectorViewMut<T>{ from_eigen, qpresults.y });
       ruiz.scale_dual_in_place_in(VectorViewMut<T>{ from_eigen, qpresults.z });
     }
-    if (is_primal_feasible) {
-
-      if (dual_feasibility_lhs >=
-            qpsettings.refactor_dual_feasibility_threshold &&
-          qpresults.info.rho != qpsettings.refactor_rho_threshold) {
-
-        T rho_new(qpsettings.refactor_rho_threshold);
-
-        refactorize(qpmodel, qpresults, qpwork, rho_new);
-        qpresults.info.rho_updates += 1;
-
-        qpresults.info.rho = rho_new;
-      }
-      if (is_dual_feasible) {
-        if (qpresults.info.duality_gap <=
-            qpsettings.eps_abs +
-              (qpsettings.eps_abs + qpsettings.eps_rel) * rhs_duality_gap) {
+    if (is_primal_feasible && is_dual_feasible) {
+      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;
       }
     }
     qpresults.info.iter_ext += 1; // We start a new external loop update
@@ -1170,9 +1160,13 @@ qp_solve( //
              std::max(dual_feasibility_rhs_1, qpwork.dual_feasibility_rhs_2)));
 
       if (is_dual_feasible) {
-        if (qpresults.info.duality_gap <=
-            qpsettings.eps_abs +
-              (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/dense/utils.hpp

@@ -397,13 +397,13 @@ global_dual_residual(Results<T>& qpresults,
   ruiz.unscale_dual_in_place_in(VectorViewMut<T>{ from_eigen, qpresults.z });
 
   const T zu =
-    helpers::positive_part(qpresults.z)
+    helpers::select(qpwork.active_set_up, qpresults.z, 0)
       .dot(helpers::at_most(qpmodel.u, helpers::infinite_bound<T>::value()));
   rhs_duality_gap = std::max(rhs_duality_gap, std::abs(zu));
   duality_gap += zu;
 
   const T zl =
-    helpers::negative_part(qpresults.z)
+    helpers::select(qpwork.active_set_low, qpresults.z, 0)
       .dot(helpers::at_least(qpmodel.l, -helpers::infinite_bound<T>::value()));
   rhs_duality_gap = std::max(rhs_duality_gap, std::abs(zl));
   duality_gap += zl;

include/proxsuite/proxqp/results.hpp

@@ -69,7 +69,6 @@ struct Results
   sparse::Vec<T> x;
   sparse::Vec<T> y;
   sparse::Vec<T> z;
-  proxsuite::linalg::veg::Vec<bool> active_constraints;
 
   Info<T> info;
 

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/fwd.hpp

@@ -1,5 +1,5 @@
 //
-// Copyright (c) 2022 INRIA
+// Copyright (c) 2022-2023 INRIA
 //
 /** \file */
 #ifndef PROXSUITE_PROXQP_SPARSE_FWD_HPP
@@ -40,6 +40,16 @@ using Vec = Eigen::Matrix<T, DYN, 1>;
 template<typename T, typename I>
 using Mat = Eigen::SparseMatrix<T, Eigen::ColMajor, I>;
 // using Mat = Eigen::SparseMatrix<T, Eigen::RowMajor, I>;
+using VecBool = Eigen::Matrix<bool, DYN, 1>;
+
+///
+/// @brief This class defines the workspace of the sparse solver.
+///
+/*!
+ * Workspace class of the sparse solver.
+ */
+template<typename T, typename I>
+struct Workspace;
 
 } // namespace sparse
 } // namespace proxqp

include/proxsuite/proxqp/sparse/helpers.hpp

@@ -148,10 +148,22 @@ qp_setup(QpView<T, I> qp,
     results.z.resize(n_in);
     results.z.setZero();
   }
-  if (results.active_constraints.len() != n_in) {
-    results.active_constraints.resize(n_in);
+  if (work.active_inequalities.len() != n_in) {
+    work.active_inequalities.resize(n_in);
     for (isize i = 0; i < n_in; ++i) {
-      results.active_constraints[i] = false;
+      work.active_inequalities[i] = false;
+    }
+  }
+  if (work.active_set_up.rows() != n_in) {
+    work.active_set_up.resize(n_in);
+    for (isize i = 0; i < n_in; ++i) {
+      work.active_set_up[i] = false;
+    }
+  }
+  if (work.active_set_low.rows() != n_in) {
+    work.active_set_low.resize(n_in);
+    for (isize i = 0; i < n_in; ++i) {
+      work.active_set_low[i] = false;
     }
   }
   bool execute_preconditioner_or_not = false;