From 19565dc1eacb2dd43b2f69f042675099b5549710 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Fri, 6 Sep 2024 14:05:57 +0200 Subject: [PATCH] optimize gabw for sparse polytopes --- include/convex_bodies/hpolytope.h | 40 +++++-- .../gaussian_accelerated_billiard_walk.hpp | 102 ++++++++++++------ .../uniform_accelerated_billiard_walk.hpp | 2 +- test/sampling_test.cpp | 22 +++- 4 files changed, 119 insertions(+), 47 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index f2648aceb..03be9bbbc 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -1010,9 +1010,10 @@ class HPolytope { // updates the velocity vector v and the position vector p after a reflection // the real value of p is given by p + moved_dist * v + // MT must be sparse, in RowMajor format template auto compute_reflection(Point &v, Point &p, update_parameters const& params) const - -> std::enable_if_t> && !std::is_same_v, void> { // MT must be in RowMajor format + -> std::enable_if_t> && !std::is_same_v, void> { NT* v_data = v.pointerToData(); NT* p_data = p.pointerToData(); for(Eigen::SparseMatrix::InnerIterator it(A, params.facet_prev); it; ++it) { @@ -1021,15 +1022,38 @@ class HPolytope { } } - template - NT compute_reflection(Point &v, const Point &, DenseMT const &AE, VT const &AEA, NT &vEv, update_parameters const ¶ms) const { + // function to compute reflection for GaBW random walk + // compatible when the polytope is both dense or sparse + template + NT compute_reflection(Point &v, Point &p, AE_type const &AE, VT const &AEA, NT &vEv, update_parameters const ¶ms) const { + + NT new_vEv; + if constexpr (!std::is_same_v>) { + Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev)); + VT x = v.getCoefficients(); + new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev)); + v += a; + } else { + + if constexpr(!std::is_same_v>) { + VT x = v.getCoefficients(); + new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev)); + } else { + new_vEv = vEv + 4.0 * params.inner_vi_ak * params.inner_vi_ak * AEA(params.facet_prev); + NT* v_data = v.pointerToData(); + for(Eigen::SparseMatrix::InnerIterator it(AE, params.facet_prev); it; ++it) { + new_vEv -= 4.0 * params.inner_vi_ak * it.value() * *(v_data + it.col()); + } + } - Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev)); - VT x = v.getCoefficients(); - NT new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev)); - v += a; + NT* v_data = v.pointerToData(); + NT* p_data = p.pointerToData(); + for(Eigen::SparseMatrix::InnerIterator it(A, params.facet_prev); it; ++it) { + *(v_data + it.col()) += (-2.0 * params.inner_vi_ak) * it.value(); + *(p_data + it.col()) -= (-2.0 * params.inner_vi_ak * params.moved_dist) * it.value(); + } + } NT coeff = std::sqrt(vEv / new_vEv); - // v = v * coeff; vEv = new_vEv; return coeff; } diff --git a/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp index 5176e3a78..b2d1b9644 100644 --- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp +++ b/include/random_walks/gaussian_accelerated_billiard_walk.hpp @@ -47,12 +47,13 @@ struct GaussianAcceleratedBilliardWalk struct update_parameters { update_parameters() - : facet_prev(0), hit_ball(false), inner_vi_ak(0.0), ball_inner_norm(0.0) + : facet_prev(0), hit_ball(false), inner_vi_ak(0.0), ball_inner_norm(0.0), moved_dist(0.0) {} int facet_prev; bool hit_ball; double inner_vi_ak; double ball_inner_norm; + double moved_dist; }; parameters param; @@ -67,9 +68,15 @@ struct GaussianAcceleratedBilliardWalk struct Walk { typedef typename Polytope::PointType Point; - typedef typename Polytope::DenseMT DenseMT; + typedef typename Polytope::MT MT; + typedef typename Eigen::Matrix DenseMT; typedef typename Polytope::VT VT; typedef typename Point::FT NT; + using AA_type = std::conditional_t< std::is_same_v>, typename Eigen::SparseMatrix, DenseMT >; + // AA is sparse colMajor if MT is sparse rowMajor, and Dense otherwise + using AE_type = std::conditional_t< std::is_same_v> && std::is_base_of, E_type>::value, typename Eigen::SparseMatrix, DenseMT >; + // ( ( ( )) ( ( ) ) ( ) ) + // AE is sparse rowMajor if (MT is sparse rowMajor and E is sparse), and Dense otherwise void computeLcov(const E_type E) { @@ -110,7 +117,11 @@ struct GaussianAcceleratedBilliardWalk _L = compute_diameter::template compute(P); computeLcov(E); _E = E; - _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose()); + if constexpr (std::is_same>::value) { + _AA = (P.get_mat() * P.get_mat().transpose()); + } else { + _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose()); + } _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) initialize(P, p, rng); } @@ -131,7 +142,11 @@ struct GaussianAcceleratedBilliardWalk ::template compute(P); computeLcov(E); _E = E; - _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose()); + if constexpr (std::is_same>::value) { + _AA = (P.get_mat() * P.get_mat().transpose()); + } else { + _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose()); + } _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) initialize(P, p, rng); } @@ -148,6 +163,12 @@ struct GaussianAcceleratedBilliardWalk int it; NT coef = 1.0; NT vEv; + typename Point::Coeff b; + NT* b_data; + if constexpr (std::is_same>::value) { + b = P.get_vec(); + b_data = b.data(); + } for (auto j=0u; j pbpair; - if(!was_reset) { - pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _update_parameters); - } else { - pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters); - was_reset = false; - } + std::pair pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters); if (T <= pbpair.first) { _p += (T * _v); @@ -177,7 +192,18 @@ struct GaussianAcceleratedBilliardWalk } _lambda_prev = dl * pbpair.first; - _p += (_lambda_prev * _v); + if constexpr (std::is_same>::value) { + _update_parameters.moved_dist = _lambda_prev; + NT* Ar_data = _lambdas.data(); + NT* Av_data = _Av.data(); + for(int i = 0; i < P.num_of_hyperplanes(); ++i) { + _distances_set.vec[i].first = ( *(b_data + i) - (*(Ar_data + i)) ) / (*(Av_data + i)); + } + // rebuild the heap with the new values of (b - Ar) / Av + _distances_set.rebuild(_update_parameters.moved_dist); + } else { + _p += (_lambda_prev * _v); + } T -= _lambda_prev; T = T * coef; @@ -188,11 +214,14 @@ struct GaussianAcceleratedBilliardWalk while (it < _rho) { - std::pair pbpair - = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters); - //_Av *= coef; - //_update_parameters.inner_vi_ak *= coef; - //pbpair.first /= coef; + std::pair pbpair; + if constexpr (std::is_same>::value) { + pbpair = P.line_positive_intersect(_p, _lambdas, _Av, _lambda_prev, + _distances_set, _AA, _update_parameters); + } else { + pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, + _AA, _update_parameters); + } if (T <= pbpair.first) { _p += (T * _v); @@ -200,7 +229,11 @@ struct GaussianAcceleratedBilliardWalk break; } _lambda_prev = dl * pbpair.first; - _p += (_lambda_prev * _v); + if constexpr (std::is_same>::value) { + _update_parameters.moved_dist += _lambda_prev; + } else { + _p += (_lambda_prev * _v); + } T -= _lambda_prev; T = T * coef; @@ -209,9 +242,10 @@ struct GaussianAcceleratedBilliardWalk it++; } + _p += _update_parameters.moved_dist * _v; + _update_parameters.moved_dist = 0.0; if (it == _rho){ _p = p0; - was_reset = true; } } @@ -227,21 +261,24 @@ struct GaussianAcceleratedBilliardWalk { unsigned int n = P.dimension(); const NT dl = 0.995; - was_reset = false; _lambdas.setZero(P.num_of_hyperplanes()); _Av.setZero(P.num_of_hyperplanes()); _p = p; - _AE.noalias() = (DenseMT)(P.get_mat() * _E); - _AEA = _AE.cwiseProduct((DenseMT)P.get_mat()).rowwise().sum(); - /* - _AEA.resize(P.num_of_hyperplanes()); - for(int i = 0; i < P.num_of_hyperplanes(); ++i) { - _AEA(i) = _AE.row(i).dot(P.get_mat().row(i)); - }*/ + DenseMT temp_AE; + temp_AE.noalias() = (DenseMT)(P.get_mat() * _E); + _AEA = temp_AE.cwiseProduct((DenseMT)P.get_mat()).rowwise().sum(); + if constexpr (std::is_same_v) + _AE = temp_AE; + else + _AE = temp_AE.sparseView(); + } + _distances_set = BoundaryOracleHeap(P.num_of_hyperplanes()); + _v = GetDirection::apply(n, rng, false); _v = Point(_L_cov.template triangularView() * _v.getCoefficients()); + NT T = -std::log(rng.sample_urdist()) * _L; Point p0 = _p; @@ -268,9 +305,6 @@ struct GaussianAcceleratedBilliardWalk { std::pair pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters); - //_Av *= coef; - //_update_parameters.inner_vi_ak *= coef; - //pbpair.first /= coef; if (T <= pbpair.first) { _p += (T * _v); @@ -288,7 +322,7 @@ struct GaussianAcceleratedBilliardWalk T = T * coef; coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); T = T / coef; - + it++; } } @@ -297,16 +331,16 @@ struct GaussianAcceleratedBilliardWalk Point _p; Point _v; NT _lambda_prev; - DenseMT _AA; + AA_type _AA; E_type _L_cov; // LL' = inv(E) - DenseMT _AE; + AE_type _AE; E_type _E; VT _AEA; unsigned int _rho; update_parameters _update_parameters; typename Point::Coeff _lambdas; typename Point::Coeff _Av; - bool was_reset; + BoundaryOracleHeap _distances_set; }; }; diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index f91abee4b..77660c745 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -275,7 +275,7 @@ struct AcceleratedBilliardWalk it++; while (it < _rho) - { + { std::pair pbpair; if constexpr (std::is_same>::value) { pbpair = P.line_positive_intersect(_p, _lambdas, _Av, _lambda_prev, diff --git a/test/sampling_test.cpp b/test/sampling_test.cpp index a15f40102..ba71d1398 100644 --- a/test/sampling_test.cpp +++ b/test/sampling_test.cpp @@ -320,8 +320,11 @@ void call_test_gabw(){ Point StartingPoint(d); std::list randPoints; + std::chrono::time_point start, stop; + start = std::chrono::high_resolution_clock::now(); + std::cout << "--- Testing Gaussian Accelerated Billiard Walk for Skinny-H-cube10" << std::endl; - P = generate_skinny_cube(10); + P = generate_skinny_cube(d); Point p = P.ComputeInnerBall().first; @@ -340,6 +343,11 @@ void call_test_gabw(){ RandomPointGenerator::apply(P, p, E, numpoints, 1, randPoints, push_back_policy, rng); + stop = std::chrono::high_resolution_clock::now(); + + std::chrono::duration total_time = stop - start; + std::cout << "Done in " << total_time.count() << '\n'; + MT samples(d, numpoints); unsigned int jj = 0; for (typename std::list::iterator rpit = randPoints.begin(); rpit != randPoints.end(); rpit++, jj++) @@ -352,11 +360,11 @@ void call_test_gabw(){ RNGType Srng(d); typedef Eigen::SparseMatrix SparseMT; - typedef HPolytope SparseHpoly; + typedef HPolytope> SparseHpoly; std::list Points; SparseHpoly SP; - SP = generate_skinny_cube(10); + SP = generate_skinny_cube(d); std::cout << "--- Testing Gaussian Accelerated Billiard Walk for Sparse Skinny-H-cube10" << std::endl; @@ -371,15 +379,21 @@ void call_test_gabw(){ > sparsewalk; typedef MultivariateGaussianRandomPointGenerator SparseRandomPointGenerator; + start = std::chrono::high_resolution_clock::now(); ellipsoid = compute_inscribed_ellipsoid - (SP.get_mat(), SP.get_vec(), p.getCoefficients(), 500, std::pow(10, -6.0), std::pow(10, -4.0)); + ((SparseMT)SP.get_mat(), SP.get_vec(), p.getCoefficients(), 500, std::pow(10, -6.0), std::pow(10, -4.0)); const SparseMT SE = get<0>(ellipsoid).sparseView(); SparseRandomPointGenerator::apply(SP, p, SE, numpoints, 1, Points, push_back_policy, Srng); + stop = std::chrono::high_resolution_clock::now(); + + total_time = stop - start; + std::cout << "Done in " << total_time.count() << '\n'; + jj = 0; MT Ssamples(d, numpoints); for (typename std::list::iterator rpit = Points.begin(); rpit != Points.end(); rpit++, jj++)