diff --git a/include/preprocess/analytic_center_linear_ineq.h b/include/preprocess/analytic_center_linear_ineq.h deleted file mode 100644 index 80d731f6c..000000000 --- a/include/preprocess/analytic_center_linear_ineq.h +++ /dev/null @@ -1,134 +0,0 @@ -// VolEsti (volume computation and sampling library) - -// Copyright (c) 2024 Vissarion Fisikopoulos -// Copyright (c) 2024 Apostolos Chalkis -// Copyright (c) 2024 Elias Tsigaridas - -// Licensed under GNU LGPL.3, see LICENCE file - - -#ifndef ANALYTIC_CENTER_H -#define ANALYTIC_CENTER_H - -#include - -#include "preprocess/max_inscribed_ball.hpp" -#include "preprocess/feasible_point.hpp" -#include "preprocess/mat_computational_operators.h" - -template -NT get_max_step(VT const& Ad, VT const& b_Ax) -{ - const int m = Ad.size(); - NT max_element = std::numeric_limits::lowest(), max_element_temp; - for (int i = 0; i < m; i++) { - max_element_temp = Ad.coeff(i) / b_Ax.coeff(i); - if (max_element_temp > max_element) { - max_element = max_element_temp; - } - } - - return NT(1) / max_element; -} - -template -void get_hessian_grad_logbarrier(MT const& A, MT const& A_trans, VT const& b, - VT const& x, VT const& Ax, MT &H, VT &grad, VT &b_Ax) -{ - b_Ax.noalias() = b - Ax; - VT s = b_Ax.cwiseInverse(); - VT s_sq = s.cwiseProduct(s); - // Gradient of the log-barrier function - grad.noalias() = A_trans * s; - // Hessian of the log-barrier function - update_Atrans_Diag_A(H, A_trans, A, s_sq.asDiagonal()); -} - -/* - This implementation computes the analytic center of a polytope given - as a set of linear inequalities P = {x | Ax <= b}. The analytic center - is the minimizer of the log barrier function i.e., the optimal solution - of the following optimization problem (Convex Optimization, Boyd and Vandenberghe, Section 8.5.3), - - \min -\sum \log(b_i - a_i^Tx), where a_i is the i-th row of A. - - The function solves the problem by using the Newton method. - - Input: (i) Matrix A, vector b such that the polytope P = {x | Ax<=b} - (ii) The number of maximum iterations, max_iters - (iii) Tolerance parameter grad_err_tol to bound the L2-norm of the gradient - (iv) Tolerance parameter rel_pos_err_tol to check the relative progress in each iteration - - Output: (i) The Hessian computed on the analytic center - (ii) The analytic center of the polytope - (iii) A boolean variable that declares convergence - - Note: Using MT as to deal with both dense and sparse matrices, MT_dense will be the type of result matrix -*/ -template -std::tuple analytic_center_linear_ineq(MT const& A, VT const& b, VT const& x0, - unsigned int const max_iters = 500, - NT const grad_err_tol = 1e-08, - NT const rel_pos_err_tol = 1e-12) -{ - // Initialization - VT x = x0; - VT Ax = A * x; - const int n = A.cols(), m = A.rows(); - MT H(n, n), A_trans = A.transpose(); - VT grad(n), d(n), Ad(m), b_Ax(m), step_d(n), x_prev; - NT grad_err, rel_pos_err, rel_pos_err_temp, step; - unsigned int iter = 0; - bool converged = false; - const NT tol_bnd = NT(0.01); - - auto llt = initialize_chol(A_trans, A); - get_hessian_grad_logbarrier(A, A_trans, b, x, Ax, H, grad, b_Ax); - - do { - iter++; - // Compute the direction - d.noalias() = - solve_vec(llt, H, grad); - Ad.noalias() = A * d; - // Compute the step length - step = std::min((NT(1) - tol_bnd) * get_max_step(Ad, b_Ax), NT(1)); - step_d.noalias() = step*d; - x_prev = x; - x += step_d; - Ax.noalias() += step*Ad; - - // Compute the max_i\{ |step*d_i| ./ |x_i| \} - rel_pos_err = std::numeric_limits::lowest(); - for (int i = 0; i < n; i++) - { - rel_pos_err_temp = std::abs(step_d.coeff(i) / x_prev.coeff(i)); - if (rel_pos_err_temp > rel_pos_err) - { - rel_pos_err = rel_pos_err_temp; - } - } - - get_hessian_grad_logbarrier(A, A_trans, b, x, Ax, H, grad, b_Ax); - grad_err = grad.norm(); - - if (iter >= max_iters || grad_err <= grad_err_tol || rel_pos_err <= rel_pos_err_tol) - { - converged = true; - break; - } - } while (true); - - return std::make_tuple(MT_dense(H), x, converged); -} - -template -std::tuple analytic_center_linear_ineq(MT const& A, VT const& b, - unsigned int const max_iters = 500, - NT const grad_err_tol = 1e-08, - NT const rel_pos_err_tol = 1e-12) -{ - VT x0 = compute_feasible_point(A, b); - return analytic_center_linear_ineq(A, b, x0, max_iters, grad_err_tol, rel_pos_err_tol); -} - -#endif diff --git a/include/preprocess/barrier_center_ellipsoid.hpp b/include/preprocess/barrier_center_ellipsoid.hpp new file mode 100644 index 000000000..0ab8c0457 --- /dev/null +++ b/include/preprocess/barrier_center_ellipsoid.hpp @@ -0,0 +1,115 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2024 Vissarion Fisikopoulos +// Copyright (c) 2024 Apostolos Chalkis +// Copyright (c) 2024 Elias Tsigaridas + +// Licensed under GNU LGPL.3, see LICENCE file + + +#ifndef BARRIER_CENTER_ELLIPSOID_HPP +#define BARRIER_CENTER_ELLIPSOID_HPP + +#include + +#include "preprocess/max_inscribed_ball.hpp" +#include "preprocess/feasible_point.hpp" +#include "preprocess/rounding_util_functions.hpp" + +/* + This implementation computes the analytic or the volumetric center of a polytope given + as a set of linear inequalities P = {x | Ax <= b}. The analytic center is the tminimizer + of the log barrier function, i.e., the optimal solution + of the following optimization problem (Convex Optimization, Boyd and Vandenberghe, Section 8.5.3), + + \min -\sum \log(b_i - a_i^Tx), where a_i is the i-th row of A. + + 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 + + The function solves the problems by using the Newton method. + + Input: (i) Matrix A, vector b such that the polytope P = {x | Ax<=b} + (ii) The number of maximum iterations, max_iters + (iii) Tolerance parameter grad_err_tol to bound the L2-norm of the gradient + (iv) Tolerance parameter rel_pos_err_tol to check the relative progress in each iteration + + Output: (i) The Hessian of the barrier function + (ii) The analytic/volumetric center of the polytope + (iii) A boolean variable that declares convergence + + Note: Using MT as to deal with both dense and sparse matrices, MT_dense will be the type of result matrix +*/ +template +std::tuple barrier_center_ellipsoid_linear_ineq(MT const& A, VT const& b, VT const& x0, + unsigned int const max_iters = 500, + NT const grad_err_tol = 1e-08, + NT const rel_pos_err_tol = 1e-12) +{ + // Initialization + VT x = x0; + VT Ax = A * x; + const int n = A.cols(), m = A.rows(); + MT H(n, n), A_trans = A.transpose(); + VT grad(n), d(n), Ad(m), b_Ax(m), step_d(n), x_prev; + NT grad_err, rel_pos_err, rel_pos_err_temp, step, obj_val, obj_val_prev; + unsigned int iter = 0; + bool converged = false; + const NT tol_bnd = NT(0.01), tol_obj = NT(1e-06); + + auto [step_iter, max_step_multiplier] = init_step(); + auto llt = initialize_chol(A_trans, A); + get_barrier_hessian_grad(A, A_trans, b, x, Ax, llt, + H, grad, b_Ax, obj_val); + do { + iter++; + // Compute the direction + d.noalias() = - solve_vec(llt, H, grad); + Ad.noalias() = A * d; + // Compute the step length + step = std::min(max_step_multiplier * get_max_step(Ad, b_Ax), step_iter); + step_d.noalias() = step*d; + x_prev = x; + x += step_d; + Ax.noalias() += step*Ad; + + // Compute the max_i\{ |step*d_i| ./ |x_i| \} + rel_pos_err = std::numeric_limits::lowest(); + for (int i = 0; i < n; i++) + { + rel_pos_err_temp = std::abs(step_d.coeff(i) / x_prev.coeff(i)); + if (rel_pos_err_temp > rel_pos_err) + { + rel_pos_err = rel_pos_err_temp; + } + } + + obj_val_prev = obj_val; + get_barrier_hessian_grad(A, A_trans, b, x, Ax, llt, + H, grad, b_Ax, obj_val); + grad_err = grad.norm(); + + if (iter >= max_iters || grad_err <= grad_err_tol || rel_pos_err <= rel_pos_err_tol) + { + converged = true; + break; + } + get_step_next_iteration(obj_val_prev, obj_val, tol_obj, step_iter); + } while (true); + + return std::make_tuple(MT_dense(H), x, converged); +} + +template +std::tuple barrier_center_ellipsoid_linear_ineq(MT const& A, VT const& b, + unsigned int const max_iters = 500, + NT const grad_err_tol = 1e-08, + NT const rel_pos_err_tol = 1e-12) +{ + VT x0 = compute_feasible_point(A, b); + return barrier_center_ellipsoid_linear_ineq(A, b, x0, max_iters, grad_err_tol, rel_pos_err_tol); +} + +#endif // BARRIER_CENTER_ELLIPSOID_HPP diff --git a/include/preprocess/inscribed_ellipsoid_rounding.hpp b/include/preprocess/inscribed_ellipsoid_rounding.hpp index dd20d533d..37468b9d9 100644 --- a/include/preprocess/inscribed_ellipsoid_rounding.hpp +++ b/include/preprocess/inscribed_ellipsoid_rounding.hpp @@ -12,14 +12,9 @@ #define INSCRIBED_ELLIPSOID_ROUNDING_HPP #include "preprocess/max_inscribed_ellipsoid.hpp" -#include "preprocess/analytic_center_linear_ineq.h" +#include "preprocess/barrier_center_ellipsoid.hpp" #include "preprocess/feasible_point.hpp" -enum EllipsoidType -{ - MAX_ELLIPSOID = 1, - LOG_BARRIER = 2 -}; template inline static std::tuple @@ -30,9 +25,10 @@ compute_inscribed_ellipsoid(Custom_MT A, VT b, VT const& x0, if constexpr (ellipsoid_type == EllipsoidType::MAX_ELLIPSOID) { return max_inscribed_ellipsoid(A, b, x0, maxiter, tol, reg); - } else if constexpr (ellipsoid_type == EllipsoidType::LOG_BARRIER) + } else if constexpr (ellipsoid_type == EllipsoidType::LOG_BARRIER || + ellipsoid_type == EllipsoidType::VOLUMETRIC_BARRIER) { - return analytic_center_linear_ineq(A, b, x0); + return barrier_center_ellipsoid_linear_ineq(A, b, x0); } else { std::runtime_error("Unknown rounding method."); diff --git a/include/preprocess/max_inscribed_ball.hpp b/include/preprocess/max_inscribed_ball.hpp index 573813dd7..5f95970ba 100644 --- a/include/preprocess/max_inscribed_ball.hpp +++ b/include/preprocess/max_inscribed_ball.hpp @@ -11,7 +11,7 @@ #ifndef MAX_INSCRIBED_BALL_HPP #define MAX_INSCRIBED_BALL_HPP -#include "preprocess/mat_computational_operators.h" +#include "preprocess/rounding_util_functions.hpp" /* This implmentation computes the largest inscribed ball in a given convex polytope P. diff --git a/include/preprocess/max_inscribed_ellipsoid.hpp b/include/preprocess/max_inscribed_ellipsoid.hpp index f44f60ff0..cc65006b9 100644 --- a/include/preprocess/max_inscribed_ellipsoid.hpp +++ b/include/preprocess/max_inscribed_ellipsoid.hpp @@ -15,7 +15,7 @@ #include #include -#include "preprocess/mat_computational_operators.h" +#include "preprocess/rounding_util_functions.hpp" /* diff --git a/include/preprocess/mat_computational_operators.h b/include/preprocess/rounding_util_functions.hpp similarity index 76% rename from include/preprocess/mat_computational_operators.h rename to include/preprocess/rounding_util_functions.hpp index af200d221..6b94ace8b 100644 --- a/include/preprocess/mat_computational_operators.h +++ b/include/preprocess/rounding_util_functions.hpp @@ -1,14 +1,15 @@ // VolEsti (volume computation and sampling library) -// Copyright (c) 2024 Vissarion Fisikopoulos -// Copyright (c) 2024 Apostolos Chalkis -// Copyright (c) 2024 Elias Tsigaridas +// Copyright (c) 2012-2024 Vissarion Fisikopoulos +// Copyright (c) 2018-2024 Apostolos Chalkis + +//Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program. // Licensed under GNU LGPL.3, see LICENCE file -#ifndef MAT_COMPUTATIONAL_OPERATORS_H -#define MAT_COMPUTATIONAL_OPERATORS_H +#ifndef ROUNDING_UTIL_FUNCTIONS_HPP +#define ROUNDING_UTIL_FUNCTIONS_HPP #include @@ -17,9 +18,35 @@ #include "Spectra/include/Spectra/MatOp/SparseSymMatProd.h" +enum EllipsoidType +{ + MAX_ELLIPSOID = 1, + LOG_BARRIER = 2, + VOLUMETRIC_BARRIER = 3 +}; + +template +struct AssertBarrierFalseType : std::false_type {}; + template struct AssertFalseType : std::false_type {}; +template +auto get_max_step(VT const& Ad, VT const& b_Ax) +{ + using NT = typename VT::Scalar; + const int m = Ad.size(); + NT max_element = std::numeric_limits::lowest(), max_element_temp; + for (int i = 0; i < m; i++) { + max_element_temp = Ad.coeff(i) / b_Ax.coeff(i); + if (max_element_temp > max_element) { + max_element = max_element_temp; + } + } + + return NT(1) / max_element; +} + template inline static auto initialize_chol(MT const& mat) @@ -312,5 +339,64 @@ update_Bmat(MT &B, VT const& AtDe, VT const& d, } } +template +std::tuple init_step() +{ + if constexpr (BarrierType == EllipsoidType::LOG_BARRIER) + { + return {NT(1), NT(0.99)}; + } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER) + { + return {NT(0.5), NT(0.4)}; + } else { + static_assert(AssertBarrierFalseType::value, + "Barrier type is not supported."); + } +} + +template +void get_barrier_hessian_grad(MT const& A, MT const& A_trans, VT const& b, + VT const& x, VT const& Ax, llt_type const& llt, + MT &H, VT &grad, VT &b_Ax, NT &obj_val) +{ + b_Ax.noalias() = b - Ax; + VT s = b_Ax.cwiseInverse(); + VT s_sq = s.cwiseProduct(s); + // Hessian of the log-barrier function + update_Atrans_Diag_A(H, A_trans, A, s_sq.asDiagonal()); + 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 { + static_assert(AssertBarrierFalseType::value, + "Barrier type is not supported."); + } +} + +template +void get_step_next_iteration(NT const obj_val_prev, NT const obj_val, + NT const tol_obj, NT &step_iter) +{ + if constexpr (BarrierType == EllipsoidType::LOG_BARRIER) + { + step_iter = NT(1); + } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER) + { + step_iter *= (obj_val_prev <= obj_val - tol_obj) ? NT(0.9) : NT(0.999); + } else { + static_assert(AssertBarrierFalseType::value, + "Barrier type is not supported."); + } +} -#endif // MAT_COMPUTATIONAL_OPERATORS_H +#endif // ROUNDING_UTIL_FUNCTIONS_HPP \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index cdda6087c..b524403ce 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -300,6 +300,8 @@ add_test(NAME test_round_log_barrier_test COMMAND rounding_test -tc=round_log_barrier_test) 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) @@ -363,7 +365,9 @@ add_test(NAME test_analytic_center COMMAND test_internal_points -tc=test_analytic_center) 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) + set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING") diff --git a/test/rounding_test.cpp b/test/rounding_test.cpp index 7e15c99a5..3817cebb3 100644 --- a/test/rounding_test.cpp +++ b/test/rounding_test.cpp @@ -198,6 +198,36 @@ void rounding_log_barrier_test(Polytope &HP, test_values(volume, expectedBilliard, exact); } +template +void rounding_volumetric_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); + std::pair InnerBall = HP.ComputeInnerBall(); + std::tuple res = inscribed_ellipsoid_rounding(HP, InnerBall.first); + // 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, @@ -284,6 +314,18 @@ void call_test_log_barrier() { rounding_log_barrier_test(P, 0, 3070.64, 3188.25, 3262.77, 3200.0); } +template +void call_test_volumetric_barrier() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "\n--- Testing volumetric barrier rounding of H-skinny_cube5" << std::endl; + P = generate_skinny_cube(5); + rounding_volumetric_barrier_test(P, 0, 3070.64, 3188.25, 3262.77, 3200.0); +} + template void call_test_svd() { @@ -314,6 +356,10 @@ TEST_CASE("round_log_barrier_test") { call_test_log_barrier(); } +TEST_CASE("round_volumetric_barrier_test") { + call_test_volumetric_barrier(); +} + TEST_CASE("round_svd") { call_test_svd(); } diff --git a/test/test_internal_points.cpp b/test/test_internal_points.cpp index 8d4df0240..50965128c 100644 --- a/test/test_internal_points.cpp +++ b/test/test_internal_points.cpp @@ -17,7 +17,7 @@ #include "convex_bodies/hpolytope.h" #include "preprocess/max_inscribed_ball.hpp" -#include "preprocess/analytic_center_linear_ineq.h" +#include "preprocess/barrier_center_ellipsoid.hpp" #include "generators/known_polytope_generators.h" #include "generators/h_polytopes_generator.h" @@ -129,10 +129,10 @@ void call_test_analytic_center() { P = skinny_random_hpoly(3, 15, pre_rounding, max_min_eig_ratio, 127); P.normalize(); - auto [Hessian, analytic_center, converged] = analytic_center_linear_ineq(P.get_mat(), P.get_vec()); - SpMT Asp = P.get_mat().sparseView(); + auto [Hessian, analytic_center, converged] = barrier_center_ellipsoid_linear_ineq(P.get_mat(), P.get_vec()); - auto [Hessian_sp, analytic_center2, converged2] = analytic_center_linear_ineq(Asp, P.get_vec()); + SpMT Asp = P.get_mat().sparseView(); + auto [Hessian_sp, analytic_center2, converged2] = barrier_center_ellipsoid_linear_ineq(Asp, P.get_vec()); CHECK(P.is_in(Point(analytic_center)) == -1); CHECK(converged); @@ -149,6 +149,41 @@ void call_test_analytic_center() { CHECK((Hessian - Hessian_sp).norm() < 1e-12); } +template +void call_test_volumetric_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 volumetric 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, volumetric_center, converged] = barrier_center_ellipsoid_linear_ineq(P.get_mat(), P.get_vec()); + SpMT Asp = P.get_mat().sparseView(); + auto [Hessian_sp, volumetric_center2, converged2] = barrier_center_ellipsoid_linear_ineq(Asp, P.get_vec()); + CHECK(P.is_in(Point(volumetric_center)) == -1); + CHECK(converged); + CHECK(std::abs(volumetric_center(0) + 1.49031) < 1e-04); + CHECK(std::abs(volumetric_center(1) + 1.51709) < 1e-04); + CHECK(std::abs(volumetric_center(2) - 2.49381) < 1e-04); + + CHECK(P.is_in(Point(volumetric_center2)) == -1); + CHECK(converged2); + CHECK(std::abs(volumetric_center(0) - volumetric_center2(0)) < 1e-12); + CHECK(std::abs(volumetric_center(1) - volumetric_center2(1)) < 1e-12); + CHECK(std::abs(volumetric_center(2) - volumetric_center2(2)) < 1e-12); + + CHECK((Hessian - Hessian_sp).norm() < 1e-12); +} + TEST_CASE("test_max_ball") { call_test_max_ball(); } @@ -161,6 +196,10 @@ TEST_CASE("test_analytic_center") { call_test_analytic_center(); } +TEST_CASE("test_volumetric_center") { + call_test_volumetric_center(); +} + TEST_CASE("test_max_ball_sparse") { call_test_max_ball_sparse(); }