Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Sparse optimizations for GaBW sampling #331

Open
wants to merge 3 commits into
base: develop
Choose a base branch
from
Open

Sparse optimizations for GaBW sampling #331

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ce1739a
simplify gabw sampling
lucaperju Sep 3, 2024
19565dc
optimize gabw for sparse polytopes
lucaperju Sep 6, 2024
4fc2b3f
minor changes
lucaperju Sep 14, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 33 additions & 8 deletions include/convex_bodies/hpolytope.h
View file
Open in desktop
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> {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is update_parameters here? I do not understand !std::is_same_v<update_parameters could you please explain?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

basically there's another compute_reflection function above which takes an integer (just the facet) as the 3rd argument, and for some reason, if I don't have that condition the compiler decides to call this function assuming that the typename of update_parameters is int. There might maybe be better ways of dealing with these issues, but I remember I tried to solve them for some time and this is the best I could do, I couldn't at all understand how the compiler chooses which function to use when there's multiple ones that match

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems that we are too generic here and that complicates the interface a lot! compute_reflection is called by several walks. Instead of creating all those complicated overloads why not simply create a new name. Especially if this is the case that this function is only called by a single walk (i.e. accelerated billiard walk). Moreover, update_parameters should not be a template but the struct defined in billiard walk, in all other cases this code will not compile since all other "update_parameters" does not have a moved_dist field.

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,39 @@ class HPolytope {
}
}

template <typename update_parameters>
NT compute_reflection(Point &v, const Point &, DenseMT const &AE, VT const &AEA, NT const &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 {
lucaperju marked this conversation as resolved.
Show resolved Hide resolved
lucaperju marked this conversation as resolved.
Show resolved Hide resolved

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
115 changes: 82 additions & 33 deletions include/random_walks/gaussian_accelerated_billiard_walk.hpp
View file
Open in desktop
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 >;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

similar comment to AE_type, is there a better naming?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hmm, I'll think of one, I'm not really sure what name I could give it but I'll see if I can come up with a better name.

// AA is sparse colMajor if MT is sparse rowMajor, and Dense otherwise
lucaperju marked this conversation as resolved.
Show resolved Hide resolved
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 >;
// ( ( ( )) ( ( ) ) ( ) )
lucaperju marked this conversation as resolved.
Show resolved Hide resolved
// 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());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should we explicitly cast it to DenseMT?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure what you mean, for the optimizations I need it to be in colmajor SparseMatrix format

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the question is: "Where is the problem with _AA.noalias() = 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) {
lucaperju marked this conversation as resolved.
Show resolved Hide resolved
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,33 +192,60 @@ 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);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not inserting the new values in the heap (in O(logn)) instead of rebuilding (in O(n))?

Copy link
Contributor Author

@lucaperju lucaperju Sep 13, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

here, this happens after we set a new direction, in which case now all the values have changed, so it's quicker to rebuild (O(n)) rather than insert each one (O(nlogn)). I'm not entirely sure if it makes a difference, but I think it does, since afterwards I never do O(nlogn) things, just O(non_zeroes * logn), so this O(nlogn) could be a bottleneck

} else {
_p += (_lambda_prev * _v);
}
T -= _lambda_prev;

T = T * coef;
coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
T = T / coef;

it++;

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;
coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
T = T / coef;

it++;
}
_p += _update_parameters.moved_dist * _v;
_update_parameters.moved_dist = 0.0;
if (it == _rho){
_p = p0;
was_reset = true;
}
}

Expand All @@ -219,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 @@ -251,15 +296,15 @@ struct GaussianAcceleratedBilliardWalk
_lambda_prev = dl * pbpair.first;
_p += (_lambda_prev * _v);
T -= _lambda_prev;

T = T * coef;
coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
T = T / coef;

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;

if (T <= pbpair.first) {
_p += (T * _v);
Expand All @@ -273,7 +318,11 @@ struct GaussianAcceleratedBilliardWalk
_lambda_prev = dl * pbpair.first;
_p += (_lambda_prev * _v);
T -= _lambda_prev;

T = T * coef;
coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
T = T / coef;

it++;
}
}
Expand All @@ -282,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;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is defined in uniform ABW, it should be better if it is defined in a separate file and both walks include it.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yeah, I was thinking about that too, any suggestions for that file name?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

boundary_oracle_heap?

};

};
Expand Down
2 changes: 1 addition & 1 deletion include/random_walks/uniform_accelerated_billiard_walk.hpp
View file
Open in desktop
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
Open in desktop
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
Loading