diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index bee0fd369..9f4083bd7 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -116,13 +116,13 @@ class HPolytope { if (_inner_ball.second <= NT(0)) { - NT const tol = 0.00000001; - std::tuple inner_ball = max_inscribed_ball(A, b, 150, tol); + NT const tol = 1e-08; + std::tuple inner_ball = max_inscribed_ball(A, b, 5000, tol); // check if the solution is feasible - if (is_in(Point(std::get<0>(inner_ball))) == 0 || std::get<1>(inner_ball) < NT(0) || + if (is_in(Point(std::get<0>(inner_ball))) == 0 || std::get<1>(inner_ball) < tol/2.0 || std::isnan(std::get<1>(inner_ball)) || std::isinf(std::get<1>(inner_ball)) || - !std::get<2>(inner_ball) || is_inner_point_nan_inf(std::get<0>(inner_ball))) + is_inner_point_nan_inf(std::get<0>(inner_ball))) { _inner_ball.second = -1.0; } else diff --git a/include/generators/h_polytopes_generator.h b/include/generators/h_polytopes_generator.h index dae44bdca..e9638b62f 100644 --- a/include/generators/h_polytopes_generator.h +++ b/include/generators/h_polytopes_generator.h @@ -9,7 +9,11 @@ #define H_POLYTOPES_GEN_H #include +#include +#include +#include +#include "preprocess/max_inscribed_ellipsoid_rounding.hpp" #ifndef isnan using std::isnan; @@ -19,17 +23,17 @@ /// @tparam Polytope Type of returned polytope /// @tparam RNGType RNGType Type template -Polytope random_hpoly(unsigned int dim, unsigned int m, double seed = std::numeric_limits::signaling_NaN()) { +Polytope random_hpoly(unsigned int dim, unsigned int m, int seed = std::numeric_limits::signaling_NaN()) { typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; typedef typename Polytope::NT NT; typedef typename Polytope::PointType Point; - unsigned rng_seed = std::chrono::system_clock::now().time_since_epoch().count(); + int rng_seed = std::chrono::system_clock::now().time_since_epoch().count(); RNGType rng(rng_seed); if (!isnan(seed)) { - unsigned rng_seed = seed; + int rng_seed = seed; rng.seed(rng_seed); } @@ -58,4 +62,76 @@ Polytope random_hpoly(unsigned int dim, unsigned int m, double seed = std::numer return Polytope(dim, A, b); } +/// This function generates a transformation that maps the unit ball to a skinny ellipsoid +/// with given ratio between the lengths of its maximum and minimum axis +template +MT get_skinny_transformation(const int d, NT const eig_ratio, int const seed) +{ + boost::normal_distribution<> gdist(0, 1); + RNGType rng(seed); + + MT W(d, d); + for (int i = 0; i < d; i++) { + for (int j = 0; j < d; j++) { + W(i, j) = gdist(rng); + } + } + + Eigen::HouseholderQR qr(W); + MT Q = qr.householderQ(); + + VT diag(d); + const NT eig_min = NT(1), eig_max = eig_ratio; + diag(0) = eig_min; + diag(d-1) = eig_max; + boost::random::uniform_real_distribution udist(NT(0), NT(1)); + NT rand; + for (int i = 1; i < d-1; i++) { + rand = udist(rng); + diag(i) = rand * eig_max + (NT(1)-rand) * eig_min; + } + std::sort(diag.begin(), diag.end()); + MT cov = Q * diag.asDiagonal() * Q.transpose(); + + return cov; +} + +/// This function generates a skinny random H-polytope of given dimension and number of hyperplanes $m$ +/// @tparam Polytope Type of returned polytope +/// @tparam NT Number type +/// @tparam RNGType RNGType Type +template +Polytope skinny_random_hpoly(unsigned int dim, unsigned int m, const bool pre_rounding = false, + const NT eig_ratio = NT(1000.0), int seed = std::numeric_limits::signaling_NaN()) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::PointType Point; + + int rng_seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(rng_seed); + if (!isnan(seed)) { + int rng_seed = seed; + rng.seed(rng_seed); + } + + Polytope P = random_hpoly(dim, m, seed); + + // rounding the polytope before applying the skinny transformation + if (pre_rounding) { + Point x0(dim); + // run only one iteration + max_inscribed_ellipsoid_rounding(P, x0, 1); + } + + MT cov = get_skinny_transformation(dim, eig_ratio, seed); + Eigen::LLT lltOfA(cov); + MT L = lltOfA.matrixL(); + P.linear_transformIt(L.inverse()); + + return P; +} + + + #endif diff --git a/include/preprocess/analytic_center_linear_ineq.h b/include/preprocess/analytic_center_linear_ineq.h new file mode 100644 index 000000000..c7918635a --- /dev/null +++ b/include/preprocess/analytic_center_linear_ineq.h @@ -0,0 +1,134 @@ +// 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 "max_inscribed_ball.hpp" + +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) +{ + const int m = A.rows(); + VT s(m); + + b_Ax.noalias() = b - Ax; + NT *s_data = s.data(); + for (int i = 0; i < m; i++) + { + *s_data = NT(1) / b_Ax.coeff(i); + s_data++; + } + + VT s_sq = s.cwiseProduct(s); + // Gradient of the log-barrier function + grad.noalias() = A_trans * s; + // Hessian of the log-barrier function + H.noalias() = A_trans * s_sq.asDiagonal() * A; +} + +/* + 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 analytic center of the polytope + (ii) A boolean variable that declares convergence +*/ +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 x; + bool feasibility_only = true, converged; + // Compute a feasible point + std::tie(x, std::ignore, converged) = max_inscribed_ball(A, b, max_iters, 1e-08, feasibility_only); + VT Ax = A * x; + if (!converged || (Ax.array() > b.array()).any()) + { + std::runtime_error("The computation of the analytic center failed."); + } + // Initialization + 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; + converged = false; + const NT tol_bnd = NT(0.01); + + get_hessian_grad_logbarrier(A, A_trans, b, x, Ax, H, grad, b_Ax); + + do { + iter++; + // Compute the direction + d.noalias() = - H.llt().solve(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(x, converged); +} + +#endif diff --git a/include/preprocess/max_inscribed_ball.hpp b/include/preprocess/max_inscribed_ball.hpp index 6fe70e9a0..672be24e4 100644 --- a/include/preprocess/max_inscribed_ball.hpp +++ b/include/preprocess/max_inscribed_ball.hpp @@ -27,7 +27,7 @@ */ template -void calcstep(MT const& A, MT const& A_trans, MT const& R, VT &s, +void calcstep(MT const& A, MT const& A_trans, Eigen::LLT const& lltOfB, VT &s, VT &y, VT &r1, VT const& r2, NT const& r3, VT &r4, VT &dx, VT &ds, NT &dt, VT &dy, VT &tmp, VT &rhs) { @@ -41,7 +41,8 @@ void calcstep(MT const& A, MT const& A_trans, MT const& R, VT &s, rhs.block(0,0,n,1).noalias() = r2 + A_trans * tmp; rhs(n) = r3 + tmp.sum(); - VT dxdt = R.colPivHouseholderQr().solve(R.transpose().colPivHouseholderQr().solve(rhs)); + + VT dxdt = lltOfB.solve(rhs); dx = dxdt.block(0,0,n,1); dt = dxdt(n); @@ -57,10 +58,12 @@ void calcstep(MT const& A, MT const& A_trans, MT const& R, VT &s, template -std::tuple max_inscribed_ball(MT const& A, VT const& b, unsigned int maxiter, NT tol) +std::tuple max_inscribed_ball(MT const& A, VT const& b, + unsigned int maxiter, NT tol, + const bool feasibility_only = false) { int m = A.rows(), n = A.cols(); - bool converge; + bool converge = false; NT bnrm = b.norm(); VT o_m = VT::Zero(m), o_n = VT::Zero(n), e_m = VT::Ones(m); @@ -81,7 +84,7 @@ std::tuple max_inscribed_ball(MT const& A, VT const& b, unsigned NT const tau0 = 0.995, power_num = 5.0 * std::pow(10.0, 15.0); NT *vec_iter1, *vec_iter2, *vec_iter3, *vec_iter4; - MT B(n + 1, n + 1), AtD(n, m), R(n + 1, n + 1), + MT B(n + 1, n + 1), AtD(n, m), eEye = std::pow(10.0, -14.0) * MT::Identity(n + 1, n + 1), A_trans = A.transpose(); @@ -105,16 +108,17 @@ std::tuple max_inscribed_ball(MT const& A, VT const& b, unsigned total_err = std::max(total_err, rgap); // progress output & check stopping - if (total_err < tol || ( t > 0 && ( (std::abs(t - t_prev) <= tol * t && std::abs(t - t_prev) <= tol * t_prev) - || (t_prev >= (1.0 - tol) * t && i > 0) - || (t <= (1.0 - tol) * t_prev && i > 0) ) ) ) + if ( (total_err < tol && t > 0) || + ( t > 0 && ( (std::abs(t - t_prev) <= tol * std::min(std::abs(t), std::abs(t_prev)) || + std::abs(t - t_prev) <= tol) && i > 10) ) || + (feasibility_only && t > tol/2.0 && i > 0) ) { //converged converge = true; break; } - if (dt > 1000.0 * bnrm || t > 1000000.0 * bnrm) + if ((dt > 10000.0 * bnrm || t > 10000000.0 * bnrm) && i > 20) { //unbounded converge = false; @@ -127,11 +131,11 @@ std::tuple max_inscribed_ball(MT const& A, VT const& b, unsigned vec_iter2 = y.data(); for (int j = 0; j < m; ++j) { *vec_iter1 = std::min(power_num, (*vec_iter2) / (*vec_iter3)); - AtD.col(j).noalias() = A_trans.col(j) * (*vec_iter1); vec_iter1++; vec_iter3++; vec_iter2++; } + AtD.noalias() = A_trans*d.asDiagonal(); AtDe.noalias() = AtD * e_m; B.block(0, 0, n, n).noalias() = AtD * A; @@ -141,11 +145,10 @@ std::tuple max_inscribed_ball(MT const& A, VT const& b, unsigned B.noalias() += eEye; // Cholesky decomposition - Eigen::LLT lltOfB(B); - R = lltOfB.matrixL().transpose(); + Eigen::LLT lltOfB(B); // predictor step & length - calcstep(A, A_trans, R, s, y, r1, r2, r3, r4, dx, ds, dt, dy, tmp, rhs); + calcstep(A, A_trans, lltOfB, s, y, r1, r2, r3, r4, dx, ds, dt, dy, tmp, rhs); alphap = -1.0; alphad = -1.0; @@ -172,7 +175,7 @@ std::tuple max_inscribed_ball(MT const& A, VT const& b, unsigned // corrector and combined step & length mu_ds_dy.noalias() = e_m * mu - ds.cwiseProduct(dy); - calcstep(A, A_trans, R, s, y, o_m, o_n, 0.0, mu_ds_dy, dxc, dsc, dtc, dyc, tmp, rhs); + calcstep(A, A_trans, lltOfB, s, y, o_m, o_n, 0.0, mu_ds_dy, dxc, dsc, dtc, dyc, tmp, rhs); dx += dxc; ds += dsc; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index a3f8f03d1..e17848fc6 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -297,6 +297,8 @@ add_test(NAME test_round_max_ellipsoid add_test(NAME test_round_svd COMMAND new_rounding_test -tc=round_svd) + + add_executable (logconcave_sampling_test logconcave_sampling_test.cpp $) add_test(NAME logconcave_sampling_test_hmc COMMAND logconcave_sampling_test -tc=hmc) @@ -348,6 +350,14 @@ add_test(NAME test_new_rdhr_uniform_MT COMMAND matrix_sampling_test -tc=new_rdhr add_test(NAME test_new_billiard_uniform_MT COMMAND matrix_sampling_test -tc=new_billiard_uniform_MT) add_test(NAME test_new_accelerated_billiard_uniform_MT COMMAND matrix_sampling_test -tc=new_accelerated_billiard_uniform_MT) +add_executable (test_internal_points test_internal_points.cpp $) +add_test(NAME test_max_ball + COMMAND test_internal_points -tc=test_max_ball) +add_test(NAME test_feasibility_point + COMMAND test_internal_points -tc=test_feasibility_point) +add_test(NAME test_analytic_center + COMMAND test_internal_points -tc=test_analytic_center) + set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING") #set_target_properties(benchmarks_crhmc @@ -388,3 +398,4 @@ TARGET_LINK_LIBRARIES(logconcave_sampling_test lp_solve ${IFOPT} ${IFOPT_IPOPT} TARGET_LINK_LIBRARIES(crhmc_sampling_test lp_solve ${IFOPT} ${IFOPT_IPOPT} ${PTHREAD} ${GMP} ${MPSOLVE} ${FFTW3} ${MKL_LINK} QD_LIB coverage_config) TARGET_LINK_LIBRARIES(order_polytope lp_solve coverage_config) TARGET_LINK_LIBRARIES(matrix_sampling_test lp_solve ${MKL_LINK} coverage_config) +TARGET_LINK_LIBRARIES(test_internal_points lp_solve ${MKL_LINK} coverage_config) diff --git a/test/test_internal_points.cpp b/test/test_internal_points.cpp new file mode 100644 index 000000000..fa44e9baf --- /dev/null +++ b/test/test_internal_points.cpp @@ -0,0 +1,107 @@ +// 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 + +#include "doctest.h" +#include +#include + +#include + +#include "misc/misc.h" +#include "cartesian_geom/cartesian_kernel.h" +#include "convex_bodies/hpolytope.h" + +#include "preprocess/max_inscribed_ball.hpp" +#include "preprocess/analytic_center_linear_ineq.h" + +#include "generators/known_polytope_generators.h" +#include "generators/h_polytopes_generator.h" + +template +void call_test_max_ball() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef HPolytope Hpolytope; + typedef boost::mt19937 PolyRNGType; + Hpolytope P; + + std::cout << "\n--- Testing Chebychev ball for skinny H-polytope" << std::endl; + bool pre_rounding = true; // round random polytope before applying the skinny transformation + NT max_min_eig_ratio = NT(2000); + P = skinny_random_hpoly(4, 180, pre_rounding, max_min_eig_ratio, 127); + P.normalize(); + std::pair InnerBall = P.ComputeInnerBall(); + + NT tol = 1e-08; + unsigned int maxiter = 500; + auto [center, radius, converged] = max_inscribed_ball(P.get_mat(), P.get_vec(), maxiter, tol); + + CHECK(P.is_in(Point(center)) == -1); + CHECK(std::abs(radius - InnerBall.second) <= 1e-03); + CHECK(converged); +} + +template +void call_test_max_ball_feasibility() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef HPolytope Hpolytope; + typedef boost::mt19937 PolyRNGType; + Hpolytope P; + + std::cout << "\n--- Testing feasibility point for skinny H-polytope" << std::endl; + bool pre_rounding = true; // round random polytope before applying the skinny transformation + NT max_min_eig_ratio = NT(2000); + P = skinny_random_hpoly(50, 500, pre_rounding, max_min_eig_ratio, 127); + P.normalize(); + + bool feasibility_only = true; // compute only a feasible point + NT tol = 1e-08; + unsigned int maxiter = 500; + auto [center, radius, converged] = max_inscribed_ball(P.get_mat(), P.get_vec(), maxiter, tol, feasibility_only); + + CHECK(P.is_in(Point(center)) == -1); + CHECK(converged); +} + +template +void call_test_analytic_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; + Hpolytope P; + + std::cout << "\n--- Testing analytic 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 [analytic_center, converged] = analytic_center_linear_ineq(P.get_mat(), P.get_vec()); + + CHECK(P.is_in(Point(analytic_center)) == -1); + CHECK(converged); + CHECK(std::abs(analytic_center(0) + 4.75912) < 1e-04); + CHECK(std::abs(analytic_center(1) + 4.28762) < 1e-04); + CHECK(std::abs(analytic_center(2) - 7.54156) < 1e-04); +} + +TEST_CASE("test_max_ball") { + call_test_max_ball(); +} + +TEST_CASE("test_feasibility_point") { + call_test_max_ball_feasibility(); +} + +TEST_CASE("test_analytic_center") { + call_test_analytic_center(); +}