diff --git a/bindings/python/src/expose-helpers.hpp b/bindings/python/src/expose-helpers.hpp index b9734b755..6dc8768e5 100644 --- a/bindings/python/src/expose-helpers.hpp +++ b/bindings/python/src/expose-helpers.hpp @@ -20,25 +20,35 @@ 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.")); + m.def( + "estimate_minimal_eigen_value_of_symmetric_matrix", + +[](const MatRef& H, + EigenValueEstimateMethodOption estimate_method_option, + T power_iteration_accuracy, + isize nb_power_iteration) { + return dense::estimate_minimal_eigen_value_of_symmetric_matrix( + H, + estimate_method_option, + power_iteration_accuracy, + nb_power_iteration); + }, + "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 diff --git a/examples/cpp/estimate_nonconvex_eigenvalue.cpp b/examples/cpp/estimate_nonconvex_eigenvalue.cpp index 2c46ec4a2..658077bc1 100644 --- a/examples/cpp/estimate_nonconvex_eigenvalue.cpp +++ b/examples/cpp/estimate_nonconvex_eigenvalue.cpp @@ -29,7 +29,7 @@ main() 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( + 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 diff --git a/include/proxsuite/proxqp/dense/helpers.hpp b/include/proxsuite/proxqp/dense/helpers.hpp index d58857bf9..b41b649a8 100644 --- a/include/proxsuite/proxqp/dense/helpers.hpp +++ b/include/proxsuite/proxqp/dense/helpers.hpp @@ -21,32 +21,39 @@ namespace proxsuite { namespace proxqp { namespace dense { -template +template T -power_iteration(MatRef H, - VecRefMut dw, - VecRefMut rhs, - VecRefMut err_v, +power_iteration(const Eigen::MatrixBase& H, + const Eigen::MatrixBase& dw, + const Eigen::MatrixBase& rhs, + const Eigen::MatrixBase& err_v, T power_iteration_accuracy, isize nb_power_iteration) { + auto& dw_cc = dw.const_cast_derived(); + auto& rhs_cc = rhs.const_cast_derived(); + auto& err_v_cc = err_v.const_cast_derived(); // computes maximal eigen value of the bottom right matrix of the LDLT isize dim = H.rows(); - rhs.setZero(); + rhs_cc.setZero(); // stores eigenvector - rhs.array() += 1. / std::sqrt(dim); + rhs_cc.array() += 1. / std::sqrt(dim); // stores Hx - dw.noalias() = H.template selfadjointView() * rhs; // Hx + dw_cc.noalias() = H.template selfadjointView() * rhs_cc; // Hx T eig = 0; for (isize i = 0; i < nb_power_iteration; i++) { - rhs = dw / dw.norm(); - dw.noalias() = (H.template selfadjointView() * rhs); + rhs_cc = dw_cc / dw_cc.norm(); + dw_cc.noalias() = (H.template selfadjointView() * rhs_cc); // calculate associated eigenvalue - eig = rhs.dot(dw); + eig = rhs.dot(dw_cc); // calculate associated error - err_v = dw - eig * rhs; - T err = proxsuite::proxqp::dense::infty_norm(err_v); + err_v_cc = dw_cc - eig * rhs_cc; + T err = proxsuite::proxqp::dense::infty_norm(err_v_cc); // std::cout << "power iteration max: i " << i << " err " << err << // std::endl; if (err <= power_iteration_accuracy) { @@ -55,37 +62,46 @@ power_iteration(MatRef H, } return eig; } -template +template T -min_eigen_value_via_modified_power_iteration(MatRef H, - VecRefMut dw, - VecRefMut rhs, - VecRefMut err_v, - T max_eigen_value, - T power_iteration_accuracy, - isize nb_power_iteration) +min_eigen_value_via_modified_power_iteration( + const Eigen::MatrixBase& H, + const Eigen::MatrixBase& dw, + const Eigen::MatrixBase& rhs, + const Eigen::MatrixBase& 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 + auto& dw_cc = dw.const_cast_derived(); + auto& rhs_cc = rhs.const_cast_derived(); + auto& err_v_cc = err_v.const_cast_derived(); isize dim = H.rows(); - rhs.setZero(); + rhs_cc.setZero(); // stores eigenvector - rhs.array() += 1. / std::sqrt(dim); + rhs_cc.array() += 1. / std::sqrt(dim); // stores Hx - dw.noalias() = -(H.template selfadjointView() * rhs); // Hx - dw += max_eigen_value * rhs; + dw_cc.noalias() = + -(H.template selfadjointView() * rhs_cc); // Hx + dw_cc += max_eigen_value * rhs_cc; 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; + rhs_cc = dw_cc / dw_cc.norm(); + dw_cc.noalias() = -(H.template selfadjointView() * rhs_cc); + dw_cc += max_eigen_value * rhs_cc; // calculate associated eigenvalue - eig = rhs.dot(dw); + eig = rhs_cc.dot(dw_cc); // calculate associated error - err_v = dw - eig * rhs; - T err = proxsuite::proxqp::dense::infty_norm(err_v); + err_v_cc = dw_cc - eig * rhs_cc; + T err = proxsuite::proxqp::dense::infty_norm(err_v_cc); // std::cout << "power iteration min: i " << i << " err " << err << // std::endl; if (err <= power_iteration_accuracy) { @@ -104,10 +120,10 @@ min_eigen_value_via_modified_power_iteration(MatRef H, * @param nb_power_iteration maximal number of power iteration executed * */ -template +template T estimate_minimal_eigen_value_of_symmetric_matrix( - MatRef H, + const Eigen::MatrixBase& H, EigenValueEstimateMethodOption estimate_method_option, T power_iteration_accuracy, isize nb_power_iteration) @@ -129,16 +145,16 @@ estimate_minimal_eigen_value_of_symmetric_matrix( Vec dw(dim); Vec rhs(dim); Vec err(dim); - T dominant_eigen_value = power_iteration( + 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); + 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: { diff --git a/test/src/dense_qp_wrapper.cpp b/test/src/dense_qp_wrapper.cpp index 00fa5b247..2c826821a 100644 --- a/test/src/dense_qp_wrapper.cpp +++ b/test/src/dense_qp_wrapper.cpp @@ -7229,7 +7229,7 @@ TEST_CASE("ProxQP::dense: estimate of minimal eigenvalues using Eigen") qp_random.H.diagonal().tail(1).setConstant(-1.); T estimate_minimal_eigen_value = - dense::estimate_minimal_eigen_value_of_symmetric_matrix( + 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); @@ -7266,7 +7266,7 @@ TEST_CASE("ProxQP::dense: estimate of minimal eigenvalues using Eigen") T minimal_eigenvalue = qp_random.H.diagonal().minCoeff(); T estimate_minimal_eigen_value = - dense::estimate_minimal_eigen_value_of_symmetric_matrix( + 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); @@ -7303,7 +7303,7 @@ TEST_CASE("ProxQP::dense: estimate of minimal eigenvalues using Eigen") T minimal_eigenvalue = T(es.eigenvalues().minCoeff()); T estimate_minimal_eigen_value = - dense::estimate_minimal_eigen_value_of_symmetric_matrix( + 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); @@ -7457,7 +7457,7 @@ TEST_CASE( qp_random.H.diagonal().tail(1).setConstant(-0.5); T estimate_minimal_eigen_value = - dense::estimate_minimal_eigen_value_of_symmetric_matrix( + dense::estimate_minimal_eigen_value_of_symmetric_matrix( qp_random.H, EigenValueEstimateMethodOption::PowerIteration, 1.E-6, @@ -7497,7 +7497,7 @@ TEST_CASE( T minimal_eigenvalue = qp_random.H.diagonal().minCoeff(); T estimate_minimal_eigen_value = - dense::estimate_minimal_eigen_value_of_symmetric_matrix( + dense::estimate_minimal_eigen_value_of_symmetric_matrix( qp_random.H, EigenValueEstimateMethodOption::PowerIteration, 1.E-6, @@ -7538,7 +7538,7 @@ TEST_CASE( T minimal_eigenvalue = T(es.eigenvalues().minCoeff()); T estimate_minimal_eigen_value = - dense::estimate_minimal_eigen_value_of_symmetric_matrix( + dense::estimate_minimal_eigen_value_of_symmetric_matrix( qp_random.H, EigenValueEstimateMethodOption::PowerIteration, 1.E-6, @@ -7588,4 +7588,30 @@ DOCTEST_TEST_CASE("check that model.is_valid function for symmetric matrices " // model.is_valid() that performs the check above proxqp::dense::QP qp(3, 0, 0); qp.init(symmetric_mat, nullopt, nullopt, nullopt, nullopt, nullopt, nullopt); +} + +TEST_CASE("ProxQP::dense: test memory allocation when estimating biggest " + "eigenvalue with 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); + Eigen::Matrix H; + Eigen::VectorXd dw(2), rhs(2), err_v(2); + // trivial test + ::proxsuite::proxqp::utils::rand::set_seed(1234); + 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); + H = qp_random.H; + PROXSUITE_EIGEN_MALLOC_NOT_ALLOWED(); + dense::power_iteration(H, dw, rhs, err_v, 1.E-6, 10000); + PROXSUITE_EIGEN_MALLOC_ALLOWED(); } \ No newline at end of file