Improve max volume ellipsoid computation (GeomScale#309)

TolisChal · TolisChal · commit 6e39fba231a8 · 2024-07-06T20:54:38.000-06:00
* 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
diff --git 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 <utility>
 #include <Eigen/Eigen>
+#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 <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;
 
@@ -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);
 
@@ -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();
 
@@ -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;
@@ -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;
@@ -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();
@@ -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();
@@ -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);
@@ -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
+
diff --git 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<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;
 
@@ -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;
@@ -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++;
     }
 
diff --git 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 <typename NT>
 NT factorial(NT n)
@@ -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);
 
@@ -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);
 }
 
 
