Estimate minimal eigenvalue of quadratic cost hessian

README.md

@@ -42,6 +42,7 @@ We are ready to integrate **ProxSuite** within other optimization ecosystems.
    - Python and Julia bindings for easy code prototyping without sacrificing performance.
 
 **Proxsuite** has a dedicated feature for solving batches of QPs.
+**Proxsuite** has a dedicated feature for solving nonconvex QPs.
 **Proxsuite** has a dedicated feature for solving the closest feasible QPs if they appear to be primal infeasible.
 **Proxsuite** is extensible.
 **Proxsuite** is reliable and extensively tested, showing the best performances on the hardest problems of the literature.

bindings/python/src/algorithms.hpp

@@ -11,6 +11,7 @@
 #include "expose-qpobject.hpp"
 #include "expose-qpvector.hpp"
 #include "expose-solve.hpp"
+#include "expose-helpers.hpp"
 #ifdef PROXSUITE_PYTHON_INTERFACE_WITH_OPENMP
 #include "expose-parallel.hpp"
 #endif

bindings/python/src/expose-all.cpp

@@ -35,6 +35,7 @@ exposeSparseAlgorithms(pybind11::module_ m)
   sparse::python::exposeQpObjectSparse<T, I>(m);
   sparse::python::exposeQPVectorSparse<T, I>(m);
   sparse::python::solveSparseQp<T, I>(m);
+  sparse::python::exposeSparseHelpers<T, I>(m);
 }
 
 template<typename T>
@@ -45,6 +46,7 @@ exposeDenseAlgorithms(pybind11::module_ m)
   dense::python::exposeQpObjectDense<T>(m);
   dense::python::exposeQPVectorDense<T>(m);
   dense::python::solveDenseQp<T>(m);
+  dense::python::exposeDenseHelpers<T>(m);
 }
 
 #ifdef PROXSUITE_PYTHON_INTERFACE_WITH_OPENMP

bindings/python/src/expose-helpers.hpp

@@ -0,0 +1,72 @@
+//
+// Copyright (c) 2022 INRIA
+//
+
+#include <pybind11/pybind11.h>
+#include <pybind11/eigen.h>
+#include <pybind11/stl.h>
+
+#include <proxsuite/proxqp/dense/helpers.hpp>
+#include <proxsuite/proxqp/sparse/helpers.hpp>
+
+namespace proxsuite {
+namespace proxqp {
+
+namespace dense {
+
+namespace python {
+
+template<typename T>
+void
+exposeDenseHelpers(pybind11::module_ m)
+{
+  m.def("estimate_minimal_eigen_value_of_symmetric_matrix",
+        &dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>,
+        "Function for estimating the minimal eigenvalue of a dense symmetric "
+        "matrix. "
+        "Two options are available: an exact method using "
+        "SelfAdjointEigenSolver from Eigen, "
+        "or a Power Iteration algorithm (with parameters : "
+        "power_iteration_accuracy and nb_power_iteration).",
+        pybind11::arg("H"),
+        pybind11::arg_v("estimate_method_option",
+                        EigenValueEstimateMethodOption::ExactMethod,
+                        "Two options are available for "
+                        "estimating smallest eigenvalue: either a power "
+                        "iteration algorithm, or an exact method from Eigen."),
+        pybind11::arg_v(
+          "power_iteration_accuracy", T(1.E-3), "power iteration accuracy."),
+        pybind11::arg_v("nb_power_iteration",
+                        1000,
+                        "maximal number of power iteration executed."));
+}
+} // namespace python
+} // namespace dense
+
+namespace sparse {
+
+namespace python {
+
+template<typename T, typename I>
+void
+exposeSparseHelpers(pybind11::module_ m)
+{
+  m.def("estimate_minimal_eigen_value_of_symmetric_matrix",
+        &sparse::estimate_minimal_eigen_value_of_symmetric_matrix<T, I>,
+        "Function for estimating the minimal eigenvalue of a sparse symmetric "
+        "matrix, "
+        " using aPower Iteration algorithm (with parameters : "
+        "power_iteration_accuracy and nb_power_iteration).",
+        pybind11::arg("H"),
+        pybind11::arg_v(
+          "power_iteration_accuracy", T(1.E-3), "power iteration accuracy."),
+        pybind11::arg_v("nb_power_iteration",
+                        1000,
+                        "maximal number of power iteration executed."));
+}
+
+} // namespace python
+} // namespace sparse
+
+} // namespace proxqp
+} // namespace proxsuite

bindings/python/src/expose-qpobject.hpp

@@ -88,6 +88,7 @@ exposeQpObjectDense(pybind11::module_ m)
                                             bool compute_preconditioner,
                                             optional<T>,
                                             optional<T>,
+                                            optional<T>,
                                             optional<T>)>(&dense::QP<T>::init),
          "function for initialize the QP model.",
          pybind11::arg_v("H", nullopt, "quadratic cost"),
@@ -105,7 +106,11 @@ exposeQpObjectDense(pybind11::module_ m)
          pybind11::arg_v(
            "mu_eq", nullopt, "dual equality constraint proximal parameter"),
          pybind11::arg_v(
-           "mu_in", nullopt, "dual inequality constraint proximal parameter"))
+           "mu_in", nullopt, "dual inequality constraint proximal parameter"),
+         pybind11::arg_v("manual_minimal_H_eigenvalue",
+                         nullopt,
+                         "manual minimal H eigenvalue proposed to regularize H"
+                         " in case it is non convex."))
     .def("init",
          static_cast<void (dense::QP<T>::*)(optional<dense::MatRef<T>>,
                                             optional<dense::VecRef<T>>,
@@ -119,6 +124,7 @@ exposeQpObjectDense(pybind11::module_ m)
                                             bool compute_preconditioner,
                                             optional<T>,
                                             optional<T>,
+                                            optional<T>,
                                             optional<T>)>(&dense::QP<T>::init),
          "function for initialize the QP model.",
          pybind11::arg_v("H", nullopt, "quadratic cost"),
@@ -140,7 +146,11 @@ exposeQpObjectDense(pybind11::module_ m)
          pybind11::arg_v(
            "mu_eq", nullopt, "dual equality constraint proximal parameter"),
          pybind11::arg_v(
-           "mu_in", nullopt, "dual inequality constraint proximal parameter"))
+           "mu_in", nullopt, "dual inequality constraint proximal parameter"),
+         pybind11::arg_v("manual_minimal_H_eigenvalue",
+                         nullopt,
+                         "manual minimal H eigenvalue proposed to regularize H"
+                         " in case it is non convex."))
     .def("solve",
          static_cast<void (dense::QP<T>::*)()>(&dense::QP<T>::solve),
          "function used for solving the QP problem, using default parameters.")
@@ -163,6 +173,7 @@ exposeQpObjectDense(pybind11::module_ m)
                                          bool update_preconditioner,
                                          optional<T>,
                                          optional<T>,
+                                         optional<T>,
                                          optional<T>)>(&dense::QP<T>::update),
       "function used for updating matrix or vector entry of the model using "
       "dense matrix entries.",
@@ -183,7 +194,11 @@ exposeQpObjectDense(pybind11::module_ m)
       pybind11::arg_v(
         "mu_eq", nullopt, "dual equality constraint proximal parameter"),
       pybind11::arg_v(
-        "mu_in", nullopt, "dual inequality constraint proximal parameter"))
+        "mu_in", nullopt, "dual inequality constraint proximal parameter"),
+      pybind11::arg_v("manual_minimal_H_eigenvalue",
+                      nullopt,
+                      "manual minimal H eigenvalue proposed to regularize H"
+                      " in case it is non convex."))
     .def(
       "update",
       static_cast<void (dense::QP<T>::*)(optional<dense::MatRef<T>>,
@@ -198,6 +213,7 @@ exposeQpObjectDense(pybind11::module_ m)
                                          bool update_preconditioner,
                                          optional<T>,
                                          optional<T>,
+                                         optional<T>,
                                          optional<T>)>(&dense::QP<T>::update),
       "function used for updating matrix or vector entry of the model using "
       "dense matrix entries.",
@@ -222,7 +238,11 @@ exposeQpObjectDense(pybind11::module_ m)
       pybind11::arg_v(
         "mu_eq", nullopt, "dual equality constraint proximal parameter"),
       pybind11::arg_v(
-        "mu_in", nullopt, "dual inequality constraint proximal parameter"))
+        "mu_in", nullopt, "dual inequality constraint proximal parameter"),
+      pybind11::arg_v("manual_minimal_H_eigenvalue",
+                      nullopt,
+                      "manual minimal H eigenvalue proposed to regularize H"
+                      " in case it is non convex."))
     .def("cleanup",
          &dense::QP<T>::cleanup,
          "function used for cleaning the workspace and result "
@@ -297,7 +317,11 @@ exposeQpObjectSparse(pybind11::module_ m)
          pybind11::arg_v(
            "mu_eq", nullopt, "dual equality constraint proximal parameter"),
          pybind11::arg_v(
-           "mu_in", nullopt, "dual inequality constraint proximal parameter"))
+           "mu_in", nullopt, "dual inequality constraint proximal parameter"),
+         pybind11::arg_v("manual_minimal_H_eigenvalue",
+                         nullopt,
+                         "manual minimal H eigenvalue proposed to regularize H"
+                         " in case it is non convex."))
 
     .def("update",
          &sparse::QP<T, I>::update,
@@ -319,7 +343,11 @@ exposeQpObjectSparse(pybind11::module_ m)
          pybind11::arg_v(
            "mu_eq", nullopt, "dual equality constraint proximal parameter"),
          pybind11::arg_v(
-           "mu_in", nullopt, "dual inequality constraint proximal parameter"))
+           "mu_in", nullopt, "dual inequality constraint proximal parameter"),
+         pybind11::arg_v("manual_minimal_H_eigenvalue",
+                         nullopt,
+                         "manual minimal H eigenvalue proposed to regularize H"
+                         " in case it is non convex."))
     .def("solve",
          static_cast<void (sparse::QP<T, I>::*)()>(&sparse::QP<T, I>::solve),
          "function used for solving the QP problem, using default parameters.")

bindings/python/src/expose-results.hpp

@@ -51,7 +51,12 @@ exposeResults(pybind11::module_ m)
     .def_readwrite("sparse_backend",
                    &Info<T>::sparse_backend,
                    "Sparse backend used to solve the qp, either SparseCholesky "
-                   "or MatrixFree.");
+                   "or MatrixFree.")
+    .def_readwrite("minimal_H_eigenvalue_estimate",
+                   &Info<T>::minimal_H_eigenvalue_estimate,
+                   "By default it equals 0, in order to get an estimate, set "
+                   "appropriately the setting option "
+                   "find_H_minimal_eigenvalue.");
 
   ::pybind11::class_<Results<T>>(m, "Results", pybind11::module_local())
     .def(::pybind11::init<i64, i64, i64>(),

bindings/python/src/expose-settings.hpp

@@ -41,6 +41,11 @@ exposeSettings(pybind11::module_ m)
     .value("MatrixFree", SparseBackend::MatrixFree)
     .value("SparseCholesky", SparseBackend::SparseCholesky)
     .export_values();
+  ::pybind11::enum_<EigenValueEstimateMethodOption>(
+    m, "EigenValueEstimateMethodOption", pybind11::module_local())
+    .value("PowerIteration", EigenValueEstimateMethodOption::PowerIteration)
+    .value("ExactMethod", EigenValueEstimateMethodOption::ExactMethod)
+    .export_values();
 
   ::pybind11::class_<Settings<T>>(m, "Settings", pybind11::module_local())
     .def(::pybind11::init(), "Default constructor.") // constructor
@@ -87,6 +92,8 @@ exposeSettings(pybind11::module_ m)
                    &Settings<T>::primal_infeasibility_solving)
     .def_readwrite("frequence_infeasibility_check",
                    &Settings<T>::frequence_infeasibility_check)
+    .def_readwrite("default_H_eigenvalue_estimate",
+                   &Settings<T>::default_H_eigenvalue_estimate)
     .def(pybind11::self == pybind11::self)
     .def(pybind11::self != pybind11::self)
     .def(pybind11::pickle(

bindings/python/src/expose-solve.hpp

@@ -43,7 +43,8 @@ solveDenseQp(pybind11::module_ m)
                             bool,
                             optional<T>,
                             optional<T>,
-                            bool>(&dense::solve<T>),
+                            bool,
+                            optional<T>>(&dense::solve<T>),
     "Function for solving a QP problem using PROXQP sparse backend directly "
     "without defining a QP object. It is possible to set up some of the solver "
     "parameters (warm start, initial guess option, proximal step sizes, "
@@ -103,7 +104,10 @@ solveDenseQp(pybind11::module_ m)
     pybind11::arg_v("primal_infeasibility_solving",
                     false,
                     "solves the closest feasible problem in L2 sense "
-                    "if the QP problem appears to be infeasible."));
+                    "if the QP problem appears to be infeasible."),
+    pybind11::arg_v("default_H_eigenvalue_estimate",
+                    0.,
+                    "Default estimate of the minimal eigen value of H."));
 
   m.def(
     "solve",
@@ -132,7 +136,8 @@ solveDenseQp(pybind11::module_ m)
                             bool,
                             optional<T>,
                             optional<T>,
-                            bool>(&dense::solve<T>),
+                            bool,
+                            optional<T>>(&dense::solve<T>),
     "Function for solving a QP problem using PROXQP sparse backend directly "
     "without defining a QP object. It is possible to set up some of the solver "
     "parameters (warm start, initial guess option, proximal step sizes, "
@@ -194,7 +199,10 @@ solveDenseQp(pybind11::module_ m)
     pybind11::arg_v("primal_infeasibility_solving",
                     false,
                     "solves the closest feasible problem in L2 sense "
-                    "if the QP problem appears to be infeasible."));
+                    "if the QP problem appears to be infeasible."),
+    pybind11::arg_v("default_H_eigenvalue_estimate",
+                    0.,
+                    "Default estimate of the minimal eigen value of H."));
 }
 
 } // namespace python
@@ -269,7 +277,10 @@ solveSparseQp(pybind11::module_ m)
     pybind11::arg_v("primal_infeasibility_solving",
                     false,
                     "solves the closest feasible problem in L2 sense "
-                    "if the QP problem appears to be infeasible."));
+                    "if the QP problem appears to be infeasible."),
+    pybind11::arg_v("default_H_eigenvalue_estimate",
+                    0.,
+                    "Default estimate of the minimal eigen value of H."));
 }
 
 } // namespace python