From 992d20321be2b0e64541594a9ed71d16982cb116 Mon Sep 17 00:00:00 2001 From: Apostolos Chalkis Date: Tue, 16 Jul 2024 20:45:27 -0600 Subject: [PATCH] Rounding with vaidya barrier (#316) * implement vaidya minimizer and hessian computation * implement unit tests and update documentation * resolve PR commits --- .../preprocess/barrier_center_ellipsoid.hpp | 7 ++- .../inscribed_ellipsoid_rounding.hpp | 3 +- .../preprocess/rounding_util_functions.hpp | 40 +++++++++++++--- test/CMakeLists.txt | 4 ++ test/rounding_test.cpp | 46 +++++++++++++++++++ test/test_internal_points.cpp | 40 ++++++++++++++++ 6 files changed, 132 insertions(+), 8 deletions(-) diff --git a/include/preprocess/barrier_center_ellipsoid.hpp b/include/preprocess/barrier_center_ellipsoid.hpp index 0ab8c0457..9f440f222 100644 --- a/include/preprocess/barrier_center_ellipsoid.hpp +++ b/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. diff --git a/include/preprocess/inscribed_ellipsoid_rounding.hpp b/include/preprocess/inscribed_ellipsoid_rounding.hpp index 37468b9d9..4b027f61b 100644 --- a/include/preprocess/inscribed_ellipsoid_rounding.hpp +++ b/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(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(A, b, x0); } else diff --git a/include/preprocess/rounding_util_functions.hpp b/include/preprocess/rounding_util_functions.hpp index 6b94ace8b..a189890e3 100644 --- a/include/preprocess/rounding_util_functions.hpp +++ b/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 @@ -345,7 +346,8 @@ std::tuple 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(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(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(Hvol, A_trans, A, s_sq.cwiseProduct(sigma).asDiagonal()); + H += Hvol; + obj_val -= s.array().log().sum(); } else { static_assert(AssertBarrierFalseType::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::value, "Barrier type is not supported."); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b524403ce..615791148 100644 --- a/test/CMakeLists.txt +++ b/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) diff --git a/test/rounding_test.cpp b/test/rounding_test.cpp index 3817cebb3..5beae0415 100644 --- a/test/rounding_test.cpp +++ b/test/rounding_test.cpp @@ -228,6 +228,37 @@ void rounding_volumetric_barrier_test(Polytope &HP, test_values(volume, expectedBilliard, exact); } +template +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 RNGType; + RNGType rng(d); + + VT x0(d); + x0 << 0.7566, 0.6374, 0.3981, 0.9248, 0.9828; + std::tuple res = inscribed_ellipsoid_rounding(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(HP, e, walk_len).second; + test_values(volume, expectedBilliard, exact); +} template 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 +void call_test_vaidya_barrier() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "\n--- Testing vaidya barrier rounding of H-skinny_cube5" << std::endl; + P = generate_skinny_cube(5); + rounding_vaidya_barrier_test(P, 0, 3070.64, 3188.25, 3262.77, 3200.0); +} template void call_test_svd() { @@ -360,6 +402,10 @@ TEST_CASE("round_volumetric_barrier_test") { call_test_volumetric_barrier(); } +TEST_CASE("round_vaidya_barrier_test") { + call_test_vaidya_barrier(); +} + TEST_CASE("round_svd") { call_test_svd(); } diff --git a/test/test_internal_points.cpp b/test/test_internal_points.cpp index 50965128c..d0ee44008 100644 --- a/test/test_internal_points.cpp +++ b/test/test_internal_points.cpp @@ -184,6 +184,42 @@ void call_test_volumetric_center() { CHECK((Hessian - Hessian_sp).norm() < 1e-12); } +template +void call_test_vaidya_center() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef HPolytope Hpolytope; + typedef typename Hpolytope::MT MT; + typedef typename Hpolytope::VT VT; + typedef boost::mt19937 PolyRNGType; + typedef Eigen::SparseMatrix 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(3, 15, pre_rounding, max_min_eig_ratio, 127); + P.normalize(); + + auto [Hessian, vaidya_center, converged] = barrier_center_ellipsoid_linear_ineq(P.get_mat(), P.get_vec()); + SpMT Asp = P.get_mat().sparseView(); + auto [Hessian_sp, vaidya_center2, converged2] = barrier_center_ellipsoid_linear_ineq(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(); } @@ -200,6 +236,10 @@ TEST_CASE("test_volumetric_center") { call_test_volumetric_center(); } +TEST_CASE("test_vaidya_center") { + call_test_vaidya_center(); +} + TEST_CASE("test_max_ball_sparse") { call_test_max_ball_sparse(); }