diff --git a/include/preprocess/max_inscribed_ellipsoid.hpp b/include/preprocess/max_inscribed_ellipsoid.hpp index ebdb261e0..7a0c19a6f 100644 --- a/include/preprocess/max_inscribed_ellipsoid.hpp +++ b/include/preprocess/max_inscribed_ellipsoid.hpp @@ -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. @@ -15,6 +15,9 @@ #include #include +#include "Spectra/include/Spectra/SymEigsSolver.h" +#include "Spectra/include/Spectra/Util/SelectionRule.h" +#include "Spectra/include/Spectra/MatOp/DenseSymMatProd.h" /* @@ -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 std::pair, 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 Diagonal_MT; @@ -49,8 +53,8 @@ std::pair, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT last_r1 = std::numeric_limits::lowest(), last_r2 = std::numeric_limits::lowest(), prev_obj = std::numeric_limits::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); @@ -74,9 +78,10 @@ std::pair, 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 llt(E2); - Q.noalias() = A * E2 * A_trans; + Q.noalias() = A * llt.solve(A_trans); h = Q.diagonal(); h = h.cwiseSqrt(); @@ -115,19 +120,38 @@ std::pair, 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 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 op(E2); + // The value of ncv is chosen empirically + Spectra::SymEigsSolver> 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 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; @@ -138,9 +162,10 @@ std::pair, 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::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::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; @@ -151,7 +176,7 @@ std::pair, 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(); @@ -165,10 +190,9 @@ std::pair, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT G.diagonal() += y2h_z; h_z = h + z; + Eigen::PartialPivLU 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(); @@ -184,11 +208,11 @@ std::pair, 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); @@ -228,9 +252,13 @@ std::pair, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT i++; } - return std::pair, bool>(std::pair(E2, x), converged); + if (!converged) { + x += x0; + } + return std::pair, bool>(std::pair(E2, x), converged); } #endif + diff --git a/include/preprocess/max_inscribed_ellipsoid_rounding.hpp b/include/preprocess/max_inscribed_ellipsoid_rounding.hpp index c13d53f70..da00a8867 100644 --- a/include/preprocess/max_inscribed_ellipsoid_rounding.hpp +++ b/include/preprocess/max_inscribed_ellipsoid_rounding.hpp @@ -22,14 +22,16 @@ template typename Point > std::tuple 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, 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; @@ -44,18 +46,28 @@ std::tuple 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 lltOfA(E); // compute the Cholesky decomposition of E + Eigen::LLT lltOfA(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of E^{-1} L = lltOfA.matrixL(); - Eigen::SelfAdjointEigenSolver 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 op(E); + // The value of ncv is chosen empirically + Spectra::SymEigsSolver> 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 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; @@ -67,6 +79,11 @@ std::tuple 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++; } diff --git a/test/new_rounding_test.cpp b/test/new_rounding_test.cpp index 87c63a32a..ab17ed939 100644 --- a/test/new_rounding_test.cpp +++ b/test/new_rounding_test.cpp @@ -26,6 +26,7 @@ #include "preprocess/svd_rounding.hpp" #include "generators/known_polytope_generators.h" +#include "generators/h_polytopes_generator.h" template NT factorial(NT n) @@ -108,7 +109,6 @@ void rounding_max_ellipsoid_test(Polytope &HP, typedef BoostRandomNumberGenerator RNGType; RNGType rng(d); - std::pair InnerBall = HP.ComputeInnerBall(); std::tuple res = max_inscribed_ellipsoid_rounding(HP, InnerBall.first); @@ -184,7 +184,7 @@ void call_test_max_ellipsoid() { std::cout << "\n--- Testing rounding of H-skinny_cube5" << std::endl; P = generate_skinny_cube(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); }