Skip to content

Commit

Permalink
optimize gabw for sparse polytopes
Browse files Browse the repository at this point in the history
  • Loading branch information
@lucaperju
lucaperju committed Sep 6, 2024
1 parent ce1739a commit 19565dc
Show file tree
Hide file tree
Showing 4 changed files with 119 additions and 47 deletions.
40 changes: 32 additions & 8 deletions include/convex_bodies/hpolytope.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <typename update_parameters>
auto compute_reflection(Point &v, Point &p, update_parameters const& params) const
-> std::enable_if_t<std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>> && !std::is_same_v<update_parameters, int>, void> { // MT must be in RowMajor format
-> std::enable_if_t<std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>> && !std::is_same_v<update_parameters, int>, void> {
NT* v_data = v.pointerToData();
NT* p_data = p.pointerToData();
for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, params.facet_prev); it; ++it) {
Expand All @@ -1021,15 +1022,38 @@ class HPolytope {
}
}

template <typename update_parameters>
NT compute_reflection(Point &v, const Point &, DenseMT const &AE, VT const &AEA, NT &vEv, update_parameters const &params) const {
// function to compute reflection for GaBW random walk
// compatible when the polytope is both dense or sparse
template <typename AE_type, typename update_parameters>
NT compute_reflection(Point &v, Point &p, AE_type const &AE, VT const &AEA, NT &vEv, update_parameters const &params) const {

NT new_vEv;
if constexpr (!std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>) {
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<AE_type, Eigen::SparseMatrix<NT, Eigen::RowMajor>>) {
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<double, Eigen::RowMajor>::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<double, Eigen::RowMajor>::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;
}
Expand Down
102 changes: 68 additions & 34 deletions include/random_walks/gaussian_accelerated_billiard_walk.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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<double, Eigen::Dynamic, Eigen::Dynamic> DenseMT;
typedef typename Polytope::VT VT;
typedef typename Point::FT NT;
using AA_type = std::conditional_t< std::is_same_v<MT, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>>, typename Eigen::SparseMatrix<NT>, DenseMT >;
// AA is sparse colMajor if MT is sparse rowMajor, and Dense otherwise
using AE_type = std::conditional_t< std::is_same_v<MT, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>> && std::is_base_of<Eigen::SparseMatrixBase<E_type>, E_type>::value, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>, DenseMT >;
// ( ( ( )) ( ( ) ) ( ) )
// AE is sparse rowMajor if (MT is sparse rowMajor and E is sparse), and Dense otherwise

void computeLcov(const E_type E)
{
Expand Down Expand Up @@ -110,7 +117,11 @@ struct GaussianAcceleratedBilliardWalk
_L = compute_diameter<GenericPolytope>::template compute<NT>(P);
computeLcov(E);
_E = E;
_AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose());
if constexpr (std::is_same<AA_type, Eigen::SparseMatrix<NT>>::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);
}
Expand All @@ -131,7 +142,11 @@ struct GaussianAcceleratedBilliardWalk
::template compute<NT>(P);
computeLcov(E);
_E = E;
_AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose());
if constexpr (std::is_same<AA_type, Eigen::SparseMatrix<NT>>::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);
}
Expand All @@ -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<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
b = P.get_vec();
b_data = b.data();
}

for (auto j=0u; j<walk_length; ++j)
{
Expand All @@ -162,13 +183,7 @@ struct GaussianAcceleratedBilliardWalk
Point p0 = _p;

it = 0;
std::pair<NT, int> 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<NT, int> pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters);

if (T <= pbpair.first) {
_p += (T * _v);
Expand All @@ -177,7 +192,18 @@ struct GaussianAcceleratedBilliardWalk
}

_lambda_prev = dl * pbpair.first;
_p += (_lambda_prev * _v);
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::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;
Expand All @@ -188,19 +214,26 @@ struct GaussianAcceleratedBilliardWalk

while (it < _rho)
{
std::pair<NT, int> 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<NT, int> pbpair;
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::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);
_lambda_prev = T;
break;
}
_lambda_prev = dl * pbpair.first;
_p += (_lambda_prev * _v);
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
_update_parameters.moved_dist += _lambda_prev;
} else {
_p += (_lambda_prev * _v);
}
T -= _lambda_prev;

T = T * coef;
Expand All @@ -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;
}
}

Expand All @@ -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_type, DenseMT>)
_AE = temp_AE;
else
_AE = temp_AE.sparseView();
}
_distances_set = BoundaryOracleHeap<NT>(P.num_of_hyperplanes());


_v = GetDirection<Point>::apply(n, rng, false);
_v = Point(_L_cov.template triangularView<Eigen::Lower>() * _v.getCoefficients());


NT T = -std::log(rng.sample_urdist()) * _L;
Point p0 = _p;
Expand All @@ -268,9 +305,6 @@ struct GaussianAcceleratedBilliardWalk
{
std::pair<NT, int> 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);
Expand All @@ -288,7 +322,7 @@ struct GaussianAcceleratedBilliardWalk
T = T * coef;
coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
T = T / coef;

it++;
}
}
Expand All @@ -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<NT> _distances_set;
};

};
Expand Down
2 changes: 1 addition & 1 deletion include/random_walks/uniform_accelerated_billiard_walk.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -275,7 +275,7 @@ struct AcceleratedBilliardWalk
it++;

while (it < _rho)
{
{
std::pair<NT, int> pbpair;
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
pbpair = P.line_positive_intersect(_p, _lambdas, _Av, _lambda_prev,
Expand Down
22 changes: 18 additions & 4 deletions test/sampling_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -320,8 +320,11 @@ void call_test_gabw(){
Point StartingPoint(d);
std::list<Point> randPoints;

std::chrono::time_point<std::chrono::high_resolution_clock> 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<Hpolytope>(10);
P = generate_skinny_cube<Hpolytope>(d);


Point p = P.ComputeInnerBall().first;
Expand All @@ -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<double> total_time = stop - start;
std::cout << "Done in " << total_time.count() << '\n';

MT samples(d, numpoints);
unsigned int jj = 0;
for (typename std::list<Point>::iterator rpit = randPoints.begin(); rpit != randPoints.end(); rpit++, jj++)
Expand All @@ -352,11 +360,11 @@ void call_test_gabw(){
RNGType Srng(d);

typedef Eigen::SparseMatrix<NT> SparseMT;
typedef HPolytope<Point, SparseMT> SparseHpoly;
typedef HPolytope<Point, Eigen::SparseMatrix<NT, Eigen::RowMajor>> SparseHpoly;
std::list<Point> Points;

SparseHpoly SP;
SP = generate_skinny_cube<SparseHpoly>(10);
SP = generate_skinny_cube<SparseHpoly>(d);


std::cout << "--- Testing Gaussian Accelerated Billiard Walk for Sparse Skinny-H-cube10" << std::endl;
Expand All @@ -371,15 +379,21 @@ void call_test_gabw(){
> sparsewalk;
typedef MultivariateGaussianRandomPointGenerator <sparsewalk> SparseRandomPointGenerator;

start = std::chrono::high_resolution_clock::now();

ellipsoid = compute_inscribed_ellipsoid<MT, EllipsoidType::MAX_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<Point>::iterator rpit = Points.begin(); rpit != Points.end(); rpit++, jj++)
Expand Down

0 comments on commit 19565dc

Please sign in to comment.