Max inscribed ellipsoid sparse

include/preprocess/max_inscribed_ellipsoid.hpp

@@ -2,16 +2,22 @@
 
 // Copyright (c) 2012-2020 Vissarion Fisikopoulos
 // Copyright (c) 2018-2020 Apostolos Chalkis
+// Copyright (c) 2021 Vaibhav Thakkar
 
 //Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
+//Contributed and/or modified by Vaibhav Thakkar, as part of Google Summer of Code 2021 program.
 
 // Licensed under GNU LGPL.3, see LICENCE file
 
 
 #ifndef MVE_COMPUTATION_HPP
 #define MVE_COMPUTATION_HPP
 
-/* 
+#include <utility>
+#include <Eigen/Eigen>
+
+
+/*
     Implementation of the interior point method to compute the largest inscribed ellipsoid in a
     given convex polytope by "Yin Zhang, An Interior-Point Algorithm for the Maximum-Volume Ellipsoid
     Problem (1999)".
@@ -28,42 +34,47 @@
             matrix V = E_transpose * E
 */
 
-template <typename MT, typename VT, typename NT>
-std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(MT A, VT b, VT const& x0,
-                                                   unsigned int const& maxiter, 
-                                                   NT const& tol, NT const& reg)
+// using Custom_MT as to deal with both dense and sparse matrices, MT will be the type of result matrix
+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)
 {
+    typedef Eigen::DiagonalMatrix<NT, Eigen::Dynamic> Diagonal_MT;
+
     int m = A.rows(), n = A.cols();
     bool converged = false;
 
     NT bnrm = b.norm(),
        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, 
+       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;
 
     NT const reg_lim = std::pow(10.0, -10.0), tau0 = 0.75, minmu = std::pow(10.0, -8.0);
 
     NT *vec_iter1, *vec_iter2, *vec_iter3;
 
     VT x = VT::Zero(n), y = VT::Ones(m), bmAx = VT::Ones(m),
-       h(m), z(m), yz(m), yh(m), R1(n), R2(m), R3(m), y2h(m), y2h_z(m), h_z(m), 
+       h(m), z(m), yz(m), yh(m), R1(n), R2(m), R3(m), y2h(m), y2h_z(m), h_z(m),
        R3Dy(m), R23(m), dx(n), Adx(m), dyDy(m), dy(m), dz(m);
 
     VT const bmAx0 = b - A * x0, ones_m = VT::Ones(m);
 
-    MT Q(m, m), Y(m, m), E2(n, n), YQ(m,m), YA(m, n), G(m,m), T(m,n), ATP(n,m), ATP_A(n,n);
+    MT Q(m, m), E2(n, n), YQ(m,m), G(m,m), T(m,n), ATP(n,m), ATP_A(n,n);
+    Diagonal_MT Y(m);
+    Custom_MT YA(m, n);
 
     A = (ones_m.cwiseProduct(bmAx0.cwiseInverse())).asDiagonal() * A, b = ones_m;
-
-    MT A_trans = A.transpose();
+    Custom_MT A_trans = A.transpose();
 
     int i = 1;
     while (i <= maxiter) {
 
         Y = y.asDiagonal();
-        E2.noalias() = (A_trans * Y * A).inverse();
+
+        E2.noalias() = MT(A_trans * Y * A).inverse();
 
         Q.noalias() = A * E2 * A_trans;
         h = Q.diagonal();
@@ -84,7 +95,7 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(MT A, VT b, VT const&
                 vec_iter3++;
             }
             Q *= (t * t);
-            Y *= (1.0 / (t * t));
+            Y = Y * (1.0 / (t * t));
         }
 
         yz = y.cwiseProduct(z);
@@ -126,7 +137,7 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(MT A, VT b, VT const&
         }
 
         // stopping criterion
-        if ((res < tol * (1.0 + bnrm) && rmu <= minmu) || 
+        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) ) ) ) {
@@ -140,7 +151,7 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(MT A, VT b, VT const&
         YQ.noalias() = Y * Q;
         G = YQ.cwiseProduct(YQ.transpose());
         y2h = 2.0 * yh;
-        YA.noalias() = Y * A;
+        YA = Y * A;
 
         vec_iter1 = y2h.data();
         vec_iter2 = z.data();
@@ -156,9 +167,9 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(MT A, VT b, VT const&
         h_z = h + z;
 
         for (int j = 0; j < n; ++j) {
-            T.col(j) = G.colPivHouseholderQr().solve(YA.col(j).cwiseProduct(h_z));
+            T.col(j) = G.colPivHouseholderQr().solve( VT(YA.col(j).cwiseProduct(h_z)) );
         }
-        ATP.noalias() = (y2h.asDiagonal()*T - YA).transpose();
+        ATP.noalias() = MT(y2h.asDiagonal()*T - YA).transpose();
 
         vec_iter1 = R3.data();
         vec_iter2 = y.data();
@@ -222,4 +233,4 @@ std::pair<std::pair<MT, VT>, bool> max_inscribed_ellipsoid(MT A, VT b, VT const&
 }
 
 
-#endif 
+#endif

include/preprocess/max_inscribed_ellipsoid_rounding.hpp

@@ -39,7 +39,7 @@ std::tuple<MT, VT, NT> max_inscribed_ellipsoid_rounding(Polytope &P,
     while (true)
     {
         // compute the largest inscribed ellipsoid in P centered at x0
-        iter_res = max_inscribed_ellipsoid(P.get_mat(), P.get_vec(), x0, maxiter, tol, reg);
+        iter_res = max_inscribed_ellipsoid<MT>(P.get_mat(), P.get_vec(), x0, maxiter, tol, reg);
         E = iter_res.first.first;
         E = (E + E.transpose()) / 2.0;
         E = E + MT::Identity(d, d)*std::pow(10, -8.0); //normalize E

test/CMakeLists.txt

@@ -252,6 +252,10 @@ else ()
   add_test(NAME new_rounding_test_round_skinny_cube
            COMMAND new_rounding_test -tc=round_skinny_cube)
 
+  add_executable (max_ellipsoid_rounding_test max_ellipsoid_rounding_test.cpp $<TARGET_OBJECTS:test_main>)
+  add_test(NAME max_ellipsoid_rounding_test_round_skinny_cube
+           COMMAND max_ellipsoid_rounding_test -tc=round_skinny_cube)
+
   add_executable (logconcave_sampling_test logconcave_sampling_test.cpp $<TARGET_OBJECTS:test_main>)
   add_test(NAME logconcave_sampling_test_hmc
           COMMAND logconcave_sampling_test -tc=hmc)
@@ -269,6 +273,7 @@ else ()
   TARGET_LINK_LIBRARIES(volume_cb_zonotopes ${LP_SOLVE})
   TARGET_LINK_LIBRARIES(volume_cb_vpoly_intersection_vpoly ${LP_SOLVE})
   TARGET_LINK_LIBRARIES(new_rounding_test ${LP_SOLVE})
+  TARGET_LINK_LIBRARIES(max_ellipsoid_rounding_test ${LP_SOLVE})
   TARGET_LINK_LIBRARIES(mcmc_diagnostics_test ${LP_SOLVE})
   TARGET_LINK_LIBRARIES(benchmarks_sob ${LP_SOLVE})
   TARGET_LINK_LIBRARIES(benchmarks_cg ${LP_SOLVE})

test/max_ellipsoid_rounding_test.cpp

@@ -0,0 +1,120 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2020 Vissarion Fisikopoulos
+// Copyright (c) 2018 Apostolos Chalkis
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#include "doctest.h"
+#include <fstream>
+#include <iostream>
+#include "misc.h"
+#include "random.hpp"
+#include "random/uniform_int.hpp"
+#include "random/normal_distribution.hpp"
+#include "random/uniform_real_distribution.hpp"
+
+#include "random_walks/random_walks.hpp"
+
+#include "volume/volume_sequence_of_balls.hpp"
+#include "volume/volume_cooling_gaussians.hpp"
+#include "volume/volume_cooling_balls.hpp"
+
+#include "preprocess/max_inscribed_ellipsoid_rounding.hpp"
+
+#include "known_polytope_generators.h"
+
+template <typename NT>
+NT factorial(NT n)
+{
+    return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
+}
+
+template <typename NT>
+void test_values(NT volume, NT expected, NT exact)
+{
+    std::cout << "Computed volume " << volume << std::endl;
+    std::cout << "Expected volume = " << expected << std::endl;
+    std::cout << "Relative error (expected) = "
+              << std::abs((volume-expected)/expected) << std::endl;
+    std::cout << "Relative error (exact) = "
+              << std::abs((volume-exact)/exact) << std::endl;
+    CHECK((std::abs((volume - exact)/exact) < 0.2 || 
+            std::abs((volume - expected)/expected) < 0.00001));
+}
+
+template <class Polytope>
+void rounding_test(Polytope &HP,
+                 double const& expectedBall,
+                 double const& expectedCDHR,
+                 double const& expectedRDHR,
+                 double const& expectedBilliard,
+                 double const& exact)
+{
+    typedef typename Polytope::PointType Point;
+    typedef typename Point::FT NT;
+    typedef typename Polytope::MT MT;
+    typedef typename Polytope::VT VT;
+
+    int d = HP.dimension();
+
+    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);
+
+    // Setup the parameters
+    int walk_len = 10;
+    NT e = 0.1;
+
+    // Estimate the volume
+    std::cout << "Number type: " << typeid(NT).name() << std::endl;
+
+
+    //TODO: low accuracy in high dimensions
+    // NT volume = res.second * volume_cooling_balls<BallWalk, RNGType>(HP, e, walk_len);
+    // test_values(volume, expectedBall, exact);
+
+    NT volume = std::get<2>(res) * volume_cooling_balls<CDHRWalk, RNGType>(HP, e, walk_len).second;
+    test_values(volume, expectedCDHR, exact);
+
+    volume = std::get<2>(res) * volume_cooling_balls<RDHRWalk, RNGType>(HP, e, 2*walk_len).second;
+    test_values(volume, expectedRDHR, exact);
+
+    volume = std::get<2>(res) * volume_cooling_balls<BilliardWalk, RNGType>(HP, e, walk_len).second;
+    test_values(volume, expectedBilliard, exact);
+}
+
+
+
+
+template <typename NT>
+void call_test_skinny_cubes() {
+    typedef Cartesian <NT> Kernel;
+    typedef typename Kernel::Point Point;
+    typedef HPolytope <Point> Hpolytope;
+    Hpolytope P;
+
+    std::cout << "\n--- Testing rounding of H-skinny_cube5" << std::endl;
+    P = generate_skinny_cube<Hpolytope>(5);
+    rounding_test(P, 0, 3070.64, 3188.25, 3140.6, 3200.0);
+
+    std::cout << "\n--- Testing rounding of H-skinny_cube10" << std::endl;
+
+    P = generate_skinny_cube<Hpolytope>(10);
+    rounding_test(P, 0, 122550, 108426, 105003.0, 102400.0);
+
+    std::cout << "\n--- Testing rounding of H-skinny_cube20" << std::endl;
+    P = generate_skinny_cube<Hpolytope>(20);
+    rounding_test(P, 0,
+                  8.26497 * std::pow(10,7),
+                  8.94948e+07,
+                  1.09218e+08,
+                  104857600.0);
+}
+
+
+TEST_CASE("round_skinny_cube") {
+    call_test_skinny_cubes<double>();
+}

test/new_rounding_test.cpp

@@ -110,7 +110,7 @@ void call_test_skinny_cubes() {
     P = generate_skinny_cube<Hpolytope>(20);
     rounding_test(P, 0,
                   8.26497 * std::pow(10,7),
-                  8.94948+07,
+                  8.94948e+07,
                   1.09218e+08,
                   104857600.0);
 }