Templating power iteration algorithm by matrix storage order

bindings/python/src/expose-helpers.hpp

@@ -20,25 +20,35 @@ 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."));
+  m.def(
+    "estimate_minimal_eigen_value_of_symmetric_matrix",
+    +[](const MatRef<T>& 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

examples/cpp/estimate_nonconvex_eigenvalue.cpp

@@ -29,7 +29,7 @@ main()
   dense::QP<T> 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<T>(
+    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

include/proxsuite/proxqp/dense/helpers.hpp

@@ -21,32 +21,39 @@ namespace proxsuite {
 namespace proxqp {
 namespace dense {
 
-template<typename T>
+template<typename T,
+         typename MatIn,
+         typename VecIn1,
+         typename VecIn2,
+         typename VecIn3>
 T
-power_iteration(MatRef<T> H,
-                VecRefMut<T> dw,
-                VecRefMut<T> rhs,
-                VecRefMut<T> err_v,
+power_iteration(const Eigen::MatrixBase<MatIn>& H,
+                const Eigen::MatrixBase<VecIn1>& dw,
+                const Eigen::MatrixBase<VecIn2>& rhs,
+                const Eigen::MatrixBase<VecIn3>& 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<Eigen::Lower>() * rhs; // Hx
+  dw_cc.noalias() = H.template selfadjointView<Eigen::Lower>() * rhs_cc; // Hx
   T eig = 0;
   for (isize i = 0; i < nb_power_iteration; i++) {
 
-    rhs = dw / dw.norm();
-    dw.noalias() = (H.template selfadjointView<Eigen::Lower>() * rhs);
+    rhs_cc = dw_cc / dw_cc.norm();
+    dw_cc.noalias() = (H.template selfadjointView<Eigen::Lower>() * 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<T> H,
   }
   return eig;
 }
-template<typename T>
+template<typename T,
+         typename MatIn,
+         typename VecIn1,
+         typename VecIn2,
+         typename VecIn3>
 T
-min_eigen_value_via_modified_power_iteration(MatRef<T> H,
-                                             VecRefMut<T> dw,
-                                             VecRefMut<T> rhs,
-                                             VecRefMut<T> err_v,
-                                             T max_eigen_value,
-                                             T power_iteration_accuracy,
-                                             isize nb_power_iteration)
+min_eigen_value_via_modified_power_iteration(
+  const Eigen::MatrixBase<MatIn>& H,
+  const Eigen::MatrixBase<VecIn1>& dw,
+  const Eigen::MatrixBase<VecIn2>& rhs,
+  const Eigen::MatrixBase<VecIn3>& 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<Eigen::Lower>() * rhs); // Hx
-  dw += max_eigen_value * rhs;
+  dw_cc.noalias() =
+    -(H.template selfadjointView<Eigen::Lower>() * 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<Eigen::Lower>() * rhs);
-    dw += max_eigen_value * rhs;
+    rhs_cc = dw_cc / dw_cc.norm();
+    dw_cc.noalias() = -(H.template selfadjointView<Eigen::Lower>() * 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<T> H,
  * @param nb_power_iteration maximal number of power iteration executed
  *
  */
-template<typename T>
+template<typename T, typename MatIn>
 T
 estimate_minimal_eigen_value_of_symmetric_matrix(
-  MatRef<T> H,
+  const Eigen::MatrixBase<MatIn>& 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<T> dw(dim);
       Vec<T> rhs(dim);
       Vec<T> err(dim);
-      T dominant_eigen_value = power_iteration<T>(
+      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<T>(
-        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: {

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<T>(
+      dense::estimate_minimal_eigen_value_of_symmetric_matrix(
         qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000);
 
     proxqp::dense::QP<T> 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<T>(
+      dense::estimate_minimal_eigen_value_of_symmetric_matrix(
         qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000);
 
     proxqp::dense::QP<T> 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<T>(
+      dense::estimate_minimal_eigen_value_of_symmetric_matrix(
         qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000);
 
     proxqp::dense::QP<T> 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<T>(
+      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<T>(
+      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<T>(
+      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<T> 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<double, 2, 2, Eigen::ColMajor> H;
+  Eigen::VectorXd dw(2), rhs(2), err_v(2);
+  // trivial test
+  ::proxsuite::proxqp::utils::rand::set_seed(1234);
+  proxqp::dense::Model<T> 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();
 }