Skip to content

Commit

Permalink
templating power_iteration
Browse files Browse the repository at this point in the history
  • Loading branch information
@quentinll @fabinsch
quentinll authored and fabinsch committed Nov 13, 2023
1 parent 4a9dd8f commit cbc5297
Show file tree
Hide file tree
Showing 4 changed files with 97 additions and 75 deletions.
48 changes: 29 additions & 19 deletions bindings/python/src/expose-helpers.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion examples/cpp/estimate_nonconvex_eigenvalue.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, Eigen::RowMajor>(
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
Expand Down
100 changes: 58 additions & 42 deletions include/proxsuite/proxqp/dense/helpers.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,32 +21,39 @@ namespace proxsuite {
namespace proxqp {
namespace dense {

template<typename T, int l>
template<typename T,
typename MatIn,
typename VecIn1,
typename VecIn2,
typename VecIn3>
T
power_iteration(MatRef<T, l> 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) {
Expand All @@ -55,37 +62,46 @@ power_iteration(MatRef<T, l> H,
}
return eig;
}
template<typename T, int l>
template<typename T,
typename MatIn,
typename VecIn1,
typename VecIn2,
typename VecIn3>
T
min_eigen_value_via_modified_power_iteration(MatRef<T, l> 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) {
Expand All @@ -104,10 +120,10 @@ min_eigen_value_via_modified_power_iteration(MatRef<T, l> H,
* @param nb_power_iteration maximal number of power iteration executed
*
*/
template<typename T, int l>
template<typename T, typename MatIn>
T
estimate_minimal_eigen_value_of_symmetric_matrix(
MatRef<T, l> H,
const Eigen::MatrixBase<MatIn>& H,
EigenValueEstimateMethodOption estimate_method_option,
T power_iteration_accuracy,
isize nb_power_iteration)
Expand All @@ -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: {
Expand Down
22 changes: 9 additions & 13 deletions test/src/dense_qp_wrapper.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,6 @@ using T = double;
using namespace proxsuite;
using namespace proxsuite::proxqp;

/*
DOCTEST_TEST_CASE(
"ProxQP::dense: sparse random strongly convex qp with inequality constraints"
"and empty equality constraints")
Expand Down Expand Up @@ -7230,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);
Expand Down Expand Up @@ -7267,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);
Expand Down Expand Up @@ -7304,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);
Expand Down Expand Up @@ -7458,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,
Expand Down Expand Up @@ -7498,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,
Expand Down Expand Up @@ -7539,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,
Expand Down Expand Up @@ -7591,10 +7590,8 @@ DOCTEST_TEST_CASE("check that model.is_valid function for symmetric matrices "
qp.init(symmetric_mat, nullopt, nullopt, nullopt, nullopt, nullopt, nullopt);
}

*/

TEST_CASE(
"ProxQP::dense: test memory allocation when estimating biggest eigenvalue with power iteration")
TEST_CASE("ProxQP::dense: test memory allocation when estimating biggest "
"eigenvalue with power iteration")
{
double sparsity_factor = 1.;
T tol = T(1e-3);
Expand All @@ -7615,7 +7612,6 @@ TEST_CASE(
qp_random.H.diagonal().tail(1).setConstant(-0.5);
H = qp_random.H;
PROXSUITE_EIGEN_MALLOC_NOT_ALLOWED();
dense::power_iteration<T, Eigen::ColMajor>(
H, dw, rhs, err_v, 1.E-6, 10000);
dense::power_iteration(H, dw, rhs, err_v, 1.E-6, 10000);
PROXSUITE_EIGEN_MALLOC_ALLOWED();
}

0 comments on commit cbc5297

Please sign in to comment.