simplifying API logic for estimating minimal H eigenvalue

Author: Bambade

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-settings.hpp

@@ -41,13 +41,10 @@ exposeSettings(pybind11::module_ m)
     .value("MatrixFree", SparseBackend::MatrixFree)
     .value("SparseCholesky", SparseBackend::SparseCholesky)
     .export_values();
-  ::pybind11::enum_<HessianCostRegularization>(
-    m, "HessianCostRegularization", pybind11::module_local())
-    .value("NoRegularization", HessianCostRegularization::NoRegularization)
-    .value("Manual", HessianCostRegularization::Manual)
-    .value("PowerIteration", HessianCostRegularization::PowerIteration)
-    .value("EigenRegularization",
-           HessianCostRegularization::EigenRegularization)
+  ::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())
@@ -95,15 +92,8 @@ exposeSettings(pybind11::module_ m)
                    &Settings<T>::primal_infeasibility_solving)
     .def_readwrite("frequence_infeasibility_check",
                    &Settings<T>::frequence_infeasibility_check)
-    .def_readwrite("nb_power_iteration", &Settings<T>::nb_power_iteration)
-    .def_readwrite("power_iteration_accuracy",
-                   &Settings<T>::power_iteration_accuracy)
-    .def_readwrite("find_minimal_H_eigenvalue",
-                   &Settings<T>::find_minimal_H_eigenvalue)
     .def_readwrite("default_H_eigenvalue_estimate",
                    &Settings<T>::default_H_eigenvalue_estimate)
-    .def_readwrite("rho_regularization_scaling",
-                   &Settings<T>::rho_regularization_scaling)
     .def(pybind11::self == pybind11::self)
     .def(pybind11::self != pybind11::self)
     .def(pybind11::pickle(

bindings/python/src/expose-solve.hpp

@@ -44,10 +44,6 @@ solveDenseQp(pybind11::module_ m)
                             optional<T>,
                             optional<T>,
                             bool,
-                            optional<isize>,
-                            optional<T>,
-                            HessianCostRegularization,
-                            optional<T>,
                             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 "
@@ -109,32 +105,9 @@ solveDenseQp(pybind11::module_ m)
                     false,
                     "solves the closest feasible problem in L2 sense "
                     "if the QP problem appears to be infeasible."),
-    pybind11::arg_v("nb_power_iteration",
-                    1000,
-                    "Number of power iteration iteration used by default "
-                    "for estimating the lowest eigenvalue of H."),
-    pybind11::arg_v("power_iteration_accuracy",
-                    1.E-6,
-                    "Accuracy target of the power iteration algorithm "
-                    "for estimating the lowest eigenvalue of H."),
-    pybind11::arg_v("find_minimal_H_eigenvalue",
-                    HessianCostRegularization::NoRegularization,
-                    "Option for estimating the minimal eigen value of H "
-                    "and regularizing default_rho following "
-                    "default_rho=rho_regularization_scaling*abs(default_H_"
-                    "eigenvalue_estimate)."
-                    "This option can be used for solving non convex QPs."),
     pybind11::arg_v("default_H_eigenvalue_estimate",
                     0.,
-                    "Default estimate of the minimal eigen value of H."),
-    pybind11::arg_v("rho_regularization_scaling",
-                    1.5,
-                    "Scaling for regularizing default_rho according to the "
-                    "minimal eigen value of H "
-                    "following "
-                    "default_rho=rho_regularization_scaling*abs(default_H_"
-                    "eigenvalue_estimate)."
-                    "This option can be used for solving non convex QPs."));
+                    "Default estimate of the minimal eigen value of H."));
 
   m.def(
     "solve",
@@ -164,10 +137,6 @@ solveDenseQp(pybind11::module_ m)
                             optional<T>,
                             optional<T>,
                             bool,
-                            optional<isize>,
-                            optional<T>,
-                            HessianCostRegularization,
-                            optional<T>,
                             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 "
@@ -231,32 +200,9 @@ solveDenseQp(pybind11::module_ m)
                     false,
                     "solves the closest feasible problem in L2 sense "
                     "if the QP problem appears to be infeasible."),
-    pybind11::arg_v("nb_power_iteration",
-                    1000,
-                    "Number of power iteration iteration used by default "
-                    "for estimating the lowest eigenvalue of H."),
-    pybind11::arg_v("power_iteration_accuracy",
-                    1.E-6,
-                    "Accuracy target of the power iteration algorithm "
-                    "for estimating the lowest eigenvalue of H."),
-    pybind11::arg_v("find_minimal_H_eigenvalue",
-                    HessianCostRegularization::NoRegularization,
-                    "Option for estimating the minimal eigen value of H "
-                    "and regularizing default_rho following "
-                    "default_rho=rho_regularization_scaling*abs(default_H_"
-                    "eigenvalue_estimate)."
-                    "This option can be used for solving non convex QPs."),
     pybind11::arg_v("default_H_eigenvalue_estimate",
                     0.,
-                    "Default estimate of the minimal eigen value of H."),
-    pybind11::arg_v("rho_regularization_scaling",
-                    1.5,
-                    "Scaling for regularizing default_rho according to the "
-                    "minimal eigen value of H "
-                    "following "
-                    "default_rho=rho_regularization_scaling*abs(default_H_"
-                    "eigenvalue_estimate)."
-                    "This option can be used for solving non convex QPs."));
+                    "Default estimate of the minimal eigen value of H."));
 }
 
 } // namespace python
@@ -332,32 +278,9 @@ solveSparseQp(pybind11::module_ m)
                     false,
                     "solves the closest feasible problem in L2 sense "
                     "if the QP problem appears to be infeasible."),
-    pybind11::arg_v("nb_power_iteration",
-                    1000,
-                    "Number of power iteration iteration used by default "
-                    "for estimating the lowest eigenvalue of H."),
-    pybind11::arg_v("power_iteration_accuracy",
-                    1.E-6,
-                    "Accuracy target of the power iteration algorithm "
-                    "for estimating the lowest eigenvalue of H."),
-    pybind11::arg_v("find_minimal_H_eigenvalue",
-                    HessianCostRegularization::NoRegularization,
-                    "Option for estimating the minimal eigen value of H "
-                    "and regularizing default_rho following "
-                    "default_rho=rho_regularization_scaling*abs(default_H_"
-                    "eigenvalue_estimate)."
-                    "This option can be used for solving non convex QPs."),
     pybind11::arg_v("default_H_eigenvalue_estimate",
                     0.,
-                    "Default estimate of the minimal eigen value of H."),
-    pybind11::arg_v("rho_regularization_scaling",
-                    1.5,
-                    "Scaling for regularizing default_rho according to the "
-                    "minimal eigen value of H "
-                    "following "
-                    "default_rho=rho_regularization_scaling*abs(default_H_"
-                    "eigenvalue_estimate)."
-                    "This option can be used for solving non convex QPs."));
+                    "Default estimate of the minimal eigen value of H."));
 }
 
 } // namespace python

doc/2-PROXQP_API/2-ProxQP_api.md

@@ -415,7 +415,7 @@ In this table, you have the three columns from left to right: the name of the se
 | nb_power_iteration                  | 1000                               | Number of power iteration iteration used by default for estimating H lowest eigenvalue.
 | power_iteration_accuracy            | 1.E-6                              | If set to true, it solves the closest primal feasible problem if primal infeasibility is detected.
 | primal_infeasibility_solving        | False                              | Accuracy target of the power iteration algorithm for estimating the lowest eigenvalue of H.
-| find_minimal_H_eigenvalue           | NoRegularization                   | Option for estimating the minimal eigen value of H and regularizing default_rho  default_rho=rho_regularization_scaling*abs(default_H_eigenvalue_estimate). This option can be used for solving non convex QPs.
+| estimate_method_option           | NoRegularization                   | Option for estimating the minimal eigen value of H and regularizing default_rho  default_rho=rho_regularization_scaling*abs(default_H_eigenvalue_estimate). This option can be used for solving non convex QPs.
 | default_H_eigenvalue_estimate       | 0.                                 | Default estimate of the minimal eigen value of H.
 | rho_regularization_scaling          | 1.5                                | Scaling for regularizing default_rho according to the minimal eigen value of H.
 
@@ -436,15 +436,13 @@ If set to this option, the solver will start with no initial guess, which means
 
 \subsection OverviewEstimatingHminimalEigenValue The different options for estimating H minimal Eigenvalue
 
-The solver has four options for estimating the minimal eigenvalue of H within the struct HessianCostRegularization:
-* NoRegularization : set by default, it means the solver does not try to estimate it,
-* Manual: the user can provide an estimate of it through the init method,
+The solver environment provides an independent function for estimating the minimal eigenvalue of a dense or sparse symmetric matrix. It is named "estimate_minimal_eigen_value_of_symmetric_matrix". In the sparse case, it uses a power iteration algorithm (with two options: the maximal number of iterations and the accuracy target for the estimate). In the dense case, we provide two options within the struct EigenValueEstimateMethodOption:
 * PowerIteration: a power iteration algorithm will be used for estimating H minimal eigenvalue,
-* EigenRegularization: in case the dense backend is used, the solver make use of Eigen method for estimating it.
+* ExactMethod: in this case, an exact method from EigenSolver is used to provide an estimate.
 
-This option is particularly usefull when solving QP with non convex quadratics. Indeed, if default_rho is set to a value strictly higher than the minimal eigenvalue of H, then ProxQP is guaranteed for find a local minimum to the problem since it relies on a Proximal Method of Multipliers (for more detail for example this [work](https://arxiv.org/pdf/2010.02653.pdf) providing convergence proof of this property).
+Estimating minimal eigenvalue is particularly usefull for solving QP with non convex quadratics. Indeed, if default_rho is set to a value strictly higher than the minimal eigenvalue of H, then ProxQP is guaranteed for find a local minimum to the problem since it relies on a Proximal Method of Multipliers (for more detail for example this [work](https://arxiv.org/pdf/2010.02653.pdf) providing convergence proof of this property).
 
-More precisely, when HessianCostRegularization is set to a value different of NoRegularization, then ProxQP first estimate a minimal eigenvalue for H and then update default_rho following the rule: default_rho = rho_regularization_scaling * abs(default_H_eigenvalue_estimate), which guarantees for appropriate scaling than the proximal step-size is larger than the minimal eigenvalue of H. We provide below examples in C++ and python for using this feature appropriately with the dense backend (it is similar with the sparse one)
+More precisely, ProxQP API enables the user to provide for the init or update methods estimate of the minimal eigenvalue of H (i.e., manual_minimal_H_eigenvalue). It the values are not empty, then the values of primal proximal step size rho will be updated according to: rho = rho + abs(manual_minimal_H_eigenvalue). It guarantees that the proximal step-size is larger than the minimal eigenvalue of H and hence to converging towards a local minimum of the QP. We provide below examples in C++ and python for using this feature appropriately with the dense backend (it is similar with the sparse one)
 
 <table class="manual">
   <tr>

doc/3-ProxQP_solve.md

@@ -161,14 +161,4 @@ Note that if some elements of your QP model are not defined (for example a QP wi
   </tr>
 </table>
 
-Finally, note that you can also you ProxQP for solving QP with non convex quadratic. For doing so, you need first to estimate the smallest eigenvalue of the quadratic cost H. ProxQP has four internal options for estimating using the setting HessianCostRegularization:
-* NoRegularization : set by default, it means the solver does not try to estimate it,
-* Manual: the user can provide an estimate of it through the init method,
-* PowerIteration: a power iteration algorithm will be used for estimating H minimal eigenvalue,
-* EigenRegularization: in case the dense backend is used, the solver make use of Eigen method for estimating it.
-
-This option is particularly usefull when solving QP with non convex quadratics. Indeed, if default_rho is set to a value strictly higher than the minimal eigenvalue of H, then ProxQP is guaranteed for find a local minimum to the problem since it relies on a Proximal Method of Multipliers (for more detail for example this [work](https://arxiv.org/pdf/2010.02653.pdf) providing convergence proof of this property).
-
-More precisely, when HessianCostRegularization is set to a value different of NoRegularization, then ProxQP first estimate a minimal eigenvalue for H and then update default_rho following the rule: default_rho = rho_regularization_scaling * abs(default_H_eigenvalue_estimate), which guarantees for appropriate scaling than the proximal step-size is larger than the minimal eigenvalue of H. We provide below examples in C++ and python for using this feature appropriately with the dense backend (it is similar with the sparse one).
-
-The solve function enables using this option directly by passing accordingly the parameters "HessianCostRegularization" and "rho_regularization_scaling". You can find more details in [ProxQP API with examples](2-ProxQP_api.md) about the different other settings that can be used for setting other related parameters (e.g., for using PowerIteration algorithm with other options).
+Finally, note that you can also you ProxQP for solving QP with non convex quadratic. For doing so, you just need to provide to the solve function an estimate of the smallest eigenvalue of the quadratic cost H. The solver environment provides an independent function for estimating the minimal eigenvalue of a dense or sparse symmetric matrix. It is named "estimate_minimal_eigen_value_of_symmetric_matrix". You can find more details in [ProxQP API with examples](2-ProxQP_api.md) about the different other settings that can be used for setting other related parameters (e.g., for using a Power Iteration algorithm).