Skip to content

Commit 87d35ab

Browse files
jcarpentjcarpent
authored
Merge pull request #279 from quentinll/qll/test-colmajor
Templating power iteration algorithm by matrix storage order
2 parents b2cc776 + c5f1081 commit 87d35ab

File tree

4 files changed

+120
-68
lines changed

4 files changed

+120
-68
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>(
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>
24+
template<typename T,
25+
typename MatIn,
26+
typename VecIn1,
27+
typename VecIn2,
28+
typename VecIn3>
2529
T
26-
power_iteration(MatRef<T> 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> H,
5562
}
5663
return eig;
5764
}
58-
template<typename T>
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> 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> H,
104120
* @param nb_power_iteration maximal number of power iteration executed
105121
*
106122
*/
107-
template<typename T>
123+
template<typename T, typename MatIn>
108124
T
109125
estimate_minimal_eigen_value_of_symmetric_matrix(
110-
MatRef<T> 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: 32 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -7229,7 +7229,7 @@ TEST_CASE("ProxQP::dense: estimate of minimal eigenvalues using Eigen")
72297229
qp_random.H.diagonal().tail(1).setConstant(-1.);
72307230

72317231
T estimate_minimal_eigen_value =
7232-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>(
7232+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
72337233
qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000);
72347234

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

72687268
T estimate_minimal_eigen_value =
7269-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>(
7269+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
72707270
qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000);
72717271

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

73057305
T estimate_minimal_eigen_value =
7306-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>(
7306+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
73077307
qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000);
73087308

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

74597459
T estimate_minimal_eigen_value =
7460-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>(
7460+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
74617461
qp_random.H,
74627462
EigenValueEstimateMethodOption::PowerIteration,
74637463
1.E-6,
@@ -7497,7 +7497,7 @@ TEST_CASE(
74977497
T minimal_eigenvalue = qp_random.H.diagonal().minCoeff();
74987498

74997499
T estimate_minimal_eigen_value =
7500-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>(
7500+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
75017501
qp_random.H,
75027502
EigenValueEstimateMethodOption::PowerIteration,
75037503
1.E-6,
@@ -7538,7 +7538,7 @@ TEST_CASE(
75387538
T minimal_eigenvalue = T(es.eigenvalues().minCoeff());
75397539

75407540
T estimate_minimal_eigen_value =
7541-
dense::estimate_minimal_eigen_value_of_symmetric_matrix<T>(
7541+
dense::estimate_minimal_eigen_value_of_symmetric_matrix(
75427542
qp_random.H,
75437543
EigenValueEstimateMethodOption::PowerIteration,
75447544
1.E-6,
@@ -7588,4 +7588,30 @@ DOCTEST_TEST_CASE("check that model.is_valid function for symmetric matrices "
75887588
// model.is_valid() that performs the check above
75897589
proxqp::dense::QP<T> qp(3, 0, 0);
75907590
qp.init(symmetric_mat, nullopt, nullopt, nullopt, nullopt, nullopt, nullopt);
7591+
}
7592+
7593+
TEST_CASE("ProxQP::dense: test memory allocation when estimating biggest "
7594+
"eigenvalue with power iteration")
7595+
{
7596+
double sparsity_factor = 1.;
7597+
T tol = T(1e-3);
7598+
utils::rand::set_seed(1);
7599+
dense::isize dim = 2;
7600+
dense::isize n_eq(dim);
7601+
dense::isize n_in(dim);
7602+
T strong_convexity_factor(1.e-2);
7603+
Eigen::Matrix<double, 2, 2, Eigen::ColMajor> H;
7604+
Eigen::VectorXd dw(2), rhs(2), err_v(2);
7605+
// trivial test
7606+
::proxsuite::proxqp::utils::rand::set_seed(1234);
7607+
proxqp::dense::Model<T> qp_random = proxqp::utils::dense_strongly_convex_qp(
7608+
dim, n_eq, n_in, sparsity_factor, strong_convexity_factor);
7609+
7610+
qp_random.H.setZero();
7611+
qp_random.H.diagonal().setOnes();
7612+
qp_random.H.diagonal().tail(1).setConstant(-0.5);
7613+
H = qp_random.H;
7614+
PROXSUITE_EIGEN_MALLOC_NOT_ALLOWED();
7615+
dense::power_iteration(H, dw, rhs, err_v, 1.E-6, 10000);
7616+
PROXSUITE_EIGEN_MALLOC_ALLOWED();
75917617
}

0 commit comments

Comments
 (0)