Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve max volume ellipsoid computation #309

Merged
merged 5 commits into from
Jun 20, 2024
Merged

Improve max volume ellipsoid computation #309

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
d5cd638
improve the implementation of maximum ellipsoid computation
Jun 18, 2024
f30407e
minor fix in rounding unit-test
Jun 18, 2024
80c0b1d
fix file copyrights
Jun 18, 2024
f238039
apply termination criterion after transformation in max ellipsoid rou…
Jun 19, 2024
a5f18aa
resolve PR comments
Jun 20, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 53 additions & 25 deletions include/preprocess/max_inscribed_ellipsoid.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
// VolEsti (volume computation and sampling library)

// Copyright (c) 2012-2020 Vissarion Fisikopoulos
// Copyright (c) 2018-2020 Apostolos Chalkis
// Copyright (c) 2012-2024 Vissarion Fisikopoulos
// Copyright (c) 2018-2024 Apostolos Chalkis
// Copyright (c) 2021 Vaibhav Thakkar

//Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
Expand All @@ -15,6 +15,9 @@

#include <utility>
#include <Eigen/Eigen>
#include "Spectra/include/Spectra/SymEigsSolver.h"
#include "Spectra/include/Spectra/Util/SelectionRule.h"
#include "Spectra/include/Spectra/MatOp/DenseSymMatProd.h"


/*
Expand All @@ -31,14 +34,15 @@
tolerance parameters tol, reg

Output: center of the ellipsoid y
matrix V = E_transpose * E
matrix E2^{-1} = E_transpose * E
*/

// using Custom_MT as to deal with both dense and sparse matrices, MT will be the type of result matrix
// TODO: Change the return data structure to std::tuple
template <typename MT, typename Custom_MT, typename VT, typename NT>
std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT const& x0,
unsigned int const& maxiter,
NT const& tol, NT const& reg)
unsigned int const& maxiter,
NT const& tol, NT const& reg)
{
typedef Eigen::DiagonalMatrix<NT, Eigen::Dynamic> Diagonal_MT;

Expand All @@ -49,8 +53,8 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT
last_r1 = std::numeric_limits<NT>::lowest(),
last_r2 = std::numeric_limits<NT>::lowest(),
prev_obj = std::numeric_limits<NT>::lowest(),
gap, rmu, res, objval, r1, r2 ,r3, rel, Rel,
astep, ax, ay, az, tau;
gap, rmu, res, objval, r1, r2 ,r3, astep, ax,
ay, az, tau, logdetE2;

NT const reg_lim = std::pow(10.0, -10.0), tau0 = 0.75, minmu = std::pow(10.0, -8.0);

Expand All @@ -74,9 +78,10 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT

Y = y.asDiagonal();

E2.noalias() = MT(A_trans * Y * A).inverse();
E2.noalias() = MT(A_trans * Y * A);
Eigen::LLT<MT> llt(E2);

Q.noalias() = A * E2 * A_trans;
Q.noalias() = A * llt.solve(A_trans);
h = Q.diagonal();
h = h.cwiseSqrt();

Expand Down Expand Up @@ -115,19 +120,38 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT

res = std::max(r1, r2);
res = std::max(res, r3);
objval = std::log(E2.determinant()) / 2.0;

Eigen::SelfAdjointEigenSolver <MT> eigensolver(E2); // E2 is positive definite matrix
// computing eigenvalues of E2
rel = eigensolver.eigenvalues().minCoeff();
Rel = eigensolver.eigenvalues().maxCoeff();
logdetE2 = llt.matrixL().toDenseMatrix().diagonal().array().log().sum();
objval = logdetE2; //logdet of E2 is already divided by 2

if (i % 10 == 0) {

NT rel, Rel;

// computing eigenvalues of E2
Spectra::DenseSymMatProd<NT> op(E2);
// The value of ncv is chosen empirically
Spectra::SymEigsSolver<NT, Spectra::SELECT_EIGENVALUE::BOTH_ENDS,
Spectra::DenseSymMatProd<NT>> eigs(&op, 2, std::min(std::max(10, n/5), n));
eigs.init();
int nconv = eigs.compute();
if (eigs.info() == Spectra::COMPUTATION_INFO::SUCCESSFUL) {
Rel = 1.0 / eigs.eigenvalues().coeff(1);
rel = 1.0 / eigs.eigenvalues().coeff(0);
} else {
Eigen::SelfAdjointEigenSolver<MT> eigensolver(E2); // E2 is positive definite matrix
if (eigensolver.info() == Eigen::ComputationInfo::Success) {
Rel = 1.0 / eigensolver.eigenvalues().coeff(0);
rel = 1.0 / eigensolver.eigenvalues().template tail<1>().value();
} else {
std::runtime_error("Computations failed.");
}
}

if (std::abs((last_r1 - r1) / std::min(NT(std::abs(last_r1)), NT(std::abs(r1)))) < 0.01 &&
std::abs((last_r2 - r2) / std::min(NT(abs(last_r2)), NT(std::abs(r2)))) < 0.01 &&
Rel / rel > 100.0 &&
reg > reg_lim) {

converged = false;
//Stopped making progress
break;
Expand All @@ -138,9 +162,10 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT

// stopping criterion
if ((res < tol * (1.0 + bnrm) && rmu <= minmu) ||
(i > 4 && prev_obj != std::numeric_limits<NT>::lowest() &&
((std::abs(objval - prev_obj) <= tol * objval && std::abs(objval - prev_obj) <= tol * prev_obj) ||
(prev_obj >= (1.0 - tol) * objval || objval <= (1.0 - tol) * prev_obj) ) ) ) {
(i > 1 && prev_obj != std::numeric_limits<NT>::lowest() &&
(std::abs(objval - prev_obj) <= tol * std::min(std::abs(objval), std::abs(prev_obj)) ||
std::abs(prev_obj - objval) <= tol) ) ) {

//converged
x += x0;
converged = true;
Expand All @@ -151,7 +176,7 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT
YQ.noalias() = Y * Q;
G = YQ.cwiseProduct(YQ.transpose());
y2h = 2.0 * yh;
YA = Y * A;
YA.noalias() = Y * A;

vec_iter1 = y2h.data();
vec_iter2 = z.data();
Expand All @@ -165,10 +190,9 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT

G.diagonal() += y2h_z;
h_z = h + z;
Eigen::PartialPivLU<MT> luG(G);
T.noalias() = luG.solve(h_z.asDiagonal()*YA);

for (int j = 0; j < n; ++j) {
T.col(j) = G.colPivHouseholderQr().solve( VT(YA.col(j).cwiseProduct(h_z)) );
}
ATP.noalias() = MT(y2h.asDiagonal()*T - YA).transpose();

vec_iter1 = R3.data();
Expand All @@ -184,11 +208,11 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT
R23 = R2 - R3Dy;
ATP_A.noalias() = ATP * A;
ATP_A.diagonal() += ones_m * reg;
dx = ATP_A.colPivHouseholderQr().solve(R1 + ATP * R23); // predictor step
dx = ATP_A.lu().solve(R1 + ATP * R23); // predictor step

// corrector and combined step & length
Adx.noalias() = A * dx;
dyDy = G.colPivHouseholderQr().solve(y2h.cwiseProduct(Adx-R23));
dyDy = luG.solve(y2h.cwiseProduct(Adx-R23));

dy = y.cwiseProduct(dyDy);
dz = R3Dy - z.cwiseProduct(dyDy);
Expand Down Expand Up @@ -228,9 +252,13 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT
i++;
}

return std::pair<std::pair<MT, VT>, bool>(std::pair<MT, VT>(E2, x), converged);
if (!converged) {
x += x0;
}

return std::pair<std::pair<MT, VT>, bool>(std::pair<MT, VT>(E2, x), converged);
}


#endif

39 changes: 28 additions & 11 deletions include/preprocess/max_inscribed_ellipsoid_rounding.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,14 +22,16 @@ template
typename Point
>
std::tuple<MT, VT, NT> max_inscribed_ellipsoid_rounding(Polytope &P,
Point const& InnerPoint)
Point const& InnerPoint,
unsigned int const max_iterations = 5,
NT const max_eig_ratio = NT(6))
{
std::pair<std::pair<MT, VT>, bool> iter_res;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can add a todo here to change this to std::tuple.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, I added a comment in max_inscribed_ellipsoid.hpp

iter_res.second = false;

VT x0 = InnerPoint.getCoefficients();
MT E, L;
unsigned int maxiter = 150, iter = 1, d = P.dimension();
unsigned int maxiter = 500, iter = 1, d = P.dimension();

NT R = 100.0, r = 1.0, tol = std::pow(10, -6.0), reg = std::pow(10, -4.0), round_val = 1.0;

Expand All @@ -44,18 +46,28 @@ std::tuple<MT, VT, NT> max_inscribed_ellipsoid_rounding(Polytope &P,
E = (E + E.transpose()) / 2.0;
E = E + MT::Identity(d, d)*std::pow(10, -8.0); //normalize E

Eigen::LLT<MT> lltOfA(E); // compute the Cholesky decomposition of E
Eigen::LLT<MT> lltOfA(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of E^{-1}
L = lltOfA.matrixL();

Eigen::SelfAdjointEigenSolver <MT> eigensolver(L);
r = eigensolver.eigenvalues().minCoeff();
R = eigensolver.eigenvalues().maxCoeff();

// check the roundness of the polytope
if(((std::abs(R / r) <= 2.3 && iter_res.second) || iter >= 20) && iter>2){
break;
// computing eigenvalues of E
Spectra::DenseSymMatProd<NT> op(E);
// The value of ncv is chosen empirically
Spectra::SymEigsSolver<NT, Spectra::SELECT_EIGENVALUE::BOTH_ENDS,
Spectra::DenseSymMatProd<NT>> eigs(&op, 2, std::min(std::max(10, int(d)/5), int(d)));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

could we add a comment on the numeric values here and how we select them?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The values is chosen empirically. I added a comment.

eigs.init();
int nconv = eigs.compute();
if (eigs.info() == Spectra::COMPUTATION_INFO::SUCCESSFUL) {
R = 1.0 / eigs.eigenvalues().coeff(1);
r = 1.0 / eigs.eigenvalues().coeff(0);
} else {
Eigen::SelfAdjointEigenSolver<MT> eigensolver(E);
if (eigensolver.info() == Eigen::ComputationInfo::Success) {
R = 1.0 / eigensolver.eigenvalues().coeff(0);
r = 1.0 / eigensolver.eigenvalues().template tail<1>().value();
} else {
std::runtime_error("Computations failed.");
}
}

// shift polytope and apply the linear transformation on P
P.shift(iter_res.first.second);
shift += T * iter_res.first.second;
Expand All @@ -67,6 +79,11 @@ std::tuple<MT, VT, NT> max_inscribed_ellipsoid_rounding(Polytope &P,
P.normalize();
x0 = VT::Zero(d);

// check the roundness of the polytope
if(((std::abs(R / r) <= max_eig_ratio && iter_res.second) || iter >= max_iterations)){
break;
}

iter++;
}

Expand Down
4 changes: 2 additions & 2 deletions test/new_rounding_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
#include "preprocess/svd_rounding.hpp"

#include "generators/known_polytope_generators.h"
#include "generators/h_polytopes_generator.h"

template <typename NT>
NT factorial(NT n)
Expand Down Expand Up @@ -108,7 +109,6 @@ void rounding_max_ellipsoid_test(Polytope &HP,

typedef BoostRandomNumberGenerator<boost::mt19937, NT, 5> RNGType;
RNGType rng(d);

std::pair<Point, NT> InnerBall = HP.ComputeInnerBall();
std::tuple<MT, VT, NT> res = max_inscribed_ellipsoid_rounding<MT, VT, NT>(HP, InnerBall.first);

Expand Down Expand Up @@ -184,7 +184,7 @@ void call_test_max_ellipsoid() {

std::cout << "\n--- Testing rounding of H-skinny_cube5" << std::endl;
P = generate_skinny_cube<Hpolytope>(5);
rounding_max_ellipsoid_test(P, 0, 3070.64, 3188.25, 3140.6, 3200.0);
rounding_max_ellipsoid_test(P, 0, 3070.64, 3188.25, 3262.61, 3200.0);
}


Expand Down
Loading