From 1fb1265efe60d3c8ca55bf0b2e221a4f4d585467 Mon Sep 17 00:00:00 2001 From: "Apostolos Chalkis (TolisChal)" Date: Thu, 27 Jun 2024 02:17:47 -0600 Subject: [PATCH] Generalize the rounding loop and support sparse computations in preprocessing routines (#312) * generalize rounding loop * support sparse cholesky operator * complete sparse support in max_inscribed_ball * complete sparse support in preprocesing * add sparse tests * change main rounding function name * improve explaining comments * resolve PR comments * changing the dates in copyrights * use if constexpr instead of SNIFAE * update the examples to cpp17 * update to cpp17 order polytope example * fix templating in mat_computational_operators * fix templating errors and change header file to mat_computational_operators --- examples/EnvelopeProblemSOS/CMakeLists.txt | 4 +- .../CMakeLists.txt | 1 + .../volesti_lecount.cpp | 2 +- examples/crhmc_prepare/CMakeLists.txt | 1 + examples/crhmc_sampling/CMakeLists.txt | 2 +- examples/ellipsoid-sampling/CMakeLists.txt | 2 +- examples/hpolytope-volume/CMakeLists.txt | 1 + examples/logconcave/CMakeLists.txt | 2 +- examples/mmcs_method/CMakeLists.txt | 2 +- examples/multithread_sampling/CMakeLists.txt | 2 +- .../optimization_spectrahedra/CMakeLists.txt | 1 + examples/order-polytope-basics/CMakeLists.txt | 2 +- .../CMakeLists.txt | 2 +- .../sampler.cpp | 17 +- examples/volume_spectrahedron/CMakeLists.txt | 2 +- examples/vpolytope-volume/CMakeLists.txt | 1 + include/convex_bodies/orderpolytope.h | 6 + include/generators/h_polytopes_generator.h | 4 +- .../preprocess/analytic_center_linear_ineq.h | 62 ++-- .../estimate_L_smooth_parameter.hpp | 4 +- include/preprocess/feasible_point.hpp | 34 ++ .../inscribed_ellipsoid_rounding.hpp | 134 ++++++++ .../preprocess/mat_computational_operators.h | 316 ++++++++++++++++++ include/preprocess/max_inscribed_ball.hpp | 48 ++- .../preprocess/max_inscribed_ellipsoid.hpp | 58 ++-- .../max_inscribed_ellipsoid_rounding.hpp | 94 ------ test/CMakeLists.txt | 18 +- test/max_ellipsoid_rounding_test.cpp | 4 +- ...ew_rounding_test.cpp => rounding_test.cpp} | 117 ++++++- test/test_internal_points.cpp | 63 +++- 30 files changed, 786 insertions(+), 220 deletions(-) create mode 100644 include/preprocess/feasible_point.hpp create mode 100644 include/preprocess/inscribed_ellipsoid_rounding.hpp create mode 100644 include/preprocess/mat_computational_operators.h delete mode 100644 include/preprocess/max_inscribed_ellipsoid_rounding.hpp rename test/{new_rounding_test.cpp => rounding_test.cpp} (61%) diff --git a/examples/EnvelopeProblemSOS/CMakeLists.txt b/examples/EnvelopeProblemSOS/CMakeLists.txt index e8f5b1a61..6d4e2a2a4 100644 --- a/examples/EnvelopeProblemSOS/CMakeLists.txt +++ b/examples/EnvelopeProblemSOS/CMakeLists.txt @@ -1,7 +1,7 @@ cmake_minimum_required(VERSION 3.15) project(EnvelopeProblem) -set(CMAKE_CXX_STANDARD 14) +set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_FLAGS_DEBUG_CUSTOM "-O0 -fno-omit-frame-pointer -mno-omit-leaf-frame-pointer -DIPM_USE_DOUBLE -DIPM_DOUBLE=double") set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${CMAKE_CXX_FLAGS_DEBUG_CUSTOM}") @@ -83,4 +83,4 @@ else () target_link_directories(EnvelopeProblem PRIVATE ${OPENMP_LIBRARIES}) endif () target_link_libraries(EnvelopeProblem Python2::Python Python2::NumPy) -endif () \ No newline at end of file +endif () diff --git a/examples/count-linear-extensions-using-hpolytope/CMakeLists.txt b/examples/count-linear-extensions-using-hpolytope/CMakeLists.txt index f3e66b634..49f295303 100644 --- a/examples/count-linear-extensions-using-hpolytope/CMakeLists.txt +++ b/examples/count-linear-extensions-using-hpolytope/CMakeLists.txt @@ -98,6 +98,7 @@ else () add_compile_definitions("EIGEN_NO_DEBUG") endif () + add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") diff --git a/examples/count-linear-extensions-using-hpolytope/volesti_lecount.cpp b/examples/count-linear-extensions-using-hpolytope/volesti_lecount.cpp index 9b0e99822..eb7133155 100644 --- a/examples/count-linear-extensions-using-hpolytope/volesti_lecount.cpp +++ b/examples/count-linear-extensions-using-hpolytope/volesti_lecount.cpp @@ -10,7 +10,7 @@ #include "volume_cooling_gaussians.hpp" #include "volume_cooling_balls.hpp" -#include "preprocess/max_inscribed_ellipsoid_rounding.hpp" +#include "preprocess/inscribed_ellipsoid_rounding.hpp" #include "preprocess/min_sampling_covering_ellipsoid_rounding.hpp" #include "preprocess/svd_rounding.hpp" diff --git a/examples/crhmc_prepare/CMakeLists.txt b/examples/crhmc_prepare/CMakeLists.txt index 3575c4a46..9d99477c2 100644 --- a/examples/crhmc_prepare/CMakeLists.txt +++ b/examples/crhmc_prepare/CMakeLists.txt @@ -105,6 +105,7 @@ endif () #add_definitions(${CMAKE_CXX_FLAGS} "-g") # enable debuger +add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") # enable C++17 standard set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING") add_definitions(${CMAKE_CXX_FLAGS} "-O3 " ${ADDITIONAL_FLAGS}) # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") diff --git a/examples/crhmc_sampling/CMakeLists.txt b/examples/crhmc_sampling/CMakeLists.txt index 91fe46a7f..e3dddb905 100644 --- a/examples/crhmc_sampling/CMakeLists.txt +++ b/examples/crhmc_sampling/CMakeLists.txt @@ -113,7 +113,7 @@ else () endif () - add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard + add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") # enable C++17 standard set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING") add_definitions(${CMAKE_CXX_FLAGS} "-O3 -DTIME_KEEPING" ${ADDITIONAL_FLAGS}) # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") diff --git a/examples/ellipsoid-sampling/CMakeLists.txt b/examples/ellipsoid-sampling/CMakeLists.txt index b275ed9b9..d16704b17 100644 --- a/examples/ellipsoid-sampling/CMakeLists.txt +++ b/examples/ellipsoid-sampling/CMakeLists.txt @@ -99,7 +99,7 @@ else () endif () - add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard + add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") # enable C++11 standard add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") diff --git a/examples/hpolytope-volume/CMakeLists.txt b/examples/hpolytope-volume/CMakeLists.txt index eff88889a..ad9bb7e39 100644 --- a/examples/hpolytope-volume/CMakeLists.txt +++ b/examples/hpolytope-volume/CMakeLists.txt @@ -88,6 +88,7 @@ include_directories (BEFORE ../../include/nlp_oracles) include_directories (BEFORE ../../include/misc) include_directories (BEFORE ../../include/optimization) +add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-ldl") diff --git a/examples/logconcave/CMakeLists.txt b/examples/logconcave/CMakeLists.txt index 33725d0c9..cc62d286c 100644 --- a/examples/logconcave/CMakeLists.txt +++ b/examples/logconcave/CMakeLists.txt @@ -123,7 +123,7 @@ else () endif () - add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard + add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") # enable C++11 standard add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") diff --git a/examples/mmcs_method/CMakeLists.txt b/examples/mmcs_method/CMakeLists.txt index 01b8a040c..5ecd129f2 100644 --- a/examples/mmcs_method/CMakeLists.txt +++ b/examples/mmcs_method/CMakeLists.txt @@ -97,7 +97,7 @@ else () endif () - add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard + add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") # enable C++17 standard add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") diff --git a/examples/multithread_sampling/CMakeLists.txt b/examples/multithread_sampling/CMakeLists.txt index ded4b4dbb..81ad78285 100644 --- a/examples/multithread_sampling/CMakeLists.txt +++ b/examples/multithread_sampling/CMakeLists.txt @@ -97,7 +97,7 @@ else () endif () - #add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard + add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") # enable C++17 standard add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") diff --git a/examples/optimization_spectrahedra/CMakeLists.txt b/examples/optimization_spectrahedra/CMakeLists.txt index a35aee4c7..7c941640c 100644 --- a/examples/optimization_spectrahedra/CMakeLists.txt +++ b/examples/optimization_spectrahedra/CMakeLists.txt @@ -10,6 +10,7 @@ if(COMMAND cmake_policy) cmake_policy(SET CMP0003 NEW) endif(COMMAND cmake_policy) +add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") add_executable (read_write_sdpa_file read_write_sdpa_file.cpp) diff --git a/examples/order-polytope-basics/CMakeLists.txt b/examples/order-polytope-basics/CMakeLists.txt index b45ef374c..70026f61b 100644 --- a/examples/order-polytope-basics/CMakeLists.txt +++ b/examples/order-polytope-basics/CMakeLists.txt @@ -99,7 +99,7 @@ else () endif () - add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard + add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") # enable C++17 standard add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") diff --git a/examples/sampling-hpolytope-with-billiard-walks/CMakeLists.txt b/examples/sampling-hpolytope-with-billiard-walks/CMakeLists.txt index c10fb041c..5c47fe0ad 100644 --- a/examples/sampling-hpolytope-with-billiard-walks/CMakeLists.txt +++ b/examples/sampling-hpolytope-with-billiard-walks/CMakeLists.txt @@ -102,7 +102,7 @@ else () endif () - add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard + add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") # enable C++17 standard add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") diff --git a/examples/sampling-hpolytope-with-billiard-walks/sampler.cpp b/examples/sampling-hpolytope-with-billiard-walks/sampler.cpp index 22dd6d9e5..803d38dc3 100644 --- a/examples/sampling-hpolytope-with-billiard-walks/sampler.cpp +++ b/examples/sampling-hpolytope-with-billiard-walks/sampler.cpp @@ -73,17 +73,16 @@ void sample_using_gaussian_billiard_walk(HPOLYTOPE& HP, RNGType& rng, unsigned i unsigned int max_iter = 150; NT tol = std::pow(10, -6.0), reg = std::pow(10, -4.0); VT x0 = q.getCoefficients(); - std::pair, bool> inscribed_ellipsoid_res = max_inscribed_ellipsoid(HP.get_mat(), - HP.get_vec(), - x0, - max_iter, - tol, - reg); - if (!inscribed_ellipsoid_res.second) // not converged + MT E; + VT center; + bool converged; + std::tie(E, center, converged) = max_inscribed_ellipsoid(HP.get_mat(), + HP.get_vec(), x0, max_iter, tol, reg); + + if (!converged) // not converged throw std::runtime_error("max_inscribed_ellipsoid not converged"); - MT A_ell = inscribed_ellipsoid_res.first.first.inverse(); - EllipsoidType inscribed_ellipsoid(A_ell); + EllipsoidType inscribed_ellipsoid(E); // -------------------------------------------------------------------- Generator::apply(HP, q, inscribed_ellipsoid, num_points, walk_len, diff --git a/examples/volume_spectrahedron/CMakeLists.txt b/examples/volume_spectrahedron/CMakeLists.txt index 424c457b3..57ecd27e0 100644 --- a/examples/volume_spectrahedron/CMakeLists.txt +++ b/examples/volume_spectrahedron/CMakeLists.txt @@ -98,7 +98,7 @@ else () endif () - add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard + add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") # enable C++17 standard add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") diff --git a/examples/vpolytope-volume/CMakeLists.txt b/examples/vpolytope-volume/CMakeLists.txt index e1a0c2576..2f7020fb1 100644 --- a/examples/vpolytope-volume/CMakeLists.txt +++ b/examples/vpolytope-volume/CMakeLists.txt @@ -91,6 +91,7 @@ else () add_compile_definitions("EIGEN_NO_DEBUG") endif () +add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") diff --git a/include/convex_bodies/orderpolytope.h b/include/convex_bodies/orderpolytope.h index 2df5c6619..79f41e161 100644 --- a/include/convex_bodies/orderpolytope.h +++ b/include/convex_bodies/orderpolytope.h @@ -638,6 +638,12 @@ class OrderPolytope { return false; } + // Apply linear transformation, of square matrix T^{-1}, in H-polytope P:= Ax<=b + // This is most of the times for testing reasons because it might destroy the sparsity + void linear_transformIt(MT const& T) + { + _A = _A * T; + } // shift polytope by a point c void shift(VT const& c) diff --git a/include/generators/h_polytopes_generator.h b/include/generators/h_polytopes_generator.h index e9638b62f..0c9293baa 100644 --- a/include/generators/h_polytopes_generator.h +++ b/include/generators/h_polytopes_generator.h @@ -13,7 +13,7 @@ #include #include -#include "preprocess/max_inscribed_ellipsoid_rounding.hpp" +#include "preprocess/inscribed_ellipsoid_rounding.hpp" #ifndef isnan using std::isnan; @@ -121,7 +121,7 @@ Polytope skinny_random_hpoly(unsigned int dim, unsigned int m, const bool pre_ro if (pre_rounding) { Point x0(dim); // run only one iteration - max_inscribed_ellipsoid_rounding(P, x0, 1); + inscribed_ellipsoid_rounding(P, x0, 1); } MT cov = get_skinny_transformation(dim, eig_ratio, seed); diff --git a/include/preprocess/analytic_center_linear_ineq.h b/include/preprocess/analytic_center_linear_ineq.h index c7918635a..80d731f6c 100644 --- a/include/preprocess/analytic_center_linear_ineq.h +++ b/include/preprocess/analytic_center_linear_ineq.h @@ -12,7 +12,9 @@ #include -#include "max_inscribed_ball.hpp" +#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) @@ -33,22 +35,13 @@ 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 = 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 - H.noalias() = A_trans * s_sq.asDiagonal() * A; + update_Atrans_Diag_A(H, A_trans, A, s_sq.asDiagonal()); } /* @@ -66,39 +59,36 @@ void get_hessian_grad_logbarrier(MT const& A, MT const& A_trans, VT const& b, (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 + 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, - unsigned int const max_iters = 500, - NT const grad_err_tol = 1e-08, - NT const rel_pos_err_tol = 1e-12) +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) { - 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 + 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; - converged = false; + 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() = - H.llt().solve(grad); + 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)); @@ -128,7 +118,17 @@ std::tuple analytic_center_linear_ineq(MT const& A, VT const& b, } } while (true); - return std::make_tuple(x, converged); + 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/estimate_L_smooth_parameter.hpp b/include/preprocess/estimate_L_smooth_parameter.hpp index 8f5a1e89c..180f20947 100644 --- a/include/preprocess/estimate_L_smooth_parameter.hpp +++ b/include/preprocess/estimate_L_smooth_parameter.hpp @@ -1,7 +1,7 @@ // VolEsti (volume computation and sampling library) -// Copyright (c) 2012-2020 Vissarion Fisikopoulos -// Copyright (c) 2018-2020 Apostolos Chalkis +// 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. diff --git a/include/preprocess/feasible_point.hpp b/include/preprocess/feasible_point.hpp new file mode 100644 index 000000000..6e7f5db0a --- /dev/null +++ b/include/preprocess/feasible_point.hpp @@ -0,0 +1,34 @@ +// 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 FEASIBLE_POINT_HPP +#define FEASIBLE_POINT_HPP + +#include + +#include "preprocess/max_inscribed_ball.hpp" + +// Using MT as to deal with both dense and sparse matrices +template +VT compute_feasible_point(MT const& A, VT const& b) +{ + VT x; + bool feasibility_only = true, converged; + unsigned max_iters = 10000; + // Compute a feasible point + std::tie(x, std::ignore, converged) = max_inscribed_ball(A, b, max_iters, 1e-08, feasibility_only); + if (!converged || ((A * x).array() > b.array()).any()) + { + std::runtime_error("The computation of a feasible point failed."); + } + return x; +} + + +#endif diff --git a/include/preprocess/inscribed_ellipsoid_rounding.hpp b/include/preprocess/inscribed_ellipsoid_rounding.hpp new file mode 100644 index 000000000..dd20d533d --- /dev/null +++ b/include/preprocess/inscribed_ellipsoid_rounding.hpp @@ -0,0 +1,134 @@ +// VolEsti (volume computation and sampling library) + +// 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 INSCRIBED_ELLIPSOID_ROUNDING_HPP +#define INSCRIBED_ELLIPSOID_ROUNDING_HPP + +#include "preprocess/max_inscribed_ellipsoid.hpp" +#include "preprocess/analytic_center_linear_ineq.h" +#include "preprocess/feasible_point.hpp" + +enum EllipsoidType +{ + MAX_ELLIPSOID = 1, + LOG_BARRIER = 2 +}; + +template +inline static std::tuple +compute_inscribed_ellipsoid(Custom_MT A, VT b, VT const& x0, + unsigned int const& maxiter, + NT const& tol, NT const& reg) +{ + 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) + { + return analytic_center_linear_ineq(A, b, x0); + } else + { + std::runtime_error("Unknown rounding method."); + } + return {}; +} + +template +< + typename MT, + typename VT, + typename NT, + typename Polytope, + int ellipsoid_type = EllipsoidType::MAX_ELLIPSOID +> +std::tuple inscribed_ellipsoid_rounding(Polytope &P, + unsigned int const max_iterations = 5, + NT const max_eig_ratio = NT(6)) +{ + typedef typename Polytope::PointType Point; + VT x = compute_feasible_point(P.get_mat(), P.get_vec()); + return inscribed_ellipsoid_rounding(P, Point(x), max_iterations, max_eig_ratio); +} + +template +< + typename MT, + typename VT, + typename NT, + typename Polytope, + typename Point, + int ellipsoid_type = EllipsoidType::MAX_ELLIPSOID +> +std::tuple inscribed_ellipsoid_rounding(Polytope &P, + Point const& InnerPoint, + unsigned int const max_iterations = 5, + NT const max_eig_ratio = NT(6)) +{ + unsigned int maxiter = 500, iter = 1, d = P.dimension(); + VT x0 = InnerPoint.getCoefficients(), center, shift = VT::Zero(d); + MT E, L, T = MT::Identity(d, d); + bool converged; + NT R = 100.0, r = 1.0, tol = std::pow(10, -6.0), reg = std::pow(10, -4.0), round_val = 1.0; + + while (true) + { + // Compute the desired inscribed ellipsoid in P + std::tie(E, center, converged) = + compute_inscribed_ellipsoid(P.get_mat(), P.get_vec(), x0, maxiter, tol, reg); + + E = (E + E.transpose()) / 2.0; + E += MT::Identity(d, d)*std::pow(10, -8.0); //normalize E + + Eigen::LLT lltOfA(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of E^{-1} + L = lltOfA.matrixL(); + + // Computing eigenvalues of E + Spectra::DenseSymMatProd op(E); + // The value of ncv is chosen empirically + Spectra::SymEigsSolver> eigs(&op, 2, std::min(std::max(10, int(d)/5), int(d))); + eigs.init(); + int nconv = eigs.compute(); + if (eigs.info() == Spectra::COMPUTATION_INFO::SUCCESSFUL) { + R = 1.0 / eigs.eigenvalues().coeff(1); + r = 1.0 / eigs.eigenvalues().coeff(0); + } else { + Eigen::SelfAdjointEigenSolver eigensolver(E); + if (eigensolver.info() == Eigen::ComputationInfo::Success) { + R = 1.0 / eigensolver.eigenvalues().coeff(0); + r = 1.0 / eigensolver.eigenvalues().template tail<1>().value(); + } else { + std::runtime_error("Computations failed."); + } + } + // Shift polytope and apply the linear transformation on P + P.shift(center); + shift.noalias() += T * center; + T.applyOnTheRight(L); // T = T * L; + round_val *= L.transpose().determinant(); + P.linear_transformIt(L); + + reg = std::max(reg / 10.0, std::pow(10, -10.0)); + P.normalize(); + x0 = VT::Zero(d); + + // Check the roundness of the polytope + if(((std::abs(R / r) <= max_eig_ratio && converged) || iter >= max_iterations)) { + break; + } + + iter++; + } + + std::tuple result = std::make_tuple(T, shift, std::abs(round_val)); + return result; +} + +#endif diff --git a/include/preprocess/mat_computational_operators.h b/include/preprocess/mat_computational_operators.h new file mode 100644 index 000000000..af200d221 --- /dev/null +++ b/include/preprocess/mat_computational_operators.h @@ -0,0 +1,316 @@ +// 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 MAT_COMPUTATIONAL_OPERATORS_H +#define MAT_COMPUTATIONAL_OPERATORS_H + +#include + +#include "Spectra/include/Spectra/SymEigsSolver.h" +#include "Spectra/include/Spectra/MatOp/DenseSymMatProd.h" +#include "Spectra/include/Spectra/MatOp/SparseSymMatProd.h" + + +template +struct AssertFalseType : std::false_type {}; + +template +inline static auto +initialize_chol(MT const& mat) +{ + using DenseMT = Eigen::Matrix; + using SparseMT = Eigen::SparseMatrix; + if constexpr (std::is_same::value) + { + return std::make_unique>(); + } else if constexpr (std::is_same::value) + { + auto llt = std::make_unique>(); + llt->analyzePattern(mat); + return llt; + } else + { + static_assert(AssertFalseType::value, + "Matrix type is not supported."); + } +} + +template +inline static auto +initialize_chol(MT const& A_trans, MT const& A) +{ + using DenseMT = Eigen::Matrix; + using SparseMT = Eigen::SparseMatrix; + if constexpr (std::is_same::value) + { + return std::make_unique>(); + } else if constexpr (std::is_same::value) + { + MT mat = A_trans * A; + return initialize_chol(mat); + } else + { + static_assert(AssertFalseType::value, + "Matrix type is not supported."); + } +} + +template +inline static VT solve_vec(std::unique_ptr const& llt, + MT const& H, VT const& b) +{ + using DenseMT = Eigen::Matrix; + using SparseMT = Eigen::SparseMatrix; + if constexpr (std::is_same::value) + { + llt->compute(H); + return llt->solve(b); + } else if constexpr (std::is_same::value) + { + llt->factorize(H); + return llt->solve(b); + } else + { + static_assert(AssertFalseType::value, + "Matrix type is not supported."); + } +} + +template +inline static Eigen::Matrix +solve_mat(std::unique_ptr const& llt, + MT const& H, MT const& mat, NT &logdetE) +{ + using DenseMT = Eigen::Matrix; + using SparseMT = Eigen::SparseMatrix; + if constexpr (std::is_same::value) + { + llt->compute(H); + logdetE = llt->matrixL().toDenseMatrix().diagonal().array().log().sum(); + return llt->solve(mat); + } else if constexpr (std::is_same::value) + { + llt->factorize(H); + logdetE = llt->matrixL().nestedExpression().diagonal().array().log().sum(); + return llt->solve(mat); + } else + { + static_assert(AssertFalseType::value, + "Matrix type is not supported."); + } +} + +template +inline static void update_Atrans_Diag_A(MT &H, MT const& A_trans, + MT const& A, diag_MT const& D) +{ + using DenseMT = Eigen::Matrix; + using SparseMT = Eigen::SparseMatrix; + if constexpr (std::is_same::value) + { + H.noalias() = A_trans * D * A; + } else if constexpr (std::is_same::value) + { + H = A_trans * D * A; + } else + { + static_assert(AssertFalseType::value, + "Matrix type is not supported."); + } +} + +template +inline static void update_Diag_A(MT &H, diag_MT const& D, MT const& A) +{ + using DenseMT = Eigen::Matrix; + using SparseMT = Eigen::SparseMatrix; + if constexpr (std::is_same::value) + { + H.noalias() = D * A; + } else if constexpr (std::is_same::value) + { + H = D * A; + } else + { + static_assert(AssertFalseType::value, + "Matrix type is not supported."); + } +} + +template +inline static void update_A_Diag(MT &H, MT const& A, diag_MT const& D) +{ + using DenseMT = Eigen::Matrix; + using SparseMT = Eigen::SparseMatrix; + if constexpr (std::is_same::value) + { + H.noalias() = A * D; + } else if constexpr (std::is_same::value) + { + H = A * D; + } else + { + static_assert(AssertFalseType::value, + "Matrix type is not supported."); + } +} + +template +inline static auto +get_mat_prod_op(MT const& E) +{ + using DenseMT = Eigen::Matrix; + using SparseMT = Eigen::SparseMatrix; + if constexpr (std::is_same::value) + { + return std::make_unique>(E); + } else if constexpr (std::is_same::value) + { + return std::make_unique>(E); + } else + { + static_assert(AssertFalseType::value, + "Matrix type is not supported."); + } +} + +template +inline static auto get_eigs_solver(std::unique_ptr const& op, int const n) +{ + using DenseMatProd = Spectra::DenseSymMatProd; + using SparseMatProd = Spectra::SparseSymMatProd; + if constexpr (std::is_same::value) + { + using SymDenseEigsSolver = Spectra::SymEigsSolver + < + NT, + Spectra::SELECT_EIGENVALUE::BOTH_ENDS, + DenseMatProd + >; + // The value of ncv is chosen empirically + return std::make_unique(op.get(), 2, std::min(std::max(10, n/5), n)); + } else if constexpr (std::is_same::value) + { + using SymSparseEigsSolver = Spectra::SymEigsSolver + < + NT, + Spectra::SELECT_EIGENVALUE::BOTH_ENDS, + SparseMatProd + >; + // The value of ncv is chosen empirically + return std::make_unique(op.get(), 2, std::min(std::max(10, n/5), n)); + } else + { + static_assert(AssertFalseType::value, + "Matrix-vector multiplication multiplication is not supported."); + } +} + +template +inline static void +init_Bmat(MT &B, int const n, MT const& A_trans, MT const& A) +{ + using DenseMT = Eigen::Matrix; + using SparseMT = Eigen::SparseMatrix; + if constexpr (std::is_same::value) + { + B.resize(n+1, n+1); + } else if constexpr (std::is_same::value) + { + // Initialize the structure of matrix B + typedef Eigen::Triplet triplet; + std::vector trp; + for (int i = 0; i < n; i++) + { + trp.push_back(triplet(i, i, NT(1))); + trp.push_back(triplet(i, n, NT(1))); + trp.push_back(triplet(n, i, NT(1))); + } + trp.push_back(triplet(n, n, NT(1))); + + MT ATA = A_trans * A; + for (int k=0; k::value, + "Matrix type is not supported."); + } +} + +template +inline static void +update_Bmat(MT &B, VT const& AtDe, VT const& d, + MT const& AtD, MT const& A) +{ + using DenseMT = Eigen::Matrix; + using SparseMT = Eigen::SparseMatrix; + const int n = A.cols(); + /* + B is (n+1)x(n+1) and AtD_A is nxn. + We set B(1:n), 1:n) = AtD_A, B(n+1, :) = AtD^T, B(:, n+1) = AtD, B(n+1, n+1) = d.sum() + */ + if constexpr (std::is_same::value) + { + B.block(0, 0, n, n).noalias() = AtD * A; + B.block(0, n, n, 1).noalias() = AtDe; + B.block(n, 0, 1, n).noalias() = AtDe.transpose(); + B(n, n) = d.sum(); + B.noalias() += 1e-14 * MT::Identity(n + 1, n + 1); + } else if constexpr (std::is_same::value) + { + MT AtD_A = AtD * A; + int k = 0; + while(k < B.outerSize()) + { + typename MT::InnerIterator it2(AtD_A, k <= n-1 ? k : k-1); + for (typename MT::InnerIterator it1(B, k); it1; ++it1) + { + if (it1.row() <= n-1 && it1.col() <= n-1) + { + it1.valueRef() = it2.value(); + } + else if (it1.row() == n && it1.col() <= n-1) + { + it1.valueRef() = AtDe.coeff(it1.col()); + } + else if (it1.col() == n && it1.row() <= n-1) + { + it1.valueRef() = AtDe.coeff(it1.row()); + } + else // then, (it1.row() == n && it1.col() == n) + { + it1.valueRef() = d.sum(); + } + + if (it1.row() == it1.col()) + { + it1.valueRef() += 1e-14; + } + if (it1.row()::value, + "Matrix type is not supported."); + } +} + + +#endif // MAT_COMPUTATIONAL_OPERATORS_H diff --git a/include/preprocess/max_inscribed_ball.hpp b/include/preprocess/max_inscribed_ball.hpp index 672be24e4..573813dd7 100644 --- a/include/preprocess/max_inscribed_ball.hpp +++ b/include/preprocess/max_inscribed_ball.hpp @@ -1,15 +1,17 @@ // VolEsti (volume computation and sampling library) -// Copyright (c) 2012-2020 Vissarion Fisikopoulos -// Copyright (c) 2018-2020 Apostolos Chalkis +// 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 MAX_INNER_BALL -#define MAX_INNER_BALL +#ifndef MAX_INSCRIBED_BALL_HPP +#define MAX_INSCRIBED_BALL_HPP + +#include "preprocess/mat_computational_operators.h" /* This implmentation computes the largest inscribed ball in a given convex polytope P. @@ -26,10 +28,11 @@ radius r */ -template -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) +template +void calcstep(MT const& A, MT const& A_trans, MT const& B, + llt_type const& llt, 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) { int m = A.rows(), n = A.cols(); NT *vec_iter1 = tmp.data(), *vec_iter2 = y.data(), *vec_iter3 = s.data(), @@ -42,7 +45,7 @@ void calcstep(MT const& A, MT const& A_trans, Eigen::LLT const& lltOfB, VT & rhs.block(0,0,n,1).noalias() = r2 + A_trans * tmp; rhs(n) = r3 + tmp.sum(); - VT dxdt = lltOfB.solve(rhs); + VT dxdt = solve_vec(llt, B, rhs); dx = dxdt.block(0,0,n,1); dt = dxdt(n); @@ -56,12 +59,13 @@ void calcstep(MT const& A, MT const& A_trans, Eigen::LLT const& lltOfB, VT & } } - +// Using MT as to deal with both dense and sparse matrices template std::tuple max_inscribed_ball(MT const& A, VT const& b, unsigned int maxiter, NT tol, const bool feasibility_only = false) { + //typedef matrix_computational_operator mat_op; int m = A.rows(), n = A.cols(); bool converge = false; @@ -84,9 +88,10 @@ std::tuple max_inscribed_ball(MT const& A, VT const& b, 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), - eEye = std::pow(10.0, -14.0) * MT::Identity(n + 1, n + 1), - A_trans = A.transpose(); + MT B, AtD(n, m), A_trans = A.transpose(); + + init_Bmat(B, n, A_trans, A); + auto llt = initialize_chol(B); for (unsigned int i = 0; i < maxiter; ++i) { @@ -135,20 +140,13 @@ std::tuple max_inscribed_ball(MT const& A, VT const& b, vec_iter3++; vec_iter2++; } - AtD.noalias() = A_trans*d.asDiagonal(); + update_A_Diag(AtD, A_trans, d.asDiagonal()); // AtD = A_trans*d.asDiagonal() AtDe.noalias() = AtD * e_m; - B.block(0, 0, n, n).noalias() = AtD * A; - B.block(0, n, n, 1).noalias() = AtDe; - B.block(n, 0, 1, n).noalias() = AtDe.transpose(); - B(n, n) = d.sum(); - B.noalias() += eEye; - - // Cholesky decomposition - Eigen::LLT lltOfB(B); + update_Bmat(B, AtDe, d, AtD, A); // predictor step & length - calcstep(A, A_trans, lltOfB, s, y, r1, r2, r3, r4, dx, ds, dt, dy, tmp, rhs); + calcstep(A, A_trans, B, llt, s, y, r1, r2, r3, r4, dx, ds, dt, dy, tmp, rhs); alphap = -1.0; alphad = -1.0; @@ -175,7 +173,7 @@ std::tuple max_inscribed_ball(MT const& A, VT const& b, // corrector and combined step & length mu_ds_dy.noalias() = e_m * mu - ds.cwiseProduct(dy); - calcstep(A, A_trans, lltOfB, s, y, o_m, o_n, 0.0, mu_ds_dy, dxc, dsc, dtc, dyc, tmp, rhs); + calcstep(A, A_trans, B, llt, s, y, o_m, o_n, 0.0, mu_ds_dy, dxc, dsc, dtc, dyc, tmp, rhs); dx += dxc; ds += dsc; @@ -216,4 +214,4 @@ std::tuple max_inscribed_ball(MT const& A, VT const& b, return result; } -#endif +#endif // MAX_INSCRIBED_BALL_HPP diff --git a/include/preprocess/max_inscribed_ellipsoid.hpp b/include/preprocess/max_inscribed_ellipsoid.hpp index 7a0c19a6f..f44f60ff0 100644 --- a/include/preprocess/max_inscribed_ellipsoid.hpp +++ b/include/preprocess/max_inscribed_ellipsoid.hpp @@ -15,9 +15,7 @@ #include #include -#include "Spectra/include/Spectra/SymEigsSolver.h" -#include "Spectra/include/Spectra/Util/SelectionRule.h" -#include "Spectra/include/Spectra/MatOp/DenseSymMatProd.h" +#include "preprocess/mat_computational_operators.h" /* @@ -37,14 +35,14 @@ matrix E2^{-1} = E_transpose * E */ -// using Custom_MT as to deal with both dense and sparse matrices, MT will be the type of result matrix -// TODO: Change the return data structure to std::tuple -template -std::pair, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT const& x0, - unsigned int const& maxiter, - NT const& tol, NT const& reg) +// Using MT as to deal with both dense and sparse matrices, MT_dense will be the type of result matrix +template +std::tuple max_inscribed_ellipsoid(MT A, VT b, VT const& x0, + unsigned int const& maxiter, + NT const& tol, NT const& reg) { typedef Eigen::DiagonalMatrix Diagonal_MT; + //typedef matrix_computational_operator mat_op; int m = A.rows(), n = A.cols(); bool converged = false; @@ -66,22 +64,23 @@ std::pair, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT VT const bmAx0 = b - A * x0, ones_m = VT::Ones(m); - MT Q(m, m), E2(n, n), YQ(m,m), G(m,m), T(m,n), ATP(n,m), ATP_A(n,n); + MT_dense Q(m,m), 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); + MT YA(m, n); A = (ones_m.cwiseProduct(bmAx0.cwiseInverse())).asDiagonal() * A, b = ones_m; - Custom_MT A_trans = A.transpose(); + MT A_trans = A.transpose(), E2(n,n); + + auto llt = initialize_chol(A_trans, A); int i = 1; while (i <= maxiter) { Y = y.asDiagonal(); - E2.noalias() = MT(A_trans * Y * A); - Eigen::LLT llt(E2); - - Q.noalias() = A * llt.solve(A_trans); + update_Atrans_Diag_A(E2, A_trans, A, Y); + Q.noalias() = A * solve_mat(llt, E2, A_trans, logdetE2); + h = Q.diagonal(); h = h.cwiseSqrt(); @@ -120,7 +119,6 @@ std::pair, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT res = std::max(r1, r2); res = std::max(res, r3); - logdetE2 = llt.matrixL().toDenseMatrix().diagonal().array().log().sum(); objval = logdetE2; //logdet of E2 is already divided by 2 if (i % 10 == 0) { @@ -128,15 +126,13 @@ std::pair, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT NT rel, Rel; // computing eigenvalues of E2 - Spectra::DenseSymMatProd op(E2); - // The value of ncv is chosen empirically - Spectra::SymEigsSolver> eigs(&op, 2, std::min(std::max(10, n/5), n)); - eigs.init(); - int nconv = eigs.compute(); - if (eigs.info() == Spectra::COMPUTATION_INFO::SUCCESSFUL) { - Rel = 1.0 / eigs.eigenvalues().coeff(1); - rel = 1.0 / eigs.eigenvalues().coeff(0); + auto op = get_mat_prod_op(E2); + auto eigs = get_eigs_solver(op, n); + eigs->init(); + int nconv = eigs->compute(); + if (eigs->info() == Spectra::COMPUTATION_INFO::SUCCESSFUL) { + Rel = 1.0 / eigs->eigenvalues().coeff(1); + rel = 1.0 / eigs->eigenvalues().coeff(0); } else { Eigen::SelfAdjointEigenSolver eigensolver(E2); // E2 is positive definite matrix if (eigensolver.info() == Eigen::ComputationInfo::Success) { @@ -176,7 +172,7 @@ std::pair, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT YQ.noalias() = Y * Q; G = YQ.cwiseProduct(YQ.transpose()); y2h = 2.0 * yh; - YA.noalias() = Y * A; + update_Diag_A(YA, Y, A); // YA = Y * A; vec_iter1 = y2h.data(); vec_iter2 = z.data(); @@ -190,10 +186,10 @@ std::pair, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT G.diagonal() += y2h_z; h_z = h + z; - Eigen::PartialPivLU luG(G); - T.noalias() = luG.solve(h_z.asDiagonal()*YA); + Eigen::PartialPivLU luG(G); + T.noalias() = luG.solve(MT_dense(h_z.asDiagonal()*YA)); - ATP.noalias() = MT(y2h.asDiagonal()*T - YA).transpose(); + ATP.noalias() = MT_dense(y2h.asDiagonal()*T - YA).transpose(); vec_iter1 = R3.data(); vec_iter2 = y.data(); @@ -256,7 +252,7 @@ std::pair, bool> max_inscribed_ellipsoid(Custom_MT A, VT b, VT x += x0; } - return std::pair, bool>(std::pair(E2, x), converged); + return std::make_tuple(E2, x, converged); } diff --git a/include/preprocess/max_inscribed_ellipsoid_rounding.hpp b/include/preprocess/max_inscribed_ellipsoid_rounding.hpp deleted file mode 100644 index da00a8867..000000000 --- a/include/preprocess/max_inscribed_ellipsoid_rounding.hpp +++ /dev/null @@ -1,94 +0,0 @@ -// VolEsti (volume computation and sampling library) - -// Copyright (c) 2012-2020 Vissarion Fisikopoulos -// Copyright (c) 2018-2020 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 MAX_ELLIPSOID_ROUNDING_HPP -#define MAX_ELLIPSOID_ROUNDING_HPP - -#include "max_inscribed_ellipsoid.hpp" - -template -< - typename MT, - typename VT, - typename NT, - typename Polytope, - typename Point -> -std::tuple max_inscribed_ellipsoid_rounding(Polytope &P, - Point const& InnerPoint, - unsigned int const max_iterations = 5, - NT const max_eig_ratio = NT(6)) -{ - std::pair, bool> iter_res; - iter_res.second = false; - - VT x0 = InnerPoint.getCoefficients(); - MT E, L; - unsigned int maxiter = 500, iter = 1, d = P.dimension(); - - NT R = 100.0, r = 1.0, tol = std::pow(10, -6.0), reg = std::pow(10, -4.0), round_val = 1.0; - - MT T = MT::Identity(d, d); - VT shift = VT::Zero(d); - - 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); - E = iter_res.first.first; - E = (E + E.transpose()) / 2.0; - E = E + MT::Identity(d, d)*std::pow(10, -8.0); //normalize E - - Eigen::LLT lltOfA(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of E^{-1} - L = lltOfA.matrixL(); - - // computing eigenvalues of E - Spectra::DenseSymMatProd op(E); - // The value of ncv is chosen empirically - Spectra::SymEigsSolver> eigs(&op, 2, std::min(std::max(10, int(d)/5), int(d))); - eigs.init(); - int nconv = eigs.compute(); - if (eigs.info() == Spectra::COMPUTATION_INFO::SUCCESSFUL) { - R = 1.0 / eigs.eigenvalues().coeff(1); - r = 1.0 / eigs.eigenvalues().coeff(0); - } else { - Eigen::SelfAdjointEigenSolver eigensolver(E); - if (eigensolver.info() == Eigen::ComputationInfo::Success) { - R = 1.0 / eigensolver.eigenvalues().coeff(0); - r = 1.0 / eigensolver.eigenvalues().template tail<1>().value(); - } else { - std::runtime_error("Computations failed."); - } - } - // shift polytope and apply the linear transformation on P - P.shift(iter_res.first.second); - shift += T * iter_res.first.second; - T = T * L; - round_val *= L.transpose().determinant(); - P.linear_transformIt(L); - - reg = std::max(reg / 10.0, std::pow(10, -10.0)); - P.normalize(); - x0 = VT::Zero(d); - - // check the roundness of the polytope - if(((std::abs(R / r) <= max_eig_ratio && iter_res.second) || iter >= max_iterations)){ - break; - } - - iter++; - } - - std::tuple result = std::make_tuple(T, shift, std::abs(round_val)); - return result; -} - -#endif diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e17848fc6..cdda6087c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -289,13 +289,17 @@ add_executable (volume_cb_vpoly_intersection_vpoly volume_cb_vpoly_intersection_ add_test(NAME volume_cb_vpoly_intersection_vpoly_random_vpoly_sphere COMMAND volume_cb_vpoly_intersection_vpoly -tc=random_vpoly_sphere) -add_executable (new_rounding_test new_rounding_test.cpp $) +add_executable (rounding_test rounding_test.cpp $) add_test(NAME test_round_min_ellipsoid - COMMAND new_rounding_test -tc=round_min_ellipsoid) + COMMAND rounding_test -tc=round_min_ellipsoid) add_test(NAME test_round_max_ellipsoid - COMMAND new_rounding_test -tc=round_max_ellipsoid) + COMMAND rounding_test -tc=round_max_ellipsoid) add_test(NAME test_round_svd - COMMAND new_rounding_test -tc=round_svd) + COMMAND rounding_test -tc=round_svd) +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) @@ -357,6 +361,10 @@ 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) +add_test(NAME test_max_ball_sparse + COMMAND test_internal_points -tc=test_max_ball_sparse) + + set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING") @@ -380,7 +388,7 @@ TARGET_LINK_LIBRARIES(volume_cb_vpolytope lp_solve coverage_config) TARGET_LINK_LIBRARIES(volume_cb_zonotopes lp_solve coverage_config) TARGET_LINK_LIBRARIES(volume_cb_vpoly_intersection_vpoly lp_solve coverage_config) TARGET_LINK_LIBRARIES(volume_cb_vpoly_intersection_vpoly lp_solve ${MKL_LINK} coverage_config) -TARGET_LINK_LIBRARIES(new_rounding_test lp_solve ${MKL_LINK} coverage_config) +TARGET_LINK_LIBRARIES(rounding_test lp_solve ${MKL_LINK} coverage_config) TARGET_LINK_LIBRARIES(mcmc_diagnostics_test lp_solve ${MKL_LINK} coverage_config) TARGET_LINK_LIBRARIES(sampling_test lp_solve ${MKL_LINK} coverage_config) TARGET_LINK_LIBRARIES(mmcs_test lp_solve ${MKL_LINK} coverage_config) diff --git a/test/max_ellipsoid_rounding_test.cpp b/test/max_ellipsoid_rounding_test.cpp index a9c25b8c1..edce7a7d2 100644 --- a/test/max_ellipsoid_rounding_test.cpp +++ b/test/max_ellipsoid_rounding_test.cpp @@ -22,7 +22,7 @@ #include "volume/volume_cooling_gaussians.hpp" #include "volume/volume_cooling_balls.hpp" -#include "preprocess/max_inscribed_ellipsoid_rounding.hpp" +#include "preprocess/inscribed_ellipsoid_rounding.hpp" #include "generators/known_polytope_generators.h" @@ -64,7 +64,7 @@ void rounding_test(Polytope &HP, RNGType rng(d); std::pair InnerBall = HP.ComputeInnerBall(); - std::tuple res = max_inscribed_ellipsoid_rounding(HP, InnerBall.first); + std::tuple res = inscribed_ellipsoid_rounding(HP, InnerBall.first); // Setup the parameters int walk_len = 10; diff --git a/test/new_rounding_test.cpp b/test/rounding_test.cpp similarity index 61% rename from test/new_rounding_test.cpp rename to test/rounding_test.cpp index ab17ed939..7e15c99a5 100644 --- a/test/new_rounding_test.cpp +++ b/test/rounding_test.cpp @@ -22,7 +22,7 @@ #include "volume/volume_cooling_balls.hpp" #include "preprocess/min_sampling_covering_ellipsoid_rounding.hpp" -#include "preprocess/max_inscribed_ellipsoid_rounding.hpp" +#include "preprocess/inscribed_ellipsoid_rounding.hpp" #include "preprocess/svd_rounding.hpp" #include "generators/known_polytope_generators.h" @@ -110,8 +110,83 @@ void rounding_max_ellipsoid_test(Polytope &HP, typedef BoostRandomNumberGenerator RNGType; RNGType rng(d); std::pair InnerBall = HP.ComputeInnerBall(); - std::tuple res = max_inscribed_ellipsoid_rounding(HP, InnerBall.first); + 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_max_ellipsoid_sparse_test(double const& expectedBilliard, + double const& expected) +{ + typedef typename Polytope::PointType Point; + typedef typename Point::FT NT; + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Poset::RT RT; + typedef typename Poset::RV RV; + typedef Eigen::SparseMatrix SpMT; + + // Create Poset, 4 elements, a0 <= a1, a0 <= a2, a1 <= a3 + RV poset_data{{0, 1}, {0, 2}, {1, 3}}; + Poset poset(4, poset_data); + + // Initialize order polytope from the poset + OrderPolytope OP(poset); + OP.normalize(); + SpMT Asp = OP.get_mat(); + + + NT tol = 1e-08; + unsigned int maxiter = 500; + auto [center, radius, converged] = max_inscribed_ball(Asp, OP.get_vec(), maxiter, tol); + CHECK(OP.is_in(Point(center)) == -1); + auto [E, x0, round_val] = inscribed_ellipsoid_rounding(OP, Point(center)); + + MT A = MT(OP.get_mat()); + VT b = OP.get_vec(); + int d = OP.dimension(); + + Polytope HP(d, A, b); + + typedef BoostRandomNumberGenerator RNGType; + // 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 = round_val * volume_cooling_balls(HP, e, walk_len).second; + test_values(volume, expectedBilliard, expected); +} + +template +void rounding_log_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; @@ -164,11 +239,11 @@ void call_test_min_ellipsoid() { typedef HPolytope Hpolytope; Hpolytope P; - std::cout << "\n--- Testing rounding of H-skinny_cube5" << std::endl; + std::cout << "\n--- Testing min ellipsoid rounding of H-skinny_cube5" << std::endl; P = generate_skinny_cube(5); rounding_min_ellipsoid_test(P, 0, 3070.64, 3188.25, 3140.6, 3200.0); - std::cout << "\n--- Testing rounding of H-skinny_cube10" << std::endl; + std::cout << "\n--- Testing min ellipsoid rounding of H-skinny_cube10" << std::endl; P = generate_skinny_cube(10); rounding_min_ellipsoid_test(P, 0, 122550, 108426, 105003.0, 102400.0); @@ -182,11 +257,33 @@ void call_test_max_ellipsoid() { typedef HPolytope Hpolytope; Hpolytope P; - std::cout << "\n--- Testing rounding of H-skinny_cube5" << std::endl; + std::cout << "\n--- Testing max ellipsoid rounding of H-skinny_cube5" << std::endl; P = generate_skinny_cube(5); rounding_max_ellipsoid_test(P, 0, 3070.64, 3188.25, 3262.61, 3200.0); } +template +void call_test_max_ellipsoid_sparse() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef HPolytope Hpolytope; + + std::cout << "\n--- Testing max ellipsoid rounding of sparse Order Polytope" << std::endl; + rounding_max_ellipsoid_sparse_test(0.13979, 3070.64); +} + +template +void call_test_log_barrier() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "\n--- Testing log-barrier rounding of H-skinny_cube5" << std::endl; + P = generate_skinny_cube(5); + rounding_log_barrier_test(P, 0, 3070.64, 3188.25, 3262.77, 3200.0); +} + template void call_test_svd() { @@ -195,7 +292,7 @@ void call_test_svd() { typedef HPolytope Hpolytope; Hpolytope P; - std::cout << "\n--- Testing rounding of H-skinny_cube5" << std::endl; + std::cout << "\n--- Testing SVD rounding of H-skinny_cube5" << std::endl; P = generate_skinny_cube(5); rounding_svd_test(P, 0, 3070.64, 3188.25, 3140.6, 3200.0); } @@ -209,6 +306,14 @@ TEST_CASE("round_max_ellipsoid") { call_test_max_ellipsoid(); } +TEST_CASE("round_max_ellipsoid_sparse") { + call_test_max_ellipsoid_sparse(); +} + +TEST_CASE("round_log_barrier_test") { + call_test_log_barrier(); +} + TEST_CASE("round_svd") { call_test_svd(); } diff --git a/test/test_internal_points.cpp b/test/test_internal_points.cpp index fa44e9baf..8d4df0240 100644 --- a/test/test_internal_points.cpp +++ b/test/test_internal_points.cpp @@ -22,6 +22,9 @@ #include "generators/known_polytope_generators.h" #include "generators/h_polytopes_generator.h" +#include "convex_bodies/orderpolytope.h" +#include "misc/poset.h" + template void call_test_max_ball() { typedef Cartesian Kernel; @@ -40,12 +43,52 @@ void call_test_max_ball() { 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_sparse() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef HPolytope Hpolytope; + typedef boost::mt19937 PolyRNGType; + typedef typename OrderPolytope::VT VT; + typedef typename OrderPolytope::MT MT; + typedef typename Poset::RT RT; + typedef typename Poset::RV RV; + typedef Eigen::SparseMatrix SpMT; + + // Create Poset, 4 elements, a0 <= a1, a0 <= a2, a1 <= a3 + RV poset_data{{0, 1}, {0, 2}, {1, 3}}; + Poset poset(4, poset_data); + + // Initialize order polytope from the poset + OrderPolytope OP(poset); + OP.normalize(); + SpMT Asp = OP.get_mat(); + MT A = MT(OP.get_mat()); + VT b = OP.get_vec(); + + std::cout << "\n--- Testing Chebychev ball for sparse order Polytope" << std::endl; + + NT tol = 1e-08; + unsigned int maxiter = 500; + auto [center, radius, converged] = max_inscribed_ball(Asp, b, maxiter, tol); + auto [center2, radius2, converged2] = max_inscribed_ball(A, b, maxiter, tol); + + VT center_(4); + center_ << 0.207107, 0.5, 0.593398, 0.792893; + + CHECK(OP.is_in(Point(center)) == -1); + auto [E, x0, round_val] = inscribed_ellipsoid_rounding(OP, Point(center)); + + CHECK((center - center_).norm() <= 1e-06); + CHECK(std::abs(radius - 0.207107) <= 1e-06); + CHECK(converged); +} + template void call_test_max_ball_feasibility() { typedef Cartesian Kernel; @@ -77,6 +120,7 @@ void call_test_analytic_center() { typedef typename Hpolytope::MT MT; typedef typename Hpolytope::VT VT; typedef boost::mt19937 PolyRNGType; + typedef Eigen::SparseMatrix SpMT; Hpolytope P; std::cout << "\n--- Testing analytic center for skinny H-polytope" << std::endl; @@ -85,13 +129,24 @@ void call_test_analytic_center() { 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()); + auto [Hessian, analytic_center, converged] = analytic_center_linear_ineq(P.get_mat(), P.get_vec()); + SpMT Asp = P.get_mat().sparseView(); + auto [Hessian_sp, analytic_center2, converged2] = analytic_center_linear_ineq(Asp, 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); + + CHECK(P.is_in(Point(analytic_center2)) == -1); + CHECK(converged2); + CHECK(std::abs(analytic_center(0) - analytic_center2(0)) < 1e-12); + CHECK(std::abs(analytic_center(1) - analytic_center2(1)) < 1e-12); + CHECK(std::abs(analytic_center(2) - analytic_center2(2)) < 1e-12); + + CHECK((Hessian - Hessian_sp).norm() < 1e-12); } TEST_CASE("test_max_ball") { @@ -105,3 +160,7 @@ TEST_CASE("test_feasibility_point") { TEST_CASE("test_analytic_center") { call_test_analytic_center(); } + +TEST_CASE("test_max_ball_sparse") { + call_test_max_ball_sparse(); +}