Skip to content

Commit

Permalink
Improve max volume ellipsoid computation (GeomScale#309)
Browse files Browse the repository at this point in the history 
* improve the implementation of maximum ellipsoid computation

* minor fix in rounding unit-test

* fix file copyrights

* apply termination criterion after transformation in max ellipsoid rounding

* resolve PR comments
  • Loading branch information
@TolisChal
TolisChal committed Jul 8, 2024
1 parent c3af466 commit ff24350
Show file tree
Hide file tree
Showing 3 changed files with 83 additions and 38 deletions.
78 changes: 53 additions & 25 deletions include/preprocess/max_inscribed_ellipsoid.hpp
View file
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
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;
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)));
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
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

0 comments on commit ff24350

Please sign in to comment.