Skip to content

Commit c5f1081

Browse files
quentinlljcarpent
quentinll
authored and
jcarpent
committed
templating power_iteration
1 parent f0cfea0 commit c5f1081

File tree

4 files changed

+97
-75
lines changed

4 files changed

+97
-75
lines changed

‎bindings/python/src/expose-helpers.hpp‎

Lines changed: 29 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -20,25 +20,35 @@ template<typename T>
2020
void
2121
exposeDenseHelpers(pybind11::module_ m)
2222
{
23-
m.def("estimate_minimal_eigen_value_of_symmetric_matrix",
24-
&dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>,
25-
"Function for estimating the minimal eigenvalue of a dense symmetric "
26-
"matrix. "
27-
"Two options are available: an exact method using "
28-
"SelfAdjointEigenSolver from Eigen, "
29-
"or a Power Iteration algorithm (with parameters : "
30-
"power_iteration_accuracy and nb_power_iteration).",
31-
pybind11::arg("H"),
32-
pybind11::arg_v("estimate_method_option",
33-
EigenValueEstimateMethodOption::ExactMethod,
34-
"Two options are available for "
35-
"estimating smallest eigenvalue: either a power "
36-
"iteration algorithm, or an exact method from Eigen."),
37-
pybind11::arg_v(
38-
"power_iteration_accuracy", T(1.E-3), "power iteration accuracy."),
39-
pybind11::arg_v("nb_power_iteration",
40-
1000,
41-
"maximal number of power iteration executed."));
23+
m.def(
24+
"estimate_minimal_eigen_value_of_symmetric_matrix",
25+
+[](const MatRef<T>& H,
26+
EigenValueEstimateMethodOption estimate_method_option,
27+
T power_iteration_accuracy,
28+
isize nb_power_iteration) {
29+
return dense::estimate_minimal_eigen_value_of_symmetric_matrix(
30+
H,
31+
estimate_method_option,
32+
power_iteration_accuracy,
33+
nb_power_iteration);
34+
},
35+
"Function for estimating the minimal eigenvalue of a dense symmetric "
36+
"matrix. "
37+
"Two options are available: an exact method using "
38+
"SelfAdjointEigenSolver from Eigen, "
39+
"or a Power Iteration algorithm (with parameters : "
40+
"power_iteration_accuracy and nb_power_iteration).",
41+
pybind11::arg("H"),
42+
pybind11::arg_v("estimate_method_option",
43+
EigenValueEstimateMethodOption::ExactMethod,
44+
"Two options are available for "
45+
"estimating smallest eigenvalue: either a power "
46+
"iteration algorithm, or an exact method from Eigen."),
47+
pybind11::arg_v(
48+
"power_iteration_accuracy", T(1.E-3), "power iteration accuracy."),
49+
pybind11::arg_v("nb_power_iteration",
50+
1000,
51+
"maximal number of power iteration executed."));
4252
}
4353
} // namespace python
4454
} // namespace dense

‎examples/cpp/estimate_nonconvex_eigenvalue.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ main()
2929
dense::QP<T> qp(dim, n_eq, n_in); // create the QP object
3030
// choose the option for estimating this eigenvalue
3131
T estimate_minimal_eigen_value =
32-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T, Eigen::RowMajor>(
32+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
3333
qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000);
3434
bool compute_preconditioner = false;
3535
// input the estimate for making rho appropriate

‎include/proxsuite/proxqp/dense/helpers.hpp‎

Lines changed: 58 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -21,32 +21,39 @@ namespace proxsuite {
2121
namespace proxqp {
2222
namespace dense {
2323

24-
template<typename T, int l>
24+
template<typename T,
25+
typename MatIn,
26+
typename VecIn1,
27+
typename VecIn2,
28+
typename VecIn3>
2529
T
26-
power_iteration(MatRef<T, l> H,
27-
VecRefMut<T> dw,
28-
VecRefMut<T> rhs,
29-
VecRefMut<T> err_v,
30+
power_iteration(const Eigen::MatrixBase<MatIn>& H,
31+
const Eigen::MatrixBase<VecIn1>& dw,
32+
const Eigen::MatrixBase<VecIn2>& rhs,
33+
const Eigen::MatrixBase<VecIn3>& err_v,
3034
T power_iteration_accuracy,
3135
isize nb_power_iteration)
3236
{
37+
auto& dw_cc = dw.const_cast_derived();
38+
auto& rhs_cc = rhs.const_cast_derived();
39+
auto& err_v_cc = err_v.const_cast_derived();
3340
// computes maximal eigen value of the bottom right matrix of the LDLT
3441
isize dim = H.rows();
35-
rhs.setZero();
42+
rhs_cc.setZero();
3643
// stores eigenvector
37-
rhs.array() += 1. / std::sqrt(dim);
44+
rhs_cc.array() += 1. / std::sqrt(dim);
3845
// stores Hx
39-
dw.noalias() = H.template selfadjointView<Eigen::Lower>() * rhs; // Hx
46+
dw_cc.noalias() = H.template selfadjointView<Eigen::Lower>() * rhs_cc; // Hx
4047
T eig = 0;
4148
for (isize i = 0; i < nb_power_iteration; i++) {
4249

43-
rhs = dw / dw.norm();
44-
dw.noalias() = (H.template selfadjointView<Eigen::Lower>() * rhs);
50+
rhs_cc = dw_cc / dw_cc.norm();
51+
dw_cc.noalias() = (H.template selfadjointView<Eigen::Lower>() * rhs_cc);
4552
// calculate associated eigenvalue
46-
eig = rhs.dot(dw);
53+
eig = rhs.dot(dw_cc);
4754
// calculate associated error
48-
err_v = dw - eig * rhs;
49-
T err = proxsuite::proxqp::dense::infty_norm(err_v);
55+
err_v_cc = dw_cc - eig * rhs_cc;
56+
T err = proxsuite::proxqp::dense::infty_norm(err_v_cc);
5057
// std::cout << "power iteration max: i " << i << " err " << err <<
5158
// std::endl;
5259
if (err <= power_iteration_accuracy) {
@@ -55,37 +62,46 @@ power_iteration(MatRef<T, l> H,
5562
}
5663
return eig;
5764
}
58-
template<typename T, int l>
65+
template<typename T,
66+
typename MatIn,
67+
typename VecIn1,
68+
typename VecIn2,
69+
typename VecIn3>
5970
T
60-
min_eigen_value_via_modified_power_iteration(MatRef<T, l> H,
61-
VecRefMut<T> dw,
62-
VecRefMut<T> rhs,
63-
VecRefMut<T> err_v,
64-
T max_eigen_value,
65-
T power_iteration_accuracy,
66-
isize nb_power_iteration)
71+
min_eigen_value_via_modified_power_iteration(
72+
const Eigen::MatrixBase<MatIn>& H,
73+
const Eigen::MatrixBase<VecIn1>& dw,
74+
const Eigen::MatrixBase<VecIn2>& rhs,
75+
const Eigen::MatrixBase<VecIn3>& err_v,
76+
T max_eigen_value,
77+
T power_iteration_accuracy,
78+
isize nb_power_iteration)
6779
{
6880
// performs power iteration on the matrix: max_eigen_value I - H
6981
// estimates then the minimal eigenvalue with: minimal_eigenvalue =
7082
// max_eigen_value - eig
83+
auto& dw_cc = dw.const_cast_derived();
84+
auto& rhs_cc = rhs.const_cast_derived();
85+
auto& err_v_cc = err_v.const_cast_derived();
7186
isize dim = H.rows();
72-
rhs.setZero();
87+
rhs_cc.setZero();
7388
// stores eigenvector
74-
rhs.array() += 1. / std::sqrt(dim);
89+
rhs_cc.array() += 1. / std::sqrt(dim);
7590
// stores Hx
76-
dw.noalias() = -(H.template selfadjointView<Eigen::Lower>() * rhs); // Hx
77-
dw += max_eigen_value * rhs;
91+
dw_cc.noalias() =
92+
-(H.template selfadjointView<Eigen::Lower>() * rhs_cc); // Hx
93+
dw_cc += max_eigen_value * rhs_cc;
7894
T eig = 0;
7995
for (isize i = 0; i < nb_power_iteration; i++) {
8096

81-
rhs = dw / dw.norm();
82-
dw.noalias() = -(H.template selfadjointView<Eigen::Lower>() * rhs);
83-
dw += max_eigen_value * rhs;
97+
rhs_cc = dw_cc / dw_cc.norm();
98+
dw_cc.noalias() = -(H.template selfadjointView<Eigen::Lower>() * rhs_cc);
99+
dw_cc += max_eigen_value * rhs_cc;
84100
// calculate associated eigenvalue
85-
eig = rhs.dot(dw);
101+
eig = rhs_cc.dot(dw_cc);
86102
// calculate associated error
87-
err_v = dw - eig * rhs;
88-
T err = proxsuite::proxqp::dense::infty_norm(err_v);
103+
err_v_cc = dw_cc - eig * rhs_cc;
104+
T err = proxsuite::proxqp::dense::infty_norm(err_v_cc);
89105
// std::cout << "power iteration min: i " << i << " err " << err <<
90106
// std::endl;
91107
if (err <= power_iteration_accuracy) {
@@ -104,10 +120,10 @@ min_eigen_value_via_modified_power_iteration(MatRef<T, l> H,
104120
* @param nb_power_iteration maximal number of power iteration executed
105121
*
106122
*/
107-
template<typename T, int l>
123+
template<typename T, typename MatIn>
108124
T
109125
estimate_minimal_eigen_value_of_symmetric_matrix(
110-
MatRef<T, l> H,
126+
const Eigen::MatrixBase<MatIn>& H,
111127
EigenValueEstimateMethodOption estimate_method_option,
112128
T power_iteration_accuracy,
113129
isize nb_power_iteration)
@@ -129,16 +145,16 @@ estimate_minimal_eigen_value_of_symmetric_matrix(
129145
Vec<T> dw(dim);
130146
Vec<T> rhs(dim);
131147
Vec<T> err(dim);
132-
T dominant_eigen_value = power_iteration<T>(
148+
T dominant_eigen_value = power_iteration(
133149
H, dw, rhs, err, power_iteration_accuracy, nb_power_iteration);
134-
T min_eigenvalue = min_eigen_value_via_modified_power_iteration<T>(
135-
H,
136-
dw,
137-
rhs,
138-
err,
139-
dominant_eigen_value,
140-
power_iteration_accuracy,
141-
nb_power_iteration);
150+
T min_eigenvalue =
151+
min_eigen_value_via_modified_power_iteration(H,
152+
dw,
153+
rhs,
154+
err,
155+
dominant_eigen_value,
156+
power_iteration_accuracy,
157+
nb_power_iteration);
142158
res = std::min(min_eigenvalue, dominant_eigen_value);
143159
} break;
144160
case EigenValueEstimateMethodOption::ExactMethod: {

‎test/src/dense_qp_wrapper.cpp‎

Lines changed: 9 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,6 @@ using T = double;
1313
using namespace proxsuite;
1414
using namespace proxsuite::proxqp;
1515

16-
/*
1716
DOCTEST_TEST_CASE(
1817
"ProxQP::dense: sparse random strongly convex qp with inequality constraints"
1918
"and empty equality constraints")
@@ -7230,7 +7229,7 @@ TEST_CASE("ProxQP::dense: estimate of minimal eigenvalues using Eigen")
72307229
qp_random.H.diagonal().tail(1).setConstant(-1.);
72317230

72327231
T estimate_minimal_eigen_value =
7233-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>(
7232+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
72347233
qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000);
72357234

72367235
proxqp::dense::QP<T> qp(dim, n_eq, n_in);
@@ -7267,7 +7266,7 @@ TEST_CASE("ProxQP::dense: estimate of minimal eigenvalues using Eigen")
72677266
T minimal_eigenvalue = qp_random.H.diagonal().minCoeff();
72687267

72697268
T estimate_minimal_eigen_value =
7270-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>(
7269+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
72717270
qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000);
72727271

72737272
proxqp::dense::QP<T> qp(dim, n_eq, n_in);
@@ -7304,7 +7303,7 @@ TEST_CASE("ProxQP::dense: estimate of minimal eigenvalues using Eigen")
73047303
T minimal_eigenvalue = T(es.eigenvalues().minCoeff());
73057304

73067305
T estimate_minimal_eigen_value =
7307-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>(
7306+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
73087307
qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000);
73097308

73107309
proxqp::dense::QP<T> qp(dim, n_eq, n_in);
@@ -7458,7 +7457,7 @@ TEST_CASE(
74587457
qp_random.H.diagonal().tail(1).setConstant(-0.5);
74597458

74607459
T estimate_minimal_eigen_value =
7461-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>(
7460+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
74627461
qp_random.H,
74637462
EigenValueEstimateMethodOption::PowerIteration,
74647463
1.E-6,
@@ -7498,7 +7497,7 @@ TEST_CASE(
74987497
T minimal_eigenvalue = qp_random.H.diagonal().minCoeff();
74997498

75007499
T estimate_minimal_eigen_value =
7501-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>(
7500+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
75027501
qp_random.H,
75037502
EigenValueEstimateMethodOption::PowerIteration,
75047503
1.E-6,
@@ -7539,7 +7538,7 @@ TEST_CASE(
75397538
T minimal_eigenvalue = T(es.eigenvalues().minCoeff());
75407539

75417540
T estimate_minimal_eigen_value =
7542-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>(
7541+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
75437542
qp_random.H,
75447543
EigenValueEstimateMethodOption::PowerIteration,
75457544
1.E-6,
@@ -7591,10 +7590,8 @@ DOCTEST_TEST_CASE("check that model.is_valid function for symmetric matrices "
75917590
qp.init(symmetric_mat, nullopt, nullopt, nullopt, nullopt, nullopt, nullopt);
75927591
}
75937592

7594-
*/
7595-
7596-
TEST_CASE(
7597-
"ProxQP::dense: test memory allocation when estimating biggest eigenvalue with power iteration")
7593+
TEST_CASE("ProxQP::dense: test memory allocation when estimating biggest "
7594+
"eigenvalue with power iteration")
75987595
{
75997596
double sparsity_factor = 1.;
76007597
T tol = T(1e-3);
@@ -7615,7 +7612,6 @@ TEST_CASE(
76157612
qp_random.H.diagonal().tail(1).setConstant(-0.5);
76167613
H = qp_random.H;
76177614
PROXSUITE_EIGEN_MALLOC_NOT_ALLOWED();
7618-
dense::power_iteration<T, Eigen::ColMajor>(
7619-
H, dw, rhs, err_v, 1.E-6, 10000);
7615+
dense::power_iteration(H, dw, rhs, err_v, 1.E-6, 10000);
76207616
PROXSUITE_EIGEN_MALLOC_ALLOWED();
76217617
}

0 commit comments

Comments
 (0)