Merge pull request #114 from jcarpent/devel

Enforce robustness computation of duality gap

Author: jcarpent

include/proxsuite/helpers/common.hpp

@@ -24,6 +24,38 @@ struct infinite_bound
   }
 };
 
+#define PROXSUITE_DEDUCE_RET(...)                                              \
+  noexcept(noexcept(__VA_ARGS__))                                              \
+    ->typename std::remove_const<decltype(__VA_ARGS__)>::type                  \
+  {                                                                            \
+    return __VA_ARGS__;                                                        \
+  }                                                                            \
+  static_assert(true, ".")
+
+/// @brief \brief Returns the part of the expression which is lower than value
+template<typename T, typename Scalar>
+auto
+at_most(T const& expr, const Scalar value) PROXSUITE_DEDUCE_RET(
+  (expr.array() < value).select(expr, T::Constant(expr.rows(), value)));
+
+/// @brief \brief Returns the part of the expression which is greater than value
+template<typename T, typename Scalar>
+auto
+at_least(T const& expr, const Scalar value) PROXSUITE_DEDUCE_RET(
+  (expr.array() > value).select(expr, T::Constant(expr.rows(), value)));
+
+/// @brief \brief Returns the positive part of an expression
+template<typename T>
+auto
+positive_part(T const& expr)
+  PROXSUITE_DEDUCE_RET((expr.array() > 0).select(expr, T::Zero(expr.rows())));
+
+/// @brief \brief Returns the negative part of an expression
+template<typename T>
+auto
+negative_part(T const& expr)
+  PROXSUITE_DEDUCE_RET((expr.array() < 0).select(expr, T::Zero(expr.rows())));
+
 } // helpers
 } // proxsuite
 

include/proxsuite/proxqp/dense/solver.hpp

@@ -459,8 +459,8 @@ compute_inner_loop_saddle_point(const Model<T>& qpmodel,
 {
 
   qpwork.active_part_z =
-    positive_part(qpwork.primal_residual_in_scaled_up) +
-    negative_part(qpwork.primal_residual_in_scaled_low) -
+    helpers::positive_part(qpwork.primal_residual_in_scaled_up) +
+    helpers::negative_part(qpwork.primal_residual_in_scaled_low) -
     qpresults.z * qpresults.info.mu_in; // contains now : [Cx-u+z_prev*mu_in]+
                                         // + [Cx-l+z_prev*mu_in]- - z*mu_in
 

include/proxsuite/proxqp/dense/utils.hpp

@@ -7,17 +7,20 @@
 #ifndef PROXSUITE_PROXQP_DENSE_UTILS_HPP
 #define PROXSUITE_PROXQP_DENSE_UTILS_HPP
 
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <type_traits>
+
+#include "proxsuite/helpers/common.hpp"
 #include "proxsuite/proxqp/dense/views.hpp"
 #include "proxsuite/proxqp/dense/workspace.hpp"
 #include <proxsuite/proxqp/dense/model.hpp>
 #include <proxsuite/proxqp/results.hpp>
 #include <proxsuite/proxqp/utils/prints.hpp>
 #include <proxsuite/proxqp/settings.hpp>
 #include <proxsuite/proxqp/dense/preconditioner/ruiz.hpp>
-#include <iostream>
-#include <fstream>
-#include <cmath>
-#include <type_traits>
+
 //#include <fmt/format.h>
 //#include <fmt/ostream.h>
 
@@ -110,22 +113,6 @@ save_data(const std::string& filename, const ::Eigen::MatrixBase<Derived>& mat)
   }
 }
 
-#define LDLT_DEDUCE_RET(...)                                                   \
-  noexcept(noexcept(__VA_ARGS__))                                              \
-    ->typename std::remove_const<decltype(__VA_ARGS__)>::type                  \
-  {                                                                            \
-    return __VA_ARGS__;                                                        \
-  }                                                                            \
-  static_assert(true, ".")
-template<typename T>
-auto
-positive_part(T const& expr)
-  LDLT_DEDUCE_RET((expr.array() > 0).select(expr, T::Zero(expr.rows())));
-template<typename T>
-auto
-negative_part(T const& expr)
-  LDLT_DEDUCE_RET((expr.array() < 0).select(expr, T::Zero(expr.rows())));
-
 /*!
  * Derives the global primal residual of the QP problem.
  *
@@ -180,8 +167,8 @@ global_primal_residual(const Model<T>& qpmodel,
   primal_feasibility_in_rhs_0 = infty_norm(qpwork.primal_residual_in_scaled_up);
 
   qpwork.primal_residual_in_scaled_low =
-    positive_part(qpwork.primal_residual_in_scaled_up - qpmodel.u) +
-    negative_part(qpwork.primal_residual_in_scaled_up - qpmodel.l);
+    helpers::positive_part(qpwork.primal_residual_in_scaled_up - qpmodel.u) +
+    helpers::negative_part(qpwork.primal_residual_in_scaled_up - qpmodel.l);
   qpwork.primal_residual_eq_scaled -= qpmodel.b;
 
   primal_feasibility_in_lhs = infty_norm(qpwork.primal_residual_in_scaled_low);
@@ -236,8 +223,8 @@ global_primal_residual_infeasibility(
   ruiz.unscale_dual_residual_in_place(ATdy);
   ruiz.unscale_dual_residual_in_place(CTdz);
   T eq_inf = dy.to_eigen().dot(qpwork.b_scaled);
-  T in_inf = positive_part(dz.to_eigen()).dot(qpwork.u_scaled) -
-             positive_part(-dz.to_eigen()).dot(qpwork.l_scaled);
+  T in_inf = helpers::positive_part(dz.to_eigen()).dot(qpwork.u_scaled) -
+             helpers::negative_part(dz.to_eigen()).dot(qpwork.l_scaled);
   ruiz.unscale_dual_in_place_eq(dy);
   ruiz.unscale_dual_in_place_in(dz);
 
@@ -375,7 +362,7 @@ global_dual_residual(Results<T>& qpresults,
   ruiz.unscale_primal_in_place(VectorViewMut<T>{ from_eigen, qpresults.x });
   duality_gap = (qpmodel.g).dot(qpresults.x);
   rhs_duality_gap = std::abs(duality_gap);
-  T xHx = (qpwork.CTz).dot(qpresults.x);
+  const T xHx = (qpwork.CTz).dot(qpresults.x);
   duality_gap += xHx; // contains now xHx+g.Tx
   rhs_duality_gap = std::max(rhs_duality_gap, std::abs(xHx));
 
@@ -402,19 +389,29 @@ global_dual_residual(Results<T>& qpresults,
     VectorViewMut<T>{ from_eigen, qpwork.dual_residual_scaled });
 
   ruiz.unscale_dual_in_place_eq(VectorViewMut<T>{ from_eigen, qpresults.y });
-  T by = (qpmodel.b).dot(qpresults.y);
+  const T by = (qpmodel.b).dot(qpresults.y);
   rhs_duality_gap = std::max(rhs_duality_gap, std::abs(by));
   duality_gap += by;
   ruiz.scale_dual_in_place_eq(VectorViewMut<T>{ from_eigen, qpresults.y });
+
   ruiz.unscale_dual_in_place_in(VectorViewMut<T>{ from_eigen, qpresults.z });
-  T zu = positive_part(qpresults.z).dot(qpmodel.u);
+
+  const T zu =
+    helpers::positive_part(qpresults.z)
+      .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;
-  T zl = negative_part(qpresults.z).dot(qpmodel.l);
+
+  const T zl =
+    helpers::negative_part(qpresults.z)
+      .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;
+
   ruiz.scale_dual_in_place_in(VectorViewMut<T>{ from_eigen, qpresults.z });
+
   duality_gap /= sqrt_max_dim; // in order to get an a-dimensional duality gap
+  rhs_duality_gap /= sqrt_max_dim;
 }
 
 } // namespace dense

include/proxsuite/proxqp/dense/wrapper.hpp

@@ -58,8 +58,8 @@ object Qp.solve(); // solve the problem
         // Verify solution accuracy
         T pri_res = std::max(
                         (qp.A * Qp.results.x - qp.b).lpNorm<Eigen::Infinity>(),
-                        (proxqp::dense::positive_part(qp.C * Qp.results.x -
-qp.u) + proxqp::dense::negative_part(qp.C * Qp.results.x - qp.l))
+                        (helpers::positive_part(qp.C * Qp.results.x -
+qp.u) + helpers::negative_part(qp.C * Qp.results.x - qp.l))
                                         .lpNorm<Eigen::Infinity>());
         T dua_res = (qp.H * Qp.results.x + qp.g + qp.A.transpose() *
 Qp.results.y + qp.C.transpose() * Qp.results.z) .lpNorm<Eigen::Infinity>();

include/proxsuite/proxqp/sparse/solver.hpp

@@ -1136,8 +1136,8 @@ qp_solve(Results<T>& results,
           primal_residual_in_scaled_up += alpha * Cdx;
 
           T err_in = std::max({
-            (infty_norm(detail::negative_part(primal_residual_in_scaled_lo) +
-                        detail::positive_part(primal_residual_in_scaled_up) -
+            (infty_norm(helpers::negative_part(primal_residual_in_scaled_lo) +
+                        helpers::positive_part(primal_residual_in_scaled_up) -
                         results.info.mu_in * z_e)),
             (infty_norm(primal_residual_eq_scaled)),
             (infty_norm(dual_residual_scaled)),

include/proxsuite/proxqp/sparse/utils.hpp

@@ -6,6 +6,11 @@
 #ifndef PROXSUITE_PROXQP_SPARSE_UTILS_HPP
 #define PROXSUITE_PROXQP_SPARSE_UTILS_HPP
 
+#include <iostream>
+#include <Eigen/IterativeLinearSolvers>
+#include <unsupported/Eigen/IterativeSolvers>
+
+#include "proxsuite/helpers/common.hpp"
 #include <proxsuite/linalg/dense/core.hpp>
 #include <proxsuite/linalg/sparse/core.hpp>
 #include <proxsuite/linalg/sparse/factorize.hpp>
@@ -21,10 +26,6 @@
 #include "proxsuite/proxqp/sparse/preconditioner/ruiz.hpp"
 #include "proxsuite/proxqp/sparse/preconditioner/identity.hpp"
 
-#include <iostream>
-#include <Eigen/IterativeLinearSolvers>
-#include <unsupported/Eigen/IterativeSolvers>
-
 namespace proxsuite {
 namespace proxqp {
 namespace sparse {
@@ -98,15 +99,8 @@ print_setup_header(const Settings<T>& settings,
         << std::endl;
   }
 }
+
 namespace detail {
-template<typename T>
-auto
-positive_part(T const& expr)
-  VEG_DEDUCE_RET((expr.array() > 0).select(expr, T::Zero(expr.rows())));
-template<typename T>
-auto
-negative_part(T const& expr)
-  VEG_DEDUCE_RET((expr.array() < 0).select(expr, T::Zero(expr.rows())));
 
 template<typename T, typename I>
 VEG_NO_INLINE void
@@ -500,8 +494,8 @@ global_primal_residual_infeasibility(VectorViewMut<T> ATdy,
   ruiz.unscale_dual_residual_in_place(ATdy);
   ruiz.unscale_dual_residual_in_place(CTdz);
   T eq_inf = dy.to_eigen().dot(qp_scaled.b.to_eigen());
-  T in_inf = positive_part(dz.to_eigen()).dot(qp_scaled.u.to_eigen()) -
-             positive_part(-dz.to_eigen()).dot(qp_scaled.l.to_eigen());
+  T in_inf = helpers::positive_part(dz.to_eigen()).dot(qp_scaled.u.to_eigen()) -
+             helpers::negative_part(dz.to_eigen()).dot(qp_scaled.l.to_eigen());
   ruiz.unscale_dual_in_place_eq(dy);
   ruiz.unscale_dual_in_place_in(dz);
 
@@ -657,27 +651,36 @@ unscaled_primal_dual_residual(
     precond.unscale_primal_in_place({ proxqp::from_eigen, x_e });
     results.info.duality_gap = x_e.dot(data.g); // contains gTx
     rhs_duality_gap = std::abs(results.info.duality_gap);
-    T xHx = (tmp).dot(x_e);
+
+    const T xHx = (tmp).dot(x_e);
     results.info.duality_gap += xHx;
     rhs_duality_gap = std::max(rhs_duality_gap, std::abs(xHx));
     tmp += data.g; // contains now Hx+g
     precond.scale_primal_in_place({ proxqp::from_eigen, x_e });
 
     precond.unscale_dual_in_place_eq({ proxsuite::proxqp::from_eigen, y_e });
-    T by = (data.b).dot(y_e);
+    const T by = (data.b).dot(y_e);
     results.info.duality_gap += by;
     rhs_duality_gap = std::max(rhs_duality_gap, std::abs(by));
     precond.scale_dual_in_place_eq({ proxsuite::proxqp::from_eigen, y_e });
+
     precond.unscale_dual_in_place_in({ proxsuite::proxqp::from_eigen, z_e });
-    T zl = (data.l).dot(detail::negative_part(z_e));
+
+    const T zl = helpers::negative_part(z_e).dot(
+      helpers::at_least(data.l, -helpers::infinite_bound<T>::value()));
     results.info.duality_gap += zl;
     rhs_duality_gap = std::max(rhs_duality_gap, std::abs(zl));
-    T zu = (data.u).dot(detail::positive_part(z_e));
+
+    const T zu = helpers::positive_part(z_e).dot(
+      helpers::at_most(data.u, helpers::infinite_bound<T>::value()));
     results.info.duality_gap += zu;
     rhs_duality_gap = std::max(rhs_duality_gap, std::abs(zu));
+
     precond.scale_dual_in_place_in({ proxsuite::proxqp::from_eigen, z_e });
+
     results.info.duality_gap /=
       sqrt_max_dim; // in order to get an a-dimensional duality gap
+    rhs_duality_gap /= sqrt_max_dim;
   }
 
   {
@@ -719,8 +722,8 @@ unscaled_primal_dual_residual(
   auto l = data.l;
   auto u = data.u;
   primal_residual_in_scaled_lo =
-    positive_part(primal_residual_in_scaled_up - u) +
-    negative_part(primal_residual_in_scaled_up - l);
+    helpers::positive_part(primal_residual_in_scaled_up - u) +
+    helpers::negative_part(primal_residual_in_scaled_up - l);
 
   primal_residual_eq_scaled -= b;
   T primal_feasibility_eq_lhs = infty_norm(primal_residual_eq_scaled);

include/proxsuite/proxqp/sparse/wrapper.hpp

@@ -68,8 +68,8 @@ ::proxsuite::proxqp::test::rand::vector_rand<T>(n_in); auto u = (l.array() +
         // Verify solution accuracy
         T pri_res = std::max(
                         (qp.A * Qp.results.x - qp.b).lpNorm<Eigen::Infinity>(),
-                        (proxqp::dense::positive_part(qp.C * Qp.results.x -
-qp.u) + proxqp::dense::negative_part(qp.C * Qp.results.x - qp.l))
+                        (helpers::positive_part(qp.C * Qp.results.x -
+qp.u) + helpers::negative_part(qp.C * Qp.results.x - qp.l))
                                         .lpNorm<Eigen::Infinity>());
         T dua_res = (qp.H * Qp.results.x + qp.g + qp.A.transpose() *
 Qp.results.y + qp.C.transpose() * Qp.results.z) .lpNorm<Eigen::Infinity>();

include/proxsuite/proxqp/utils/prints.hpp

@@ -31,10 +31,10 @@ print_preambule()
 {
   print_line();
   std::cout
-    << "                              ProxQP  -  Primal Dual Proximal QP "
+    << "                              ProxQP - Primal-Dual Proximal QP "
        "Solver\n"
     << "     (c) Antoine Bambade, Sarah El Kazdadi, Fabian Schramm, Adrien "
-       "Taylor, "
+       "Taylor, and "
        "Justin Carpentier\n"
     << "                                         Inria Paris 2022        \n"
     << std::endl;

test/src/cvxpy.cpp

@@ -42,8 +42,8 @@ DOCTEST_TEST_CASE("3 dim test case from cvxpy, check feasibility")
   Results<T> results = dense::solve<T>(
     H, g, nullopt, nullopt, C, l, u, nullopt, nullopt, nullopt, eps_abs, 0);
 
-  T pri_res = (dense::positive_part(C * results.x - u) +
-               dense::negative_part(C * results.x - l))
+  T pri_res = (helpers::positive_part(C * results.x - u) +
+               helpers::negative_part(C * results.x - l))
                 .lpNorm<Eigen::Infinity>();
   T dua_res =
     (H * results.x + g + C.transpose() * results.z).lpNorm<Eigen::Infinity>();
@@ -82,8 +82,8 @@ DOCTEST_TEST_CASE("simple test case from cvxpy, check feasibility")
   Results<T> results = dense::solve<T>(
     H, g, nullopt, nullopt, C, l, u, nullopt, nullopt, nullopt, eps_abs, 0);
 
-  T pri_res = (dense::positive_part(C * results.x - u) +
-               dense::negative_part(C * results.x - l))
+  T pri_res = (helpers::positive_part(C * results.x - u) +
+               helpers::negative_part(C * results.x - l))
                 .lpNorm<Eigen::Infinity>();
   T dua_res =
     (H * results.x + g + C.transpose() * results.z).lpNorm<Eigen::Infinity>();
@@ -140,8 +140,8 @@ DOCTEST_TEST_CASE("simple test case from cvxpy, init with solution, check that "
   z << 0.0;
   qp.solve(x, nullopt, z);
 
-  T pri_res = (dense::positive_part(C * qp.results.x - u) +
-               dense::negative_part(C * qp.results.x - l))
+  T pri_res = (helpers::positive_part(C * qp.results.x - u) +
+               helpers::negative_part(C * qp.results.x - l))
                 .lpNorm<Eigen::Infinity>();
   T dua_res = (H * qp.results.x + g + C.transpose() * qp.results.z)
                 .lpNorm<Eigen::Infinity>();

test/src/dense_maros_meszaros.cpp

@@ -136,8 +136,8 @@ TEST_CASE("dense maros meszaros using the api")
 
         T prim_eq = proxqp::dense::infty_norm(A * x - b);
         T prim_in =
-          proxqp::dense::infty_norm(proxqp::dense::positive_part(C * x - u) +
-                                    proxqp::dense::negative_part(C * x - l));
+          proxqp::dense::infty_norm(helpers::positive_part(C * x - u) +
+                                    helpers::negative_part(C * x - l));
         std::cout << "primal residual " << std::max(prim_eq, prim_in)
                   << std::endl;
         std::cout << "dual residual "