Rounding with vaidya barrier (#316)

* implement vaidya minimizer and hessian computation

* implement unit tests and update documentation

* resolve PR commits

include/preprocess/barrier_center_ellipsoid.hpp

@@ -27,7 +27,12 @@
     The volumetric center is the minimizer of the volumetric barrier function, i.e., the optimal
     solution of the following optimization problem,
 
-    \min logdet \nabla^2 f(x), where f(x) the log barrier function 
+    \min logdet \nabla^2 f(x), where f(x) the log barrier function
+
+    The Vaidya center is the minimizer of the Vaidya barrier function, i.e., the optimal
+    solution of the following optimization problem,
+
+    \min logdet \nabla^2 f(x) + d/m f(x), where f(x) the log barrier function.
     
     The function solves the problems by using the Newton method.
 

include/preprocess/inscribed_ellipsoid_rounding.hpp

@@ -26,7 +26,8 @@ compute_inscribed_ellipsoid(Custom_MT A, VT b, VT const& x0,
     {
         return max_inscribed_ellipsoid<MT>(A, b, x0, maxiter, tol, reg);
     } else if constexpr (ellipsoid_type == EllipsoidType::LOG_BARRIER ||
-                         ellipsoid_type == EllipsoidType::VOLUMETRIC_BARRIER)
+                         ellipsoid_type == EllipsoidType::VOLUMETRIC_BARRIER ||
+                         ellipsoid_type == EllipsoidType::VAIDYA_BARRIER)
     {
         return barrier_center_ellipsoid_linear_ineq<MT, ellipsoid_type, NT>(A, b, x0);
     } else

include/preprocess/rounding_util_functions.hpp

@@ -22,7 +22,8 @@ enum EllipsoidType
 {
   MAX_ELLIPSOID = 1,
   LOG_BARRIER = 2,
-  VOLUMETRIC_BARRIER = 3
+  VOLUMETRIC_BARRIER = 3,
+  VAIDYA_BARRIER = 4
 };
 
 template <int T>
@@ -345,7 +346,8 @@ std::tuple<NT, NT> init_step()
   if constexpr (BarrierType == EllipsoidType::LOG_BARRIER)
   {
     return {NT(1), NT(0.99)};
-  } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER)
+  } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER ||
+                       BarrierType == EllipsoidType::VAIDYA_BARRIER)
   {
     return {NT(0.5), NT(0.4)};
   } else {
@@ -362,21 +364,44 @@ void get_barrier_hessian_grad(MT const& A, MT const& A_trans, VT const& b,
   b_Ax.noalias() = b - Ax;
   VT s = b_Ax.cwiseInverse();
   VT s_sq = s.cwiseProduct(s);
+  VT sigma;
   // Hessian of the log-barrier function
   update_Atrans_Diag_A<NT>(H, A_trans, A, s_sq.asDiagonal());
+
+  if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER ||
+                BarrierType == EllipsoidType::VAIDYA_BARRIER)
+  {
+    // Computing sigma(x)_i = (a_i^T H^{-1} a_i) / (b_i - a_i^Tx)^2
+    MT_dense HA = solve_mat(llt, H, A_trans, obj_val);
+    MT_dense aiHai = HA.transpose().cwiseProduct(A);
+    sigma = (aiHai.rowwise().sum()).cwiseProduct(s_sq);
+  }
+
   if constexpr (BarrierType == EllipsoidType::LOG_BARRIER)
   {
     grad.noalias() = A_trans * s;
   } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER)
   {
-    // Computing sigma(x)_i = (a_i^T H^{-1} a_i) / (b_i - a_i^Tx)^2
-    MT_dense HA = solve_mat(llt, H, A_trans, obj_val);
-    MT_dense aiHai = HA.transpose().cwiseProduct(A);
-    VT sigma = (aiHai.rowwise().sum()).cwiseProduct(s_sq);
     // Gradient of the volumetric barrier function
     grad.noalias() = A_trans * (s.cwiseProduct(sigma));
     // Hessian of the volumetric barrier function
     update_Atrans_Diag_A<NT>(H, A_trans, A, s_sq.cwiseProduct(sigma).asDiagonal());
+  } else if constexpr (BarrierType == EllipsoidType::VAIDYA_BARRIER)
+  {
+    const int m = b.size(), d = x.size();
+    NT const d_m = NT(d) / NT(m);
+    // Weighted gradient of the log barrier function
+    grad.noalias() = A_trans * s;
+    grad *= d_m;
+    // Add the gradient of the volumetric function
+    grad.noalias() += A_trans * (s.cwiseProduct(sigma));
+    // Weighted Hessian of the log barrier function
+    H *= d_m;
+    // Add the Hessian of the volumetric function
+    MT Hvol(d, d);
+    update_Atrans_Diag_A<NT>(Hvol, A_trans, A, s_sq.cwiseProduct(sigma).asDiagonal());
+    H += Hvol;
+    obj_val -= s.array().log().sum();
   } else {
     static_assert(AssertBarrierFalseType<BarrierType>::value,
             "Barrier type is not supported.");
@@ -393,6 +418,9 @@ void get_step_next_iteration(NT const obj_val_prev, NT const obj_val,
   } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER)
   {
     step_iter *= (obj_val_prev <= obj_val - tol_obj) ? NT(0.9) : NT(0.999);
+  } else if constexpr (BarrierType == EllipsoidType::VAIDYA_BARRIER)
+  {
+    step_iter *= NT(0.999);
   } else {
     static_assert(AssertBarrierFalseType<BarrierType>::value,
             "Barrier type is not supported.");

test/CMakeLists.txt

@@ -302,6 +302,8 @@ add_test(NAME round_max_ellipsoid_sparse
           COMMAND rounding_test -tc=round_max_ellipsoid_sparse)
 add_test(NAME round_volumetric_barrier_test
           COMMAND rounding_test -tc=round_volumetric_barrier_test)
+add_test(NAME round_vaidya_barrier_test
+          COMMAND rounding_test -tc=round_vaidya_barrier_test)
 
 
 
@@ -367,6 +369,8 @@ add_test(NAME test_max_ball_sparse
           COMMAND test_internal_points -tc=test_max_ball_sparse)
 add_test(NAME test_volumetric_center
           COMMAND test_internal_points -tc=test_volumetric_center)
+add_test(NAME test_vaidya_center
+          COMMAND test_internal_points -tc=test_vaidya_center)
 
 
 

test/rounding_test.cpp

@@ -228,6 +228,37 @@ void rounding_volumetric_barrier_test(Polytope &HP,
     test_values(volume, expectedBilliard, exact);
 }
 
+template <class Polytope>
+void rounding_vaidya_barrier_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);
+
+    VT x0(d);
+    x0 << 0.7566, 0.6374, 0.3981, 0.9248, 0.9828;
+    std::tuple<MT, VT, NT> res = inscribed_ellipsoid_rounding<MT, VT, NT, Polytope, Point, EllipsoidType::VAIDYA_BARRIER>(HP, Point(x0));
+    // Setup the parameters
+    int walk_len = 1;
+    NT e = 0.1;
+
+    // Estimate the volume
+    std::cout << "Number type: " << typeid(NT).name() << std::endl;
+
+    NT volume = std::get<2>(res) * volume_cooling_balls<BilliardWalk, RNGType>(HP, e, walk_len).second;
+    test_values(volume, expectedBilliard, exact);
+}
 
 template <class Polytope>
 void rounding_svd_test(Polytope &HP,
@@ -326,6 +357,17 @@ void call_test_volumetric_barrier() {
     rounding_volumetric_barrier_test(P, 0, 3070.64, 3188.25, 3262.77, 3200.0);
 }
 
+template <typename NT>
+void call_test_vaidya_barrier() {
+    typedef Cartesian <NT> Kernel;
+    typedef typename Kernel::Point Point;
+    typedef HPolytope <Point> Hpolytope;
+    Hpolytope P;
+
+    std::cout << "\n--- Testing vaidya barrier rounding of H-skinny_cube5" << std::endl;
+    P = generate_skinny_cube<Hpolytope>(5);
+    rounding_vaidya_barrier_test(P, 0, 3070.64, 3188.25, 3262.77, 3200.0);
+}
 
 template <typename NT>
 void call_test_svd() {
@@ -360,6 +402,10 @@ TEST_CASE("round_volumetric_barrier_test") {
     call_test_volumetric_barrier<double>();
 }
 
+TEST_CASE("round_vaidya_barrier_test") {
+    call_test_vaidya_barrier<double>();
+}
+
 TEST_CASE("round_svd") {
     call_test_svd<double>();
 }

test/test_internal_points.cpp

@@ -184,6 +184,42 @@ void call_test_volumetric_center() {
     CHECK((Hessian - Hessian_sp).norm() < 1e-12);
 }
 
+template <typename NT>
+void call_test_vaidya_center() {
+    typedef Cartesian <NT> Kernel;
+    typedef typename Kernel::Point Point;
+    typedef HPolytope <Point> Hpolytope;
+    typedef typename Hpolytope::MT MT;
+    typedef typename Hpolytope::VT VT;
+    typedef boost::mt19937 PolyRNGType;
+    typedef Eigen::SparseMatrix<NT> SpMT;
+    Hpolytope P;
+
+    std::cout << "\n--- Testing vaidya center for skinny H-polytope" << std::endl;
+    bool pre_rounding = true; // round random polytope before applying the skinny transformation 
+    NT max_min_eig_ratio = NT(100);
+    P = skinny_random_hpoly<Hpolytope, NT, PolyRNGType>(3, 15, pre_rounding, max_min_eig_ratio, 127);
+    P.normalize();
+
+    auto [Hessian, vaidya_center, converged] = barrier_center_ellipsoid_linear_ineq<MT, EllipsoidType::VAIDYA_BARRIER, NT>(P.get_mat(), P.get_vec());
+    SpMT Asp = P.get_mat().sparseView();
+    auto [Hessian_sp, vaidya_center2, converged2] = barrier_center_ellipsoid_linear_ineq<MT, EllipsoidType::VAIDYA_BARRIER, NT>(Asp, P.get_vec());
+
+    CHECK(P.is_in(Point(vaidya_center)) == -1);
+    CHECK(converged);
+    CHECK(std::abs(vaidya_center(0) + 2.4076) < 1e-04);
+    CHECK(std::abs(vaidya_center(1) + 2.34072) < 1e-04);
+    CHECK(std::abs(vaidya_center(2) - 3.97138) < 1e-04);
+
+    CHECK(P.is_in(Point(vaidya_center2)) == -1);
+    CHECK(converged2);
+    CHECK(std::abs(vaidya_center(0) - vaidya_center2(0)) < 1e-12);
+    CHECK(std::abs(vaidya_center(1) - vaidya_center2(1)) < 1e-12);
+    CHECK(std::abs(vaidya_center(2) - vaidya_center2(2)) < 1e-12);
+
+    CHECK((Hessian - Hessian_sp).norm() < 1e-12);
+}
+
 TEST_CASE("test_max_ball") {
     call_test_max_ball<double>();
 }
@@ -200,6 +236,10 @@ TEST_CASE("test_volumetric_center") {
     call_test_volumetric_center<double>();
 }
 
+TEST_CASE("test_vaidya_center") {
+    call_test_vaidya_center<double>();
+}
+
 TEST_CASE("test_max_ball_sparse") {
     call_test_max_ball_sparse<double>();
 }