diff --git a/README.md b/README.md index d3d81e0d4..98205c5d0 100644 --- a/README.md +++ b/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. diff --git a/bindings/python/src/algorithms.hpp b/bindings/python/src/algorithms.hpp index 9549c2e31..766467259 100644 --- a/bindings/python/src/algorithms.hpp +++ b/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 diff --git a/bindings/python/src/expose-all.cpp b/bindings/python/src/expose-all.cpp index 95858bda4..64105733d 100644 --- a/bindings/python/src/expose-all.cpp +++ b/bindings/python/src/expose-all.cpp @@ -35,6 +35,7 @@ exposeSparseAlgorithms(pybind11::module_ m) sparse::python::exposeQpObjectSparse(m); sparse::python::exposeQPVectorSparse(m); sparse::python::solveSparseQp(m); + sparse::python::exposeSparseHelpers(m); } template @@ -45,6 +46,7 @@ exposeDenseAlgorithms(pybind11::module_ m) dense::python::exposeQpObjectDense(m); dense::python::exposeQPVectorDense(m); dense::python::solveDenseQp(m); + dense::python::exposeDenseHelpers(m); } #ifdef PROXSUITE_PYTHON_INTERFACE_WITH_OPENMP diff --git a/bindings/python/src/expose-helpers.hpp b/bindings/python/src/expose-helpers.hpp new file mode 100644 index 000000000..b9734b755 --- /dev/null +++ b/bindings/python/src/expose-helpers.hpp @@ -0,0 +1,72 @@ +// +// Copyright (c) 2022 INRIA +// + +#include +#include +#include + +#include +#include + +namespace proxsuite { +namespace proxqp { + +namespace dense { + +namespace python { + +template +void +exposeDenseHelpers(pybind11::module_ m) +{ + m.def("estimate_minimal_eigen_value_of_symmetric_matrix", + &dense::estimate_minimal_eigen_value_of_symmetric_matrix, + "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 +void +exposeSparseHelpers(pybind11::module_ m) +{ + m.def("estimate_minimal_eigen_value_of_symmetric_matrix", + &sparse::estimate_minimal_eigen_value_of_symmetric_matrix, + "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 diff --git a/bindings/python/src/expose-qpobject.hpp b/bindings/python/src/expose-qpobject.hpp index f07b8d945..31ca9513c 100644 --- a/bindings/python/src/expose-qpobject.hpp +++ b/bindings/python/src/expose-qpobject.hpp @@ -88,6 +88,7 @@ exposeQpObjectDense(pybind11::module_ m) bool compute_preconditioner, optional, optional, + optional, optional)>(&dense::QP::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::*)(optional>, optional>, @@ -119,6 +124,7 @@ exposeQpObjectDense(pybind11::module_ m) bool compute_preconditioner, optional, optional, + optional, optional)>(&dense::QP::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::*)()>(&dense::QP::solve), "function used for solving the QP problem, using default parameters.") @@ -163,6 +173,7 @@ exposeQpObjectDense(pybind11::module_ m) bool update_preconditioner, optional, optional, + optional, optional)>(&dense::QP::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::*)(optional>, @@ -198,6 +213,7 @@ exposeQpObjectDense(pybind11::module_ m) bool update_preconditioner, optional, optional, + optional, optional)>(&dense::QP::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::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::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::*)()>(&sparse::QP::solve), "function used for solving the QP problem, using default parameters.") diff --git a/bindings/python/src/expose-results.hpp b/bindings/python/src/expose-results.hpp index c431bebbb..45c2037b6 100644 --- a/bindings/python/src/expose-results.hpp +++ b/bindings/python/src/expose-results.hpp @@ -51,7 +51,12 @@ exposeResults(pybind11::module_ m) .def_readwrite("sparse_backend", &Info::sparse_backend, "Sparse backend used to solve the qp, either SparseCholesky " - "or MatrixFree."); + "or MatrixFree.") + .def_readwrite("minimal_H_eigenvalue_estimate", + &Info::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_>(m, "Results", pybind11::module_local()) .def(::pybind11::init(), diff --git a/bindings/python/src/expose-settings.hpp b/bindings/python/src/expose-settings.hpp index c481b1150..1f5949282 100644 --- a/bindings/python/src/expose-settings.hpp +++ b/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_( + m, "EigenValueEstimateMethodOption", pybind11::module_local()) + .value("PowerIteration", EigenValueEstimateMethodOption::PowerIteration) + .value("ExactMethod", EigenValueEstimateMethodOption::ExactMethod) + .export_values(); ::pybind11::class_>(m, "Settings", pybind11::module_local()) .def(::pybind11::init(), "Default constructor.") // constructor @@ -87,6 +92,8 @@ exposeSettings(pybind11::module_ m) &Settings::primal_infeasibility_solving) .def_readwrite("frequence_infeasibility_check", &Settings::frequence_infeasibility_check) + .def_readwrite("default_H_eigenvalue_estimate", + &Settings::default_H_eigenvalue_estimate) .def(pybind11::self == pybind11::self) .def(pybind11::self != pybind11::self) .def(pybind11::pickle( diff --git a/bindings/python/src/expose-solve.hpp b/bindings/python/src/expose-solve.hpp index 21a624d9b..8a0a2f091 100644 --- a/bindings/python/src/expose-solve.hpp +++ b/bindings/python/src/expose-solve.hpp @@ -43,7 +43,8 @@ solveDenseQp(pybind11::module_ m) bool, optional, optional, - bool>(&dense::solve), + bool, + optional>(&dense::solve), "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, optional, - bool>(&dense::solve), + bool, + optional>(&dense::solve), "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 diff --git a/doc/2-PROXQP_API/2-ProxQP_api.md b/doc/2-PROXQP_API/2-ProxQP_api.md index 26abc91a1..20d9302d3 100644 --- a/doc/2-PROXQP_API/2-ProxQP_api.md +++ b/doc/2-PROXQP_API/2-ProxQP_api.md @@ -73,7 +73,7 @@ $$\begin{equation}\label{eq:approx_dg_sol} where $[z]_+$ and $[z]_-$ stand for the z projection onto the positive and negative orthant. ProxQP provides the ``check_duality_gap`` option to include this duality gap in the stopping criterion. Note that it is disabled by default, as other solvers don't check this criterion in general. Enable this option if you want a stronger guarantee that your solution is optimal. ProxQP will then check the same termination condition as SCS (for more details see, e.g., SCS's [optimality conditions checks](https://www.cvxgrp.org/scs/algorithm/index.html#optimality-conditions) as well as [section 7.2](https://doi.org/10.1137/20M1366307) in the corresponding paper). The absolute and relative thresholds $\epsilon^{\text{gap}}_{\text{abs}}, \epsilon^{\text{gap}}_{\text{rel}}$ for the duality gap can differ from those $\epsilon_{\text{abs}}, \epsilon_{\text{rel}}$ for residuals because, contrary to residuals which result from an infinite norm, the duality gap scales with the square root of the problem dimension (thus it is numerically harder to achieve a given duality gap for larger problems). A recommended choice is $\epsilon^{\text{gap}}_{\text{abs}} = \epsilon_{\text{abs}} \sqrt{\max(n, n_{\text{eq}}, n_{\text{ineq}})}$. Note finally that meeting all residual and duality-gap criteria can be difficult for ill-conditioned problems. -Finally, note that ProxQP has a specific feature for handling primal infeasibility. More precisely, if the problem appears to be primal infeasible, it will solve the closest primal feasible problem in $$\ell_2$$ sense, and (x,y,z) will satisfy. +Finally, note that ProxQP has a specific feature for handling primal infeasibility. More precisely, if the problem appears to be primal infeasible, it will solve the closest primal feasible problem in $\ell_2$ sense, and (x,y,z) will satisfy. $$\begin{equation}\label{eq:approx_closest_qp_sol_rel} \begin{aligned} @@ -412,6 +412,12 @@ In this table, you have the three columns from left to right: the name of the se | preconditioner_accuracy | 1.E-3 | Accuracy level of the preconditioner. | HessianType | Dense | Defines the type of problem solved (Dense, Zero, or Diagonal). In case the Zero or Diagonal option is used, the solver exploits the Hessian structure to evaluate the Cholesky factorization efficiently. | primal_infeasibility_solving | False | If set to true, it solves the closest primal feasible problem if primal infeasibility is detected. +| 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. +| 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. \subsection OverviewInitialGuess The different initial guesses @@ -428,6 +434,35 @@ The different options will be commented below in the introduced order above. If set to this option, the solver will start with no initial guess, which means that the initial values of x, y, and z are the 0 vector. +\subsection OverviewEstimatingHminimalEigenValue The different options for estimating H minimal Eigenvalue + +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, +* ExactMethod: in this case, an exact method from EigenSolver is used to provide an estimate. + +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, 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) + + + + + + + + + + +
examples/cpp/estimate_nonconvex_eigenvalue.cppexamples/python/estimate_nonconvex_eigenvalue.py
+ \include estimate_nonconvex_eigenvalue.cpp + + \include estimate_nonconvex_eigenvalue.py +
+ +\subsubsection OverviewNoInitialGuess No initial guess + +If set to this option, the solver will start with no initial guess, which means that the initial values of x, y, and z are the 0 vector. + \subsubsection OverviewEqualityConstrainedInitialGuess Equality constrained initial guess If set to this option, the solver will solve at the beginning the following system for warm starting x and y. @@ -480,8 +515,8 @@ The result subclass is composed of the following: * x: a primal solution, * y: a Lagrange optimal multiplier for equality constraints, * z: a Lagrange optimal multiplier for inequality constraints, -* se: the optimal shift in $$\ell_2$$ with respect to equality constraints, -* si: the optimal shift in $$\ell_2$$ with respect to inequality constraints, +* se: the optimal shift in $\ell_2$ with respect to equality constraints, +* si: the optimal shift in $\ell_2$ with respect to inequality constraints, * info: a subclass which containts some information about the solver's execution. If the solver has solved the problem, the triplet (x,y,z) satisfies: @@ -500,7 +535,7 @@ $$\begin{equation}\label{eq:approx_qp_sol_rel} accordingly with the parameters eps_abs and eps_rel chosen by the user. If the problem is primal infeasible and you have enabled the solver to solve the closest feasible problem, then (x,y,z) will satisfy. -$$\begin{equation}\label{eq:approx_closest_qp_sol_rel} +$$\begin{equation}\label{eq:approx_closest_qp_sol_rel_bis} \begin{aligned} &\left\{ \begin{array}{ll} @@ -511,9 +546,9 @@ $$\begin{equation}\label{eq:approx_closest_qp_sol_rel} \end{aligned} \end{equation}$$ -(se, si) stands in this context for the optimal shifts in $$\ell_2$$ sense which enables recovering a primal feasible problem. More precisely, they are derived such that +(se, si) stands in this context for the optimal shifts in $\ell_2$ sense which enables recovering a primal feasible problem. More precisely, they are derived such that -$$\begin{equation}\label{eq:QP_primal_feasible}\tag{QP_feas} +\begin{equation}\label{eq:QP_primal_feasible}\tag{QP_feas} \begin{aligned} \min_{x\in\mathbb{R}^{d}} & \quad \frac{1}{2}x^{T}Hx+g^{T}x \\\ \text{s.t.}&\left\{ @@ -526,7 +561,7 @@ $$\begin{equation}\label{eq:QP_primal_feasible}\tag{QP_feas} \end{equation}\\\ defines a primal feasible problem. -Note that if you use the dense backend and its specific feature for handling box inequality constraints, then the first $$n_in$$ elements of z correspond to multipliers associated to the linear inequality formed with $$C$$ matrix, whereas the last $$d$$ elements correspond to multipliers associated to the box inequality constraints (see for example solve_dense_qp.cpp or solve_dense_qp.py). +Note that if you use the dense backend and its specific feature for handling box inequality constraints, then the first $n_{in}$ elements of z correspond to multipliers associated to the linear inequality formed with $C$ matrix, whereas the last $d$ elements correspond to multipliers associated to the box inequality constraints (see for example solve_dense_qp.cpp or solve_dense_qp.py). \subsection OverviewInfoClass The info subclass diff --git a/doc/3-ProxQP_solve.md b/doc/3-ProxQP_solve.md index 09410de8e..693b1038b 100644 --- a/doc/3-ProxQP_solve.md +++ b/doc/3-ProxQP_solve.md @@ -160,3 +160,5 @@ Note that if some elements of your QP model are not defined (for example a QP wi + +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). \ No newline at end of file diff --git a/examples/cpp/estimate_nonconvex_eigenvalue.cpp b/examples/cpp/estimate_nonconvex_eigenvalue.cpp new file mode 100644 index 000000000..25e51fdac --- /dev/null +++ b/examples/cpp/estimate_nonconvex_eigenvalue.cpp @@ -0,0 +1,52 @@ +#include // load the dense solver backend +#include // used for generating a random convex qp + +using namespace proxsuite; +using namespace proxsuite::proxqp; +using T = double; + +int +main() +{ + dense::isize dim = 10; + dense::isize n_eq(dim / 4); + dense::isize n_in(dim / 4); + // generate a random qp + T sparsity_factor(0.15); + T strong_convexity_factor(1.e-2); + dense::Model qp_random = utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + // make the QP nonconvex + dense::Vec diag(dim); + diag.setOnes(); + qp_random.H.diagonal().array() -= + 2. * diag.array(); // add some nonpositive values dense matrix + Eigen::SelfAdjointEigenSolver> es(qp_random.H, + Eigen::EigenvaluesOnly); + T minimal_eigenvalue = T(es.eigenvalues().minCoeff()); + // choose scaling for regularizing default_rho accordingly + dense::QP qp(dim, n_eq, n_in); // create the QP object + // choose the option for estimating this eigenvalue + T estimate_minimal_eigen_value = + dense::estimate_minimal_eigen_value_of_symmetric_matrix( + qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000); + bool compute_preconditioner = false; + // input the estimate for making rho appropriate + qp.init(qp_random.H, + qp_random.g, + qp_random.A, + qp_random.b, + qp_random.C, + qp_random.l, + qp_random.u, + compute_preconditioner, + nullopt, + nullopt, + nullopt, + estimate_minimal_eigen_value); + // print the estimates + std::cout << "ProxQP estimate " + << qp.results.info.minimal_H_eigenvalue_estimate << std::endl; + std::cout << "minimal_eigenvalue " << minimal_eigenvalue << std::endl; + std::cout << "default_rho " << qp.settings.default_rho << std::endl; +} diff --git a/examples/python/estimate_nonconvex_eigenvalue.py b/examples/python/estimate_nonconvex_eigenvalue.py new file mode 100644 index 000000000..21e2e54ef --- /dev/null +++ b/examples/python/estimate_nonconvex_eigenvalue.py @@ -0,0 +1,51 @@ +import proxsuite +import numpy as np +import scipy.sparse as spa + + +def generate_mixed_qp(n, seed=1, reg=-2.0): + # A function for generating sparse random convex qps in dense format + + np.random.seed(seed) + n_eq = int(n / 4) + n_in = int(n / 4) + m = n_eq + n_in + + P = spa.random( + n, n, density=0.075, data_rvs=np.random.randn, format="csc" + ).toarray() + P = (P + P.T) / 2.0 + + s = max(np.absolute(np.linalg.eigvals(P))) + P += (abs(s) + 1e-02) * spa.eye(n) + P = spa.coo_matrix(P) + q = np.random.randn(n) + A = spa.random(m, n, density=0.15, data_rvs=np.random.randn, format="csc").toarray() + v = np.random.randn(n) # Fictitious solution + delta = np.random.rand(m) # To get inequality + u = A @ v + l = -1.0e20 * np.ones(m) + + return P.toarray(), q, A[:n_eq, :], u[:n_eq], A[n_in:, :], u[n_in:], l[n_in:] + + +# load a qp object using qp problem dimensions +n = 10 +n_eq = 2 +n_in = 2 +qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) +# generate a random non convex QP +H, g, A, b, C, u, l = generate_mixed_qp(n) +# initialize the model of the problem to solve +estimate_minimal_eigen_value = ( + proxsuite.proxqp.dense.estimate_minimal_eigen_value_of_symmetric_matrix( + H, proxsuite.proxqp.EigenValueEstimateMethodOption.ExactMethod + ) +) +qp.init(H, g, A, b, C, l, u, manual_minimal_H_eigenvalue=estimate_minimal_eigen_value) +vals, _ = spa.linalg.eigs(H, which="SR") +min_eigenvalue = float(np.min(vals)) +# print the estimates +print(f"{min_eigenvalue=}") +print(f"{estimate_minimal_eigen_value=}") +print(f"{qp.results.info.minimal_H_eigenvalue_estimate=}") diff --git a/include/proxsuite/proxqp/dense/helpers.hpp b/include/proxsuite/proxqp/dense/helpers.hpp index e5a86c50e..cfef42082 100644 --- a/include/proxsuite/proxqp/dense/helpers.hpp +++ b/include/proxsuite/proxqp/dense/helpers.hpp @@ -15,12 +15,161 @@ #include #include #include +#include namespace proxsuite { namespace proxqp { namespace dense { +template +T +power_iteration(Mat& H, + Vec& dw, + Vec& rhs, + Vec& err_v, + T power_iteration_accuracy, + isize nb_power_iteration) +{ + // computes maximal eigen value of the bottom right matrix of the LDLT + isize dim = H.rows(); + rhs.setZero(); + // stores eigenvector + rhs.array() += 1. / std::sqrt(dim); + // stores Hx + dw.noalias() = H.template selfadjointView() * rhs; // Hx + T eig = 0; + for (isize i = 0; i < nb_power_iteration; i++) { + + rhs = dw / dw.norm(); + dw.noalias() = (H.template selfadjointView() * rhs); + // calculate associated eigenvalue + eig = rhs.dot(dw); + // calculate associated error + err_v = dw - eig * rhs; + T err = proxsuite::proxqp::dense::infty_norm(err_v); + // std::cout << "power iteration max: i " << i << " err " << err << + // std::endl; + if (err <= power_iteration_accuracy) { + break; + } + } + return eig; +} +template +T +min_eigen_value_via_modified_power_iteration(Mat& H, + Vec& dw, + Vec& rhs, + Vec& err_v, + T max_eigen_value, + T power_iteration_accuracy, + isize nb_power_iteration) +{ + // performs power iteration on the matrix: max_eigen_value I - H + // estimates then the minimal eigenvalue with: minimal_eigenvalue = + // max_eigen_value - eig + isize dim = H.rows(); + rhs.setZero(); + // stores eigenvector + rhs.array() += 1. / std::sqrt(dim); + // stores Hx + dw.noalias() = -(H.template selfadjointView() * rhs); // Hx + dw += max_eigen_value * rhs; + T eig = 0; + for (isize i = 0; i < nb_power_iteration; i++) { + + rhs = dw / dw.norm(); + dw.noalias() = -(H.template selfadjointView() * rhs); + dw += max_eigen_value * rhs; + // calculate associated eigenvalue + eig = rhs.dot(dw); + // calculate associated error + err_v = dw - eig * rhs; + T err = proxsuite::proxqp::dense::infty_norm(err_v); + // std::cout << "power iteration min: i " << i << " err " << err << + // std::endl; + if (err <= power_iteration_accuracy) { + break; + } + } + T minimal_eigenvalue = max_eigen_value - eig; + return minimal_eigenvalue; +} /////// SETUP //////// +/*! + * Estimate minimal eigenvalue of a symmetric Matrix + * @param H symmetric matrix. + * @param EigenValueEstimateMethodOption + * @param power_iteration_accuracy power iteration algorithm accuracy tracked + * @param nb_power_iteration maximal number of power iteration executed + * + */ +template +T +estimate_minimal_eigen_value_of_symmetric_matrix( + Mat& H, + EigenValueEstimateMethodOption estimate_method_option, + T power_iteration_accuracy, + isize nb_power_iteration) +{ + PROXSUITE_THROW_PRETTY((!H.isApprox(H.transpose(), 0.0)), + std::invalid_argument, + "H is not symmetric."); + if (H.size()) { + PROXSUITE_CHECK_ARGUMENT_SIZE( + H.rows(), + H.cols(), + "H has a number of rows different of the number of columns."); + } + isize dim = H.rows(); + T res(0.); + switch (estimate_method_option) { + case EigenValueEstimateMethodOption::PowerIteration: { + Vec dw(dim); + Vec rhs(dim); + Vec err(dim); + T dominant_eigen_value = power_iteration( + H, dw, rhs, err, power_iteration_accuracy, nb_power_iteration); + T min_eigenvalue = min_eigen_value_via_modified_power_iteration( + H, + dw, + rhs, + err, + dominant_eigen_value, + power_iteration_accuracy, + nb_power_iteration); + res = std::min(min_eigenvalue, dominant_eigen_value); + } break; + case EigenValueEstimateMethodOption::ExactMethod: { + Eigen::SelfAdjointEigenSolver> es(H, Eigen::EigenvaluesOnly); + res = T(es.eigenvalues()[0]); + } break; + } + return res; +} +/////// SETUP //////// +/*! + * Estimate H minimal eigenvalue + * @param settings solver settings + * @param results solver results. + * @param manual_minimal_H_eigenvalue minimal H eigenvalue estimate. + */ +template +void +update_default_rho_with_minimal_Hessian_eigen_value( + optional manual_minimal_H_eigenvalue, + Results& results, + Settings& settings) +{ + if (manual_minimal_H_eigenvalue != nullopt) { + settings.default_H_eigenvalue_estimate = + manual_minimal_H_eigenvalue.value(); + results.info.minimal_H_eigenvalue_estimate = + settings.default_H_eigenvalue_estimate; + } + settings.default_rho += std::abs(results.info.minimal_H_eigenvalue_estimate); + results.info.rho = settings.default_rho; +} /*! * Computes the equality constrained initial guess of a QP problem. * diff --git a/include/proxsuite/proxqp/dense/wrapper.hpp b/include/proxsuite/proxqp/dense/wrapper.hpp index 84db43526..22e58d82c 100644 --- a/include/proxsuite/proxqp/dense/wrapper.hpp +++ b/include/proxsuite/proxqp/dense/wrapper.hpp @@ -348,6 +348,7 @@ struct QP * @param rho proximal step size wrt primal variable. * @param mu_eq proximal step size wrt equality constrained multiplier. * @param mu_in proximal step size wrt inequality constrained multiplier. + * @param manual_minimal_H_eigenvalue manual minimal eigenvalue proposed for H */ void init(optional> H, optional> g, @@ -359,7 +360,8 @@ struct QP bool compute_preconditioner = true, optional rho = nullopt, optional mu_eq = nullopt, - optional mu_in = nullopt) + optional mu_in = nullopt, + optional manual_minimal_H_eigenvalue = nullopt) { PROXSUITE_THROW_PRETTY( box_constraints == true, @@ -467,6 +469,9 @@ struct QP } proxsuite::proxqp::dense::update_proximal_parameters( settings, results, work, rho, mu_eq, mu_in); + proxsuite::proxqp::dense:: + update_default_rho_with_minimal_Hessian_eigen_value( + manual_minimal_H_eigenvalue, results, settings); typedef optional> optional_VecRef; proxsuite::proxqp::dense::setup(H, g, @@ -509,6 +514,7 @@ struct QP * @param rho proximal step size wrt primal variable. * @param mu_eq proximal step size wrt equality constrained multiplier. * @param mu_in proximal step size wrt inequality constrained multiplier. + * @param manual_minimal_H_eigenvalue manual minimal eigenvalue proposed for H */ void init(optional> H, optional> g, @@ -522,7 +528,8 @@ struct QP bool compute_preconditioner = true, optional rho = nullopt, optional mu_eq = nullopt, - optional mu_in = nullopt) + optional mu_in = nullopt, + optional manual_minimal_H_eigenvalue = nullopt) { // dense case @@ -668,6 +675,9 @@ struct QP } proxsuite::proxqp::dense::update_proximal_parameters( settings, results, work, rho, mu_eq, mu_in); + proxsuite::proxqp::dense:: + update_default_rho_with_minimal_Hessian_eigen_value( + manual_minimal_H_eigenvalue, results, settings); proxsuite::proxqp::dense::setup(H, g, A, @@ -705,6 +715,7 @@ struct QP * @param rho proximal step size wrt primal variable. * @param mu_eq proximal step size wrt equality constrained multiplier. * @param mu_in proximal step size wrt inequality constrained multiplier. + * @param manual_minimal_H_eigenvalue manual minimal eigenvalue proposed for H * @note The init method should be called before update. If it has not been * done before, init is called depending on the is_initialized flag. */ @@ -718,7 +729,8 @@ struct QP bool update_preconditioner = false, optional rho = nullopt, optional mu_eq = nullopt, - optional mu_in = nullopt) + optional mu_in = nullopt, + optional manual_minimal_H_eigenvalue = nullopt) { PROXSUITE_THROW_PRETTY( box_constraints == true, @@ -764,7 +776,9 @@ struct QP } proxsuite::proxqp::dense::update_proximal_parameters( settings, results, work, rho, mu_eq, mu_in); - + proxsuite::proxqp::dense:: + update_default_rho_with_minimal_Hessian_eigen_value( + manual_minimal_H_eigenvalue, results, settings); typedef optional> optional_MatRef; typedef optional> optional_VecRef; proxsuite::proxqp::dense::setup(/* avoid double assignation */ @@ -809,6 +823,7 @@ struct QP * @param rho proximal step size wrt primal variable. * @param mu_eq proximal step size wrt equality constrained multiplier. * @param mu_in proximal step size wrt inequality constrained multiplier. + * @param manual_minimal_H_eigenvalue manual minimal eigenvalue proposed for H * @note The init method should be called before update. If it has not been * done before, init is called depending on the is_initialized flag. */ @@ -824,7 +839,8 @@ struct QP bool update_preconditioner = false, optional rho = nullopt, optional mu_eq = nullopt, - optional mu_in = nullopt) + optional mu_in = nullopt, + optional manual_minimal_H_eigenvalue = nullopt) { PROXSUITE_THROW_PRETTY( box_constraints == false && (l_box != nullopt || u_box != nullopt), @@ -871,6 +887,9 @@ struct QP } proxsuite::proxqp::dense::update_proximal_parameters( settings, results, work, rho, mu_eq, mu_in); + proxsuite::proxqp::dense:: + update_default_rho_with_minimal_Hessian_eigen_value( + manual_minimal_H_eigenvalue, results, settings); typedef optional> optional_MatRef; typedef optional> optional_VecRef; proxsuite::proxqp::dense::setup(/* avoid double assignation */ @@ -1004,7 +1023,8 @@ solve( bool check_duality_gap = false, optional eps_duality_gap_abs = nullopt, optional eps_duality_gap_rel = nullopt, - bool primal_infeasibility_solving = false) + bool primal_infeasibility_solving = false, + optional manual_minimal_H_eigenvalue = nullopt) { isize n(0); isize n_eq(0); @@ -1043,7 +1063,23 @@ solve( } Qp.settings.compute_timings = compute_timings; Qp.settings.primal_infeasibility_solving = primal_infeasibility_solving; - Qp.init(H, g, A, b, C, l, u, compute_preconditioner, rho, mu_eq, mu_in); + if (manual_minimal_H_eigenvalue != nullopt) { + Qp.init(H, + g, + A, + b, + C, + l, + u, + compute_preconditioner, + rho, + mu_eq, + mu_in, + manual_minimal_H_eigenvalue.value()); + } else { + Qp.init( + H, g, A, b, C, l, u, compute_preconditioner, rho, mu_eq, mu_in, nullopt); + } Qp.solve(x, y, z); return Qp.results; @@ -1119,7 +1155,8 @@ solve( bool check_duality_gap = false, optional eps_duality_gap_abs = nullopt, optional eps_duality_gap_rel = nullopt, - bool primal_infeasibility_solving = false) + bool primal_infeasibility_solving = false, + optional manual_minimal_H_eigenvalue = nullopt) { isize n(0); isize n_eq(0); @@ -1158,19 +1195,37 @@ solve( } Qp.settings.compute_timings = compute_timings; Qp.settings.primal_infeasibility_solving = primal_infeasibility_solving; - Qp.init(H, - g, - A, - b, - C, - l, - u, - l_box, - u_box, - compute_preconditioner, - rho, - mu_eq, - mu_in); + if (manual_minimal_H_eigenvalue != nullopt) { + Qp.init(H, + g, + A, + b, + C, + l, + u, + l_box, + u_box, + compute_preconditioner, + rho, + mu_eq, + mu_in, + manual_minimal_H_eigenvalue.value()); + } else { + Qp.init(H, + g, + A, + b, + C, + l, + u, + l_box, + u_box, + compute_preconditioner, + rho, + mu_eq, + mu_in, + nullopt); + } Qp.solve(x, y, z); return Qp.results; diff --git a/include/proxsuite/proxqp/results.hpp b/include/proxsuite/proxqp/results.hpp index dc97de46b..f3110dc13 100644 --- a/include/proxsuite/proxqp/results.hpp +++ b/include/proxsuite/proxqp/results.hpp @@ -52,6 +52,8 @@ struct Info T duality_gap; //// sparse backend used by solver, either CholeskySparse or MatrixFree SparseBackend sparse_backend; + //// quadratic cost minimal eigenvalue estimate + T minimal_H_eigenvalue_estimate; }; /// /// @brief This class stores all the results of PROXQP solvers with sparse and @@ -134,6 +136,7 @@ struct Results info.duality_gap = 0.; info.status = QPSolverOutput::PROXQP_NOT_RUN; info.sparse_backend = SparseBackend::Automatic; + info.minimal_H_eigenvalue_estimate = 0.; } /*! * cleanups the Result variables and set the info variables to their initial @@ -172,12 +175,15 @@ struct Results info.mu_in_inv = 1e1; info.mu_in = 1e-1; info.nu = 1.; + info.minimal_H_eigenvalue_estimate = 0.; if (settings != nullopt) { info.rho = settings.value().default_rho; info.mu_eq = settings.value().default_mu_eq; info.mu_eq_inv = T(1) / info.mu_eq; info.mu_in = settings.value().default_mu_in; info.mu_in_inv = T(1) / info.mu_in; + info.minimal_H_eigenvalue_estimate = + settings.value().default_H_eigenvalue_estimate; } cleanup_statistics(); } @@ -207,7 +213,8 @@ operator==(const Info& info1, const Info& info2) info1.solve_time == info2.solve_time && info1.run_time == info2.run_time && info1.objValue == info2.objValue && info1.pri_res == info2.pri_res && info1.dua_res == info2.dua_res && info1.duality_gap == info2.duality_gap && - info1.duality_gap == info2.duality_gap; + info1.duality_gap == info2.duality_gap && + info1.minimal_H_eigenvalue_estimate == info2.minimal_H_eigenvalue_estimate; return value; } diff --git a/include/proxsuite/proxqp/settings.hpp b/include/proxsuite/proxqp/settings.hpp index 3d8920b5b..b2e5188ec 100644 --- a/include/proxsuite/proxqp/settings.hpp +++ b/include/proxsuite/proxqp/settings.hpp @@ -43,6 +43,14 @@ enum struct HessianType Dense, // Quadratic Program Diagonal // Quadratic Program with diagonal Hessian }; +// Augment the minimal eigen value of H with some regularizer +enum struct EigenValueEstimateMethodOption +{ + PowerIteration, // Process a power iteration algorithm + ExactMethod // Use Eigen's method to estimate the minimal eigenvalue + // watch out, the last option is only available for dense + // matrices! +}; inline std::ostream& operator<<(std::ostream& os, const SparseBackend& sparse_backend) { @@ -132,6 +140,7 @@ struct Settings SparseBackend sparse_backend; bool primal_infeasibility_solving; isize frequence_infeasibility_check; + T default_H_eigenvalue_estimate; /*! * Default constructor. * @param default_rho default rho parameter of result class @@ -195,6 +204,10 @@ struct Settings * problem if activated * @param frequence_infeasibility_check frequence at which infeasibility is * checked + * @param find_H_minimal_eigenvalue track the minimal eigen value of the + * quadratic cost H + * @param default_H_eigenvalue_estimate default H eigenvalue estimate (i.e., + * if we make a model update and H does not change this one is used) */ Settings( @@ -243,7 +256,8 @@ struct Settings T alpha_gpdal = 0.95, SparseBackend sparse_backend = SparseBackend::Automatic, bool primal_infeasibility_solving = false, - isize frequence_infeasibility_check = 20) + isize frequence_infeasibility_check = 1, + T default_H_eigenvalue_estimate = 0.) : default_mu_eq(default_mu_eq) , default_mu_in(default_mu_in) , alpha_bcl(alpha_bcl) @@ -285,6 +299,7 @@ struct Settings , sparse_backend(sparse_backend) , primal_infeasibility_solving(primal_infeasibility_solving) , frequence_infeasibility_check(frequence_infeasibility_check) + , default_H_eigenvalue_estimate(default_H_eigenvalue_estimate) { switch (dense_backend) { case DenseBackend::PrimalDualLDLT: @@ -349,7 +364,9 @@ operator==(const Settings& settings1, const Settings& settings2) settings1.primal_infeasibility_solving == settings2.primal_infeasibility_solving && settings1.frequence_infeasibility_check == - settings2.frequence_infeasibility_check; + settings2.frequence_infeasibility_check && + settings1.default_H_eigenvalue_estimate == + settings2.default_H_eigenvalue_estimate; return value; } diff --git a/include/proxsuite/proxqp/sparse/helpers.hpp b/include/proxsuite/proxqp/sparse/helpers.hpp index 3bc25c21b..9674a6033 100644 --- a/include/proxsuite/proxqp/sparse/helpers.hpp +++ b/include/proxsuite/proxqp/sparse/helpers.hpp @@ -10,11 +10,165 @@ #include #include - +#include namespace proxsuite { namespace proxqp { namespace sparse { +template +T +power_iteration(SparseMat& H, + sparse::Vec& dw, + sparse::Vec& rhs, + sparse::Vec& err_v, + T power_iteration_accuracy, + isize nb_power_iteration) +{ + // computes maximal eigen value of the bottom right matrix of the LDLT + isize dim = H.rows(); + rhs.setZero(); + // stores eigenvector + rhs.array() += 1. / std::sqrt(dim); + // stores Hx + dw.setZero(); + // detail::noalias_symhiv_add( + // detail::vec_mut(dw), H, detail::vec_mut(rhs)); + dw = H * rhs; + T eig = 0; + for (isize i = 0; i < nb_power_iteration; i++) { + + rhs = dw / dw.norm(); + dw.setZero(); + // detail::noalias_symhiv_add( + // detail::vec_mut(dw), H, detail::vec_mut(rhs)); + dw = H * rhs; + // calculate associated eigenvalue + eig = rhs.dot(dw); + // calculate associated error + err_v = dw - eig * rhs; + T err = proxsuite::proxqp::dense::infty_norm(err_v); + // std::cout << "power iteration max: i " << i << " err " << err << + // std::endl; + if (err <= power_iteration_accuracy) { + break; + } + } + return eig; +} +template +T +min_eigen_value_via_modified_power_iteration(SparseMat& H, + sparse::Vec& dw, + sparse::Vec& rhs, + sparse::Vec& err_v, + T max_eigen_value, + T power_iteration_accuracy, + isize nb_power_iteration) +{ + // performs power iteration on the matrix: max_eigen_value I - H + // estimates then the minimal eigenvalue with: minimal_eigenvalue = + // max_eigen_value - eig + isize dim = H.rows(); + rhs.setZero(); + // stores eigenvector + rhs.array() += 1. / std::sqrt(dim); + // stores Hx + dw.setZero(); + // does not work below with full symmetry... + // detail::noalias_symhiv_add( + // detail::vec_mut(dw), H, detail::vec_mut(rhs)); + dw.noalias() = H * rhs; + dw.array() *= (-1.); + dw += max_eigen_value * rhs; + T eig = 0; + for (isize i = 0; i < nb_power_iteration; i++) { + + rhs = dw / dw.norm(); + dw.setZero(); + // idem + // detail::noalias_symhiv_add( + // detail::vec_mut(dw), H, detail::vec_mut(rhs)); + dw.noalias() = H * rhs; + dw.array() *= (-1.); + dw += max_eigen_value * rhs; + // calculate associated eigenvalue + eig = rhs.dot(dw); + // calculate associated error + err_v = dw - eig * rhs; + T err = proxsuite::proxqp::dense::infty_norm(err_v); + // std::cout << "power iteration min: i " << i << " err " << err << + // std::endl; + if (err <= power_iteration_accuracy) { + break; + } + } + T minimal_eigenvalue = max_eigen_value - eig; + return minimal_eigenvalue; +} +/////// SETUP //////// +/*! + * Estimate minimal eigenvalue of a symmetric Matrix via power iteration + * @param H symmetric matrix. + * @param power_iteration_accuracy power iteration algorithm accuracy tracked + * @param nb_power_iteration maximal number of power iteration executed + * + */ +template +T +estimate_minimal_eigen_value_of_symmetric_matrix(SparseMat& H, + T power_iteration_accuracy, + isize nb_power_iteration) +{ + PROXSUITE_THROW_PRETTY((!H.isApprox(H.transpose(), 0.0)), + std::invalid_argument, + "H is not symmetric."); + PROXSUITE_CHECK_ARGUMENT_SIZE( + H.rows(), + H.cols(), + "H has a number of rows different of the number of columns."); + isize dim = H.rows(); + T res(0.); + sparse::Vec dw(dim); + sparse::Vec rhs(dim); + sparse::Vec err(dim); + T dominant_eigen_value = power_iteration( + H, dw, rhs, err, power_iteration_accuracy, nb_power_iteration); + T min_eigenvalue = + min_eigen_value_via_modified_power_iteration(H, + dw, + rhs, + err, + dominant_eigen_value, + power_iteration_accuracy, + nb_power_iteration); + // std::cout << "dominant_eigen_value " << dominant_eigen_value<< " + // min_eigenvalue " << min_eigenvalue << std::endl; + res = std::min(min_eigenvalue, dominant_eigen_value); + return res; +} +/////// SETUP //////// +/*! + * Estimate H minimal eigenvalue + * @param settings solver settings + * @param results solver results. + * @param manual_minimal_H_eigenvalue minimal H eigenvalue estimate. + */ +template +void +update_default_rho_with_minimal_Hessian_eigen_value( + optional manual_minimal_H_eigenvalue, + Results& results, + Settings& settings) +{ + if (manual_minimal_H_eigenvalue != nullopt) { + settings.default_H_eigenvalue_estimate = + manual_minimal_H_eigenvalue.value(); + results.info.minimal_H_eigenvalue_estimate = + settings.default_H_eigenvalue_estimate; + } + settings.default_rho += std::abs(results.info.minimal_H_eigenvalue_estimate); + results.info.rho = settings.default_rho; +} /*! * Update the proximal parameters of the results object. * diff --git a/include/proxsuite/proxqp/sparse/wrapper.hpp b/include/proxsuite/proxqp/sparse/wrapper.hpp index 8070f0433..b9884bc7e 100644 --- a/include/proxsuite/proxqp/sparse/wrapper.hpp +++ b/include/proxsuite/proxqp/sparse/wrapper.hpp @@ -172,7 +172,8 @@ struct QP bool compute_preconditioner_ = true, optional rho = nullopt, optional mu_eq = nullopt, - optional mu_in = nullopt) + optional mu_in = nullopt, + optional manual_minimal_H_eigenvalue = nullopt) { if (settings.compute_timings) { work.timer.stop(); @@ -293,6 +294,7 @@ struct QP if (H != nullopt) { SparseMat H_triu = (H.value()).template triangularView(); + sparse::QpView qp = { { proxsuite::linalg::sparse::from_eigen, H_triu }, { proxsuite::linalg::sparse::from_eigen, model.g }, @@ -303,6 +305,9 @@ struct QP { proxsuite::linalg::sparse::from_eigen, model.u } }; qp_setup(qp, results, model, work, settings, ruiz, preconditioner_status); + proxsuite::proxqp::sparse:: + update_default_rho_with_minimal_Hessian_eigen_value( + manual_minimal_H_eigenvalue, results, settings); } else { SparseMat H_triu(model.dim, model.dim); H_triu.setZero(); @@ -354,7 +359,8 @@ struct QP bool update_preconditioner = false, optional rho = nullopt, optional mu_eq = nullopt, - optional mu_in = nullopt) + optional mu_in = nullopt, + optional manual_minimal_H_eigenvalue = nullopt) { if (!work.internal.is_initialized) { init(H, g, A, b, C, l, u, update_preconditioner, rho, mu_eq, mu_in); @@ -615,6 +621,11 @@ struct QP ruiz, preconditioner_status); // store model value + performs scaling // according to chosen options + if (H != nullopt) { + proxsuite::proxqp::sparse:: + update_default_rho_with_minimal_Hessian_eigen_value( + manual_minimal_H_eigenvalue, results, settings); + } if (settings.compute_timings) { results.info.setup_time = work.timer.elapsed().user; // in microseconds } @@ -720,7 +731,8 @@ solve( bool check_duality_gap = false, optional eps_duality_gap_abs = nullopt, optional eps_duality_gap_rel = nullopt, - bool primal_infeasibility_solving = false) + bool primal_infeasibility_solving = false, + optional manual_minimal_H_eigenvalue = nullopt) { isize n(0); @@ -761,7 +773,23 @@ solve( Qp.settings.compute_timings = compute_timings; Qp.settings.sparse_backend = sparse_backend; Qp.settings.primal_infeasibility_solving = primal_infeasibility_solving; - Qp.init(H, g, A, b, C, l, u, compute_preconditioner, rho, mu_eq, mu_in); + if (manual_minimal_H_eigenvalue != nullopt) { + Qp.init(H, + g, + A, + b, + C, + l, + u, + compute_preconditioner, + rho, + mu_eq, + mu_in, + manual_minimal_H_eigenvalue.value()); + } else { + Qp.init( + H, g, A, b, C, l, u, compute_preconditioner, rho, mu_eq, mu_in, nullopt); + } Qp.solve(x, y, z); return Qp.results; diff --git a/test/src/dense_qp_wrapper.cpp b/test/src/dense_qp_wrapper.cpp index 70f4e68c8..27a926fc7 100644 --- a/test/src/dense_qp_wrapper.cpp +++ b/test/src/dense_qp_wrapper.cpp @@ -7206,4 +7206,361 @@ TEST_CASE("ProxQP::dense: test primal infeasibility solving") DOCTEST_CHECK(pri_res <= scaled_eps); DOCTEST_CHECK(dua_res <= eps_abs); } -} \ No newline at end of file +} + +TEST_CASE("ProxQP::dense: estimate of minimal eigenvalues using Eigen") +{ + double sparsity_factor = 1.; + T tol = T(1e-6); + utils::rand::set_seed(1); + dense::isize dim = 2; + dense::isize n_eq(dim); + dense::isize n_in(dim); + T strong_convexity_factor(1.e-2); + for (isize i = 0; i < 1; ++i) { + // trivial test + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + qp_random.H.setZero(); + qp_random.H.diagonal().setOnes(); + qp_random.H.diagonal().tail(1).setConstant(-1.); + + T estimate_minimal_eigen_value = + dense::estimate_minimal_eigen_value_of_symmetric_matrix( + qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000); + + proxqp::dense::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + qp.init(qp_random.H, + qp_random.g, + qp_random.A, + qp_random.b, + qp_random.C, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + estimate_minimal_eigen_value); + + DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate + 1) <= + tol); + } + dim = 50; + n_eq = dim; + n_in = dim; + for (isize i = 0; i < 20; ++i) { + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + qp_random.H.setZero(); + dense::Vec random_diag = proxqp::utils::rand::vector_rand(dim); + qp_random.H.diagonal().array() += random_diag.array(); + T minimal_eigenvalue = qp_random.H.diagonal().minCoeff(); + + T estimate_minimal_eigen_value = + dense::estimate_minimal_eigen_value_of_symmetric_matrix( + qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000); + + proxqp::dense::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + qp.init(qp_random.H, + qp_random.g, + qp_random.A, + qp_random.b, + qp_random.C, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + estimate_minimal_eigen_value); + DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate - + minimal_eigenvalue) <= tol); + } + dim = 50; + n_eq = dim; + n_in = dim; + for (isize i = 0; i < 20; ++i) { + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + dense::Vec random_diag = proxqp::utils::rand::vector_rand(dim); + qp_random.H.diagonal().array() += 100 * random_diag.array(); + Eigen::SelfAdjointEigenSolver> es(qp_random.H, + Eigen::EigenvaluesOnly); + T minimal_eigenvalue = T(es.eigenvalues().minCoeff()); + + T estimate_minimal_eigen_value = + dense::estimate_minimal_eigen_value_of_symmetric_matrix( + qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000); + + proxqp::dense::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + qp.init(qp_random.H, + qp_random.g, + qp_random.A, + qp_random.b, + qp_random.C, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + estimate_minimal_eigen_value); + + DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate - + minimal_eigenvalue) <= tol); + } +} + +TEST_CASE( + "ProxQP::dense: test estimate of minimal eigenvalue using manual choice") +{ + double sparsity_factor = 1.; + T tol = T(1e-6); + utils::rand::set_seed(1); + dense::isize dim = 2; + dense::isize n_eq(dim); + dense::isize n_in(dim); + T strong_convexity_factor(1.e-2); + for (isize i = 0; i < 1; ++i) { + // trivial test + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + qp_random.H.setZero(); + qp_random.H.diagonal().setOnes(); + qp_random.H.diagonal().tail(1).setConstant(-1.); + + proxqp::dense::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + qp.init(qp_random.H, + qp_random.g, + qp_random.A, + qp_random.b, + qp_random.C, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + -1); + + DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate + 1) <= + tol); + } + dim = 50; + n_eq = dim; + n_in = dim; + for (isize i = 0; i < 20; ++i) { + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + qp_random.H.setZero(); + dense::Vec random_diag = proxqp::utils::rand::vector_rand(dim); + qp_random.H.diagonal().array() += random_diag.array(); + T minimal_eigenvalue = qp_random.H.diagonal().minCoeff(); + + proxqp::dense::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + qp.init(qp_random.H, + qp_random.g, + qp_random.A, + qp_random.b, + qp_random.C, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + minimal_eigenvalue); + DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate - + minimal_eigenvalue) <= tol); + } + dim = 50; + n_eq = dim; + n_in = dim; + for (isize i = 0; i < 20; ++i) { + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + dense::Vec random_diag = proxqp::utils::rand::vector_rand(dim); + qp_random.H.diagonal().array() += 100 * random_diag.array(); + Eigen::SelfAdjointEigenSolver> es(qp_random.H, + Eigen::EigenvaluesOnly); + T minimal_eigenvalue = T(es.eigenvalues().minCoeff()); + + proxqp::dense::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + qp.init(qp_random.H, + qp_random.g, + qp_random.A, + qp_random.b, + qp_random.C, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + minimal_eigenvalue); + + DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate - + minimal_eigenvalue) <= tol); + } +} + +TEST_CASE( + "ProxQP::dense: test estimate of minimal eigenvalue using power iteration") +{ + double sparsity_factor = 1.; + T tol = T(1e-3); + utils::rand::set_seed(1); + dense::isize dim = 2; + dense::isize n_eq(dim); + dense::isize n_in(dim); + T strong_convexity_factor(1.e-2); + for (isize i = 0; i < 1; ++i) { + // trivial test + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + qp_random.H.setZero(); + qp_random.H.diagonal().setOnes(); + qp_random.H.diagonal().tail(1).setConstant(-0.5); + + T estimate_minimal_eigen_value = + dense::estimate_minimal_eigen_value_of_symmetric_matrix( + qp_random.H, + EigenValueEstimateMethodOption::PowerIteration, + 1.E-6, + 10000); + + proxqp::dense::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + qp.init(qp_random.H, + qp_random.g, + qp_random.A, + qp_random.b, + qp_random.C, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + estimate_minimal_eigen_value); + + DOCTEST_CHECK( + std::abs(qp.results.info.minimal_H_eigenvalue_estimate + 0.5) <= tol); + } + dim = 50; + n_eq = dim; + n_in = dim; + for (isize i = 0; i < 20; ++i) { + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + qp_random.H.setZero(); + dense::Vec random_diag = proxqp::utils::rand::vector_rand(dim); + qp_random.H.diagonal().array() += random_diag.array(); + T minimal_eigenvalue = qp_random.H.diagonal().minCoeff(); + + T estimate_minimal_eigen_value = + dense::estimate_minimal_eigen_value_of_symmetric_matrix( + qp_random.H, + EigenValueEstimateMethodOption::PowerIteration, + 1.E-6, + 10000); + + proxqp::dense::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + qp.init(qp_random.H, + qp_random.g, + qp_random.A, + qp_random.b, + qp_random.C, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + estimate_minimal_eigen_value); + DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate - + minimal_eigenvalue) <= tol); + } + dim = 50; + n_eq = dim; + n_in = dim; + for (isize i = 0; i < 20; ++i) { + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + dense::Vec random_diag = proxqp::utils::rand::vector_rand(dim); + qp_random.H.diagonal().array() += + 100 * random_diag.array(); // add some random values to dense matrix + Eigen::SelfAdjointEigenSolver> es(qp_random.H, + Eigen::EigenvaluesOnly); + T minimal_eigenvalue = T(es.eigenvalues().minCoeff()); + + T estimate_minimal_eigen_value = + dense::estimate_minimal_eigen_value_of_symmetric_matrix( + qp_random.H, + EigenValueEstimateMethodOption::PowerIteration, + 1.E-6, + 10000); + + proxqp::dense::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + qp.init(qp_random.H, + qp_random.g, + qp_random.A, + qp_random.b, + qp_random.C, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + estimate_minimal_eigen_value); + + DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate - + minimal_eigenvalue) <= tol); + } +} diff --git a/test/src/dense_qp_wrapper.py b/test/src/dense_qp_wrapper.py index 768791fb9..7e2d63119 100644 --- a/test/src/dense_qp_wrapper.py +++ b/test/src/dense_qp_wrapper.py @@ -14,7 +14,7 @@ def normInf(x): return np.linalg.norm(x, np.inf) -def generate_mixed_qp(n, seed=1): +def generate_mixed_qp(n, seed=1, reg=0.01): """ Generate sparse problem in dense QP format """ @@ -31,7 +31,7 @@ def generate_mixed_qp(n, seed=1): P = (P + P.T) / 2.0 s = max(np.absolute(np.linalg.eigvals(P))) - P += (abs(s) + 1e-02) * spa.eye(n) + P += (abs(s) + reg) * spa.eye(n) P = spa.coo_matrix(P) print("sparsity of P : {}".format((P.nnz) / (n**2))) q = np.random.randn(n) @@ -4742,6 +4742,119 @@ def test_dense_infeasibility_solving( assert dua_res <= qp.settings.eps_abs assert pri_res <= scaled_eps + def test_minimal_eigenvalue_estimation_nonconvex_eigen_option( + self, + ): + print( + "------------------------dense non convex qp with inequality constraints, estimate minimal eigenvalue with eigen method" + ) + n = 50 + tol = 1.0e-3 + for i in range(50): + H, g, A, b, C, u, l = generate_mixed_qp(n, i, -0.01) + n_eq = A.shape[0] + n_in = C.shape[0] + qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) + qp.settings.verbose = False + qp.settings.initial_guess = proxsuite.proxqp.InitialGuess.NO_INITIAL_GUESS + estimate_minimal_eigen_value = ( + proxsuite.proxqp.dense.estimate_minimal_eigen_value_of_symmetric_matrix( + H, + proxsuite.proxqp.EigenValueEstimateMethodOption.ExactMethod, + 1.0e-6, + 10000, + ) + ) + vals, _ = spa.linalg.eigs(H, which="SR") + min_eigenvalue = float(np.min(vals)) + qp.init( + H, + np.asfortranarray(g), + A, + np.asfortranarray(b), + C, + np.asfortranarray(l), + np.asfortranarray(u), + manual_minimal_H_eigenvalue=estimate_minimal_eigen_value, + ) + assert ( + np.abs(min_eigenvalue - qp.results.info.minimal_H_eigenvalue_estimate) + <= tol + ) + + def test_minimal_eigenvalue_estimation_nonconvex_manual_option( + self, + ): + print( + "------------------------dense non convex qp with inequality constraints, estimate minimal eigenvalue with manual option" + ) + n = 50 + tol = 1.0e-3 + for i in range(50): + H, g, A, b, C, u, l = generate_mixed_qp(n, i, -0.01) + n_eq = A.shape[0] + n_in = C.shape[0] + qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) + qp.settings.verbose = False + qp.settings.initial_guess = proxsuite.proxqp.InitialGuess.NO_INITIAL_GUESS + vals, _ = spa.linalg.eigs(H, which="SR") + min_eigenvalue = float(np.min(vals)) + qp.init( + H, + np.asfortranarray(g), + A, + np.asfortranarray(b), + C, + np.asfortranarray(l), + np.asfortranarray(u), + manual_minimal_H_eigenvalue=min_eigenvalue, + ) + assert ( + np.abs(min_eigenvalue - qp.results.info.minimal_H_eigenvalue_estimate) + <= tol + ) + + def test_minimal_eigenvalue_estimation_nonconvex_power_iter_option( + self, + ): + print( + "------------------------dense non convex qp with inequality constraints, estimate minimal eigenvalue with power iter option" + ) + n = 50 + tol = 1.0e-3 + for i in range(50): + H, g, A, b, C, u, l = generate_mixed_qp(n, i, -0.01) + n_eq = A.shape[0] + n_in = C.shape[0] + + qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) + qp.settings.verbose = False + qp.settings.initial_guess = proxsuite.proxqp.InitialGuess.NO_INITIAL_GUESS + estimate_minimal_eigen_value = ( + proxsuite.proxqp.dense.estimate_minimal_eigen_value_of_symmetric_matrix( + H, + proxsuite.proxqp.EigenValueEstimateMethodOption.PowerIteration, + 1.0e-6, + 10000, + ) + ) + vals, _ = spa.linalg.eigs(H, which="SR") + min_eigenvalue = float(np.min(vals)) + qp.init( + H, + np.asfortranarray(g), + A, + np.asfortranarray(b), + C, + np.asfortranarray(l), + np.asfortranarray(u), + manual_minimal_H_eigenvalue=estimate_minimal_eigen_value, + ) + assert ( + np.abs(min_eigenvalue - qp.results.info.minimal_H_eigenvalue_estimate) + <= tol + ) + if __name__ == "__main__": unittest.main() diff --git a/test/src/sparse_qp_wrapper.cpp b/test/src/sparse_qp_wrapper.cpp index 85648a5e2..f8d07398a 100644 --- a/test/src/sparse_qp_wrapper.cpp +++ b/test/src/sparse_qp_wrapper.cpp @@ -1,5 +1,5 @@ // -// Copyright (c) 2022 INRIA +// Copyright (c) 2022-2023 INRIA // #include #include @@ -13,6 +13,7 @@ using namespace proxsuite::proxqp::utils; using T = double; using I = c_int; using namespace proxsuite::linalg::sparse::tags; + DOCTEST_TEST_CASE( "ProxQP::sparse: sparse random strongly convex qp with inequality constraints" "and empty equality constraints") @@ -6049,6 +6050,7 @@ TEST_CASE("ProxQP::sparse: init must be called before update") DOCTEST_CHECK(pri_res <= eps_abs); DOCTEST_CHECK(dua_res <= eps_abs); } + TEST_CASE("ProxQP::sparse: test primal infeasibility solving") { double sparsity_factor = 0.15; @@ -6108,4 +6110,262 @@ TEST_CASE("ProxQP::sparse: test primal infeasibility solving") DOCTEST_CHECK(pri_res <= scaled_eps); DOCTEST_CHECK(dua_res <= eps_abs); } -} \ No newline at end of file +} +// TEST_CASE("ProxQP::sparse: estimate of minimal eigenvalues using Eigen") +// { +// double sparsity_factor = 0.25; +// T tol = T(1e-6); +// utils::rand::set_seed(1); +// dense::isize dim = 2; +// dense::isize n_eq(dim); +// dense::isize n_in(dim); +// T strong_convexity_factor(1.e-2); +// dim = 50; +// n_eq = dim; +// n_in = dim; +// for (isize i = 0; i < 20; ++i) { +// ::proxsuite::proxqp::utils::rand::set_seed(i); +// proxqp::dense::Model qp_random = +// proxqp::utils::dense_strongly_convex_qp( +// dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); +// // proxqp::sparse::SparseModel qp_random = +// // utils::sparse_strongly_convex_qp( +// // dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + +// qp_random.H.setZero(); +// dense::Vec random_diag = proxqp::utils::rand::vector_rand(dim); +// qp_random.H.diagonal().array() += random_diag.array(); +// T minimal_eigenvalue = qp_random.H.diagonal().minCoeff(); + +// proxqp::sparse::QP qp(dim, n_eq, n_in); +// qp.settings.max_iter = 1; +// qp.settings.max_iter_in = 1; +// qp.settings.estimate_method_option = +// EigenValueEstimateMethodOption::EigenRegularization; +// qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; +// SparseMat H_sparse = qp_random.H.sparseView(); +// SparseMat A_sparse = qp_random.A.sparseView(); +// SparseMat C_sparse = qp_random.C.sparseView(); +// qp.init(H_sparse, +// qp_random.g, +// A_sparse, +// qp_random.b, +// C_sparse, +// qp_random.l, +// qp_random.u); +// DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate - +// minimal_eigenvalue) <= tol); +// } +// dim = 50; +// n_eq = dim; +// n_in = dim; +// for (isize i = 0; i < 20; ++i) { +// ::proxsuite::proxqp::utils::rand::set_seed(i); +// proxqp::dense::Model qp_random = +// proxqp::utils::dense_strongly_convex_qp( +// dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + +// dense::Vec random_diag = proxqp::utils::rand::vector_rand(dim); +// qp_random.H.diagonal().array() += 100 * random_diag.array(); + +// proxqp::sparse::QP qp(dim, n_eq, n_in); +// qp.settings.max_iter = 1; +// qp.settings.max_iter_in = 1; +// qp.settings.estimate_method_option = +// EigenValueEstimateMethodOption::EigenRegularization; +// qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; +// SparseMat H_sparse = qp_random.H.sparseView(); +// SparseMat A_sparse = qp_random.A.sparseView(); +// SparseMat C_sparse = qp_random.C.sparseView(); +// Eigen::SelfAdjointEigenSolver> es(H_sparse, +// Eigen::EigenvaluesOnly); +// T minimal_eigenvalue = T(es.eigenvalues().minCoeff()); +// qp.init(H_sparse, +// qp_random.g, +// A_sparse, +// qp_random.b, +// C_sparse, +// qp_random.l, +// qp_random.u); +// DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate - +// minimal_eigenvalue) <= 1.); +// } +// } +TEST_CASE("ProxQP::sparse: estimate of minimal eigenvalues using manual choice") +{ + double sparsity_factor = 0.25; + T tol = T(1e-6); + utils::rand::set_seed(1); + dense::isize dim = 2; + dense::isize n_eq(dim); + dense::isize n_in(dim); + T strong_convexity_factor(1.e-2); + dim = 50; + n_eq = dim; + n_in = dim; + for (isize i = 0; i < 20; ++i) { + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + // proxqp::sparse::SparseModel qp_random = + // utils::sparse_strongly_convex_qp( + // dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + qp_random.H.setZero(); + dense::Vec random_diag = proxqp::utils::rand::vector_rand(dim); + qp_random.H.diagonal().array() += random_diag.array(); + T minimal_eigenvalue = qp_random.H.diagonal().minCoeff(); + + proxqp::sparse::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + SparseMat H_sparse = qp_random.H.sparseView(); + SparseMat A_sparse = qp_random.A.sparseView(); + SparseMat C_sparse = qp_random.C.sparseView(); + qp.init(H_sparse, + qp_random.g, + A_sparse, + qp_random.b, + C_sparse, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + minimal_eigenvalue); + DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate - + minimal_eigenvalue) <= tol); + } + dim = 50; + n_eq = dim; + n_in = dim; + for (isize i = 0; i < 20; ++i) { + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + dense::Vec random_diag = proxqp::utils::rand::vector_rand(dim); + qp_random.H.diagonal().array() += 100 * random_diag.array(); + + proxqp::sparse::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + SparseMat H_sparse = qp_random.H.sparseView(); + SparseMat A_sparse = qp_random.A.sparseView(); + SparseMat C_sparse = qp_random.C.sparseView(); + Eigen::SelfAdjointEigenSolver> es(qp_random.H, + Eigen::EigenvaluesOnly); + T minimal_eigenvalue = T(es.eigenvalues().minCoeff()); + qp.init(H_sparse, + qp_random.g, + A_sparse, + qp_random.b, + C_sparse, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + minimal_eigenvalue); + DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate - + minimal_eigenvalue) <= 1.); + } +} + +TEST_CASE( + "ProxQP::sparse: estimate of minimal eigenvalues using Power iteration") +{ + double sparsity_factor = 0.25; + T tol = T(1e-6); + utils::rand::set_seed(1); + dense::isize dim = 2; + dense::isize n_eq(dim); + dense::isize n_in(dim); + T strong_convexity_factor(1.e-2); + dim = 50; + n_eq = dim; + n_in = dim; + for (isize i = 0; i < 20; ++i) { + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + // proxqp::sparse::SparseModel qp_random = + // utils::sparse_strongly_convex_qp( + // dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + qp_random.H.setZero(); + dense::Vec random_diag = proxqp::utils::rand::vector_rand(dim); + qp_random.H.diagonal().array() += random_diag.array(); + T minimal_eigenvalue = qp_random.H.diagonal().minCoeff(); + // std::cout << "qp_random.H" << std::endl; + // std::cout << qp_random.H << std::endl; + // std::cout << "minimal_eigenvalue " << minimal_eigenvalue << std::endl; + proxqp::sparse::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + SparseMat H_sparse = qp_random.H.sparseView(); + SparseMat A_sparse = qp_random.A.sparseView(); + SparseMat C_sparse = qp_random.C.sparseView(); + const T estimate_minimal_eigen_value = + sparse::estimate_minimal_eigen_value_of_symmetric_matrix( + H_sparse, 1.E-6, 10000); + qp.init(H_sparse, + qp_random.g, + A_sparse, + qp_random.b, + C_sparse, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + estimate_minimal_eigen_value); + DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate - + minimal_eigenvalue) <= tol); + } + dim = 50; + n_eq = dim; + n_in = dim; + for (isize i = 0; i < 20; ++i) { + ::proxsuite::proxqp::utils::rand::set_seed(i); + proxqp::dense::Model qp_random = proxqp::utils::dense_strongly_convex_qp( + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); + + dense::Vec random_diag = proxqp::utils::rand::vector_rand(dim); + qp_random.H.diagonal().array() += 100 * random_diag.array(); + + proxqp::sparse::QP qp(dim, n_eq, n_in); + qp.settings.max_iter = 1; + qp.settings.max_iter_in = 1; + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; + SparseMat H_sparse = qp_random.H.sparseView(); + SparseMat A_sparse = qp_random.A.sparseView(); + SparseMat C_sparse = qp_random.C.sparseView(); + T estimate_minimal_eigen_value = + sparse::estimate_minimal_eigen_value_of_symmetric_matrix( + H_sparse, 1.E-6, 10000); + Eigen::SelfAdjointEigenSolver> es(qp_random.H, + Eigen::EigenvaluesOnly); + const T minimal_eigenvalue = T(es.eigenvalues().minCoeff()); + qp.init(H_sparse, + qp_random.g, + A_sparse, + qp_random.b, + C_sparse, + qp_random.l, + qp_random.u, + true, + nullopt, + nullopt, + nullopt, + estimate_minimal_eigen_value); + DOCTEST_CHECK(std::abs(qp.results.info.minimal_H_eigenvalue_estimate - + minimal_eigenvalue) <= 1.); + } +} diff --git a/test/src/sparse_qp_wrapper.py b/test/src/sparse_qp_wrapper.py index 41fd9697e..466ca2f2d 100644 --- a/test/src/sparse_qp_wrapper.py +++ b/test/src/sparse_qp_wrapper.py @@ -14,7 +14,7 @@ def normInf(x): return np.linalg.norm(x, np.inf) -def generate_mixed_qp(n, seed=1): +def generate_mixed_qp(n, seed=1, reg=0.01): """ Generate problem in QP format """ @@ -31,7 +31,7 @@ def generate_mixed_qp(n, seed=1): P = (P + P.T) / 2.0 s = max(np.absolute(np.linalg.eigvals(P))) - P += (abs(s) + 1e-02) * spa.eye(n) + P += (abs(s) + reg) * spa.eye(n) P = spa.coo_matrix(P) print("sparsity of P : {}".format((P.nnz) / (n**2))) q = np.random.randn(n) @@ -46,7 +46,6 @@ def generate_mixed_qp(n, seed=1): class SparseqpWrapper(unittest.TestCase): # TESTS OF GENERAL METHODS OF THE API - """ def test_case_update_rho(self): print( "------------------------sparse random strongly convex qp with equality and inequality constraints: test update rho" @@ -4575,7 +4574,6 @@ def test_initializing_with_None(self): qp.results.info.setup_time, qp.results.info.solve_time ) ) - """ def test_sparse_infeasibility_solving( self, @@ -4626,6 +4624,116 @@ def test_sparse_infeasibility_solving( assert dua_res <= qp.settings.eps_abs assert pri_res <= scaled_eps + # def test_minimal_eigenvalue_estimation_nonconvex_eigen_option( + # self, + # ): + # print( + # "------------------------dense non convex qp with inequality constraints, estimate minimal eigenvalue with eigen method" + # ) + # n = 50 + # tol = 1.0 + # for i in range(50): + # H, g, A, b, C, u, l = generate_mixed_qp(n, i,-0.01) + # n_eq = A.shape[0] + # n_in = C.shape[0] + # qp = proxsuite.proxqp.sparse.QP(n, n_eq, n_in) + # qp.settings.verbose = False + # qp.settings.initial_guess = proxsuite.proxqp.InitialGuess.NO_INITIAL_GUESS + # qp.settings.estimate_method_option = ( + # proxsuite.proxqp.EigenValueEstimateMethodOption.EigenRegularization + # ) + # vals, _ = spa.linalg.eigs(H, which="SR") + # min_eigenvalue = float(np.min(vals)) + # qp.init( + # H, + # np.asfortranarray(g), + # A, + # np.asfortranarray(b), + # C, + # np.asfortranarray(l), + # np.asfortranarray(u), + # ) + # print(f"{min_eigenvalue=}") + # print(f"{qp.results.info.minimal_H_eigenvalue_estimate=}") + # input() + # assert ( + # np.abs(min_eigenvalue - qp.results.info.minimal_H_eigenvalue_estimate) + # <= tol + # ) + + def test_minimal_eigenvalue_estimation_nonconvex_manual_option( + self, + ): + print( + "------------------------dense non convex qp with inequality constraints, estimate minimal eigenvalue with manual option" + ) + n = 50 + tol = 1.0e-3 + for i in range(50): + H, g, A, b, C, u, l = generate_mixed_qp(n, i, -0.01) + n_eq = A.shape[0] + n_in = C.shape[0] + qp = proxsuite.proxqp.sparse.QP(n, n_eq, n_in) + qp.settings.verbose = False + qp.settings.initial_guess = proxsuite.proxqp.InitialGuess.NO_INITIAL_GUESS + vals, _ = spa.linalg.eigs(H, which="SR") + min_eigenvalue = float(np.min(vals)) + qp.init( + H, + np.asfortranarray(g), + A, + np.asfortranarray(b), + C, + np.asfortranarray(l), + np.asfortranarray(u), + manual_minimal_H_eigenvalue=min_eigenvalue, + ) + assert ( + np.abs(min_eigenvalue - qp.results.info.minimal_H_eigenvalue_estimate) + <= tol + ) + + def test_minimal_eigenvalue_estimation_nonconvex_power_iter_option( + self, + ): + print( + "------------------------sparse non convex qp with inequality constraints, estimate minimal eigenvalue with power iter option" + ) + n = 50 + tol = 1.0 + for i in range(50): + H, g, A, b, C, u, l = generate_mixed_qp(n, i, -0.01) + n_eq = A.shape[0] + n_in = C.shape[0] + qp = proxsuite.proxqp.sparse.QP(n, n_eq, n_in) + qp.settings.verbose = False + qp.settings.initial_guess = proxsuite.proxqp.InitialGuess.NO_INITIAL_GUESS + estimate_minimal_eigen_value = proxsuite.proxqp.sparse.estimate_minimal_eigen_value_of_symmetric_matrix( + H, 1.0e-10, 100000 + ) + vals, _ = spa.linalg.eigs(H, which="SR") + min_eigenvalue = float(np.min(vals)) + qp.init( + H, + np.asfortranarray(g), + A, + np.asfortranarray(b), + C, + np.asfortranarray(l), + np.asfortranarray(u), + manual_minimal_H_eigenvalue=estimate_minimal_eigen_value, + ) + # vals_bis, _ = spa.linalg.eigs(H, which="LM") + # print(f"{vals_bis}=") + # print(f"{vals}=") + # print(f"{min_eigenvalue=}") + # print(f"{qp.results.info.minimal_H_eigenvalue_estimate=}") + # input() + assert ( + np.abs(min_eigenvalue - qp.results.info.minimal_H_eigenvalue_estimate) + <= tol + ) + if __name__ == "__main__": unittest.main()