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.

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

Exact Hamiltonian Monte Carlo for spherical Gaussian sampling #144

Merged
merged 57 commits into from
Oct 3, 2021

Conversation

TolisChal
Copy link
Member

@TolisChal TolisChal commented Mar 8, 2021

This PR implements exact Hamiltonian Monte Carlo for spherical Gaussian sampling.

Algorithmic details:
-- The exact HMC for spherical Gaussian sampling exploits trigonometric trajectories in the interior of the polytope.
-- When the trajectory hits the boundary it is reflected.

Implementation details:
-- The sampler is implemented by a new structure.
-- The PR exploits the C++ interface for Gaussian sampling.
-- The sampler is integrated into the R interface of volesti.

This PR depends on #142

Copy link
Collaborator

@papachristoumarios papachristoumarios left a comment

Choose a reason for hiding this comment

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

My review is done for this PR. All in all my comments are almost identical to #142 .
Finally, it seems that there's some code in this PR that should belong to #143 .

Thank you!

VT& Av,
int& facet_prev) const
{
return std::make_pair(0, 0);
Copy link
Collaborator

Choose a reason for hiding this comment

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

Exceptions should be raised as in #142. The same applies for all classes where std::make_pair(0, 0) is called by actually no output is computed.

Copy link
Member Author

Choose a reason for hiding this comment

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

@papachristoumarios Thank you for the comments!
This one changed according to #142


// Hamiltonian Monte Carlo with leapfrog for sampling from the exponential distribution

struct HMCLeafrogExponential
Copy link
Collaborator

Choose a reason for hiding this comment

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

Should this code belong to #143 ?

Copy link
Member Author

Choose a reason for hiding this comment

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

This structure removed. The #143 is going to be updated.


// Exact HMC for sampling from the spherical Gaussian distribution

struct ExactHMCGaussianWalk
Copy link
Collaborator

Choose a reason for hiding this comment

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

Naming consistency comments as in #142 .

Copy link
Member Author

Choose a reason for hiding this comment

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

Thanks. Done.

};



Copy link
Collaborator

Choose a reason for hiding this comment

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

Too much whitespace here :)

Copy link
Member Author

Choose a reason for hiding this comment

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

Thanks. I erased it.

@@ -20,6 +21,9 @@
#include "random_walks/uniform_john_walk.hpp"
#include "random_walks/uniform_vaidya_walk.hpp"
#include "random_walks/uniform_accelerated_billiard_walk.hpp"
#include "random_walks/exponential_hamiltonian_monte_carlo_walk.hpp"
Copy link
Collaborator

Choose a reason for hiding this comment

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

Header naming consistency. 

Copy link
Member Author

Choose a reason for hiding this comment

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

Thanks. Done as mentioned in #142

@vissarion vissarion mentioned this pull request Sep 21, 2021
Copy link
Member

@vissarion vissarion left a comment

Choose a reason for hiding this comment

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

Thanks!

int m = num_of_hyperplanes();
int facet = -1;

Ar.noalias() = A * r.getCoefficients();
Copy link
Member

Choose a reason for hiding this comment

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

the two functions quadratic_positive_intersect differ only on how Ar is computed. Should you combine them?

Copy link
Member Author

Choose a reason for hiding this comment

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

I implemented a new function get_positive_quadratic_root. Both oracles are using this function now.

const double tol = 1e-10;
if (lamda1 * lamda2 < NT(0))
{
if (previous_facet == current_facet)
Copy link
Member

Choose a reason for hiding this comment

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

this can be simplified

lamda = previous_facet == current_facet 
    ? std::max(lamda1, lamda2) < NT(tol) ? std::min(lamda1, lamda2) : std::max(lamda1, lamda2)
    : std::max(lamda1, lamda2);

Copy link
Member

Choose a reason for hiding this comment

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

Also min max values can be precomputed (e.g. using std::minmax)

Copy link
Member Author

Choose a reason for hiding this comment

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

Done, thanks for the code.

sum_nom = Ar - b;
Av.noalias() = A * v.getCoefficients();

NT* sum_nom_data = sum_nom.data();
Copy link
Member

Choose a reason for hiding this comment

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

why are you using pointers (returned by data()) here?

Copy link
Member Author

Choose a reason for hiding this comment

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

This is the standard approach in all boundary oracle - function in hpolytope.h.

t1 = (2.0 * M_PI) / omega;
}

t2 = (-std::acos((*b_data) / C) - Phi) / omega;
Copy link
Member

Choose a reason for hiding this comment

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

std::acos((*b_data) / C) can be precomputed

Copy link
Member Author

Choose a reason for hiding this comment

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

done

#ifndef RANDOM_WALKS_EXPONENTIAL_EXACT_HMC_WALK_HPP
#define RANDOM_WALKS_EXPONENTIAL_EXACT_HMC_WALK_HPP

#define IN_INSIDE_BODY_TOLLERANCE 1e-10
Copy link
Member

Choose a reason for hiding this comment

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

shouldn't be named INSIDE_BODY_TOLLERANCE?

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes, I agree. Done

for (unsigned int i=0; i<rnum; ++i)
{
success = walk.template apply(P, p, walk_length, rng);
if (!success) {
Copy link
Member

Choose a reason for hiding this comment

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

should an exception be thrown here instead?

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes, done.

for (unsigned int i=0; i<rnum; ++i)
{
success = walk.template apply(P, p, walk_length, rng);
if (!success) {
Copy link
Member

Choose a reason for hiding this comment

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

similar here

Copy link
Member Author

Choose a reason for hiding this comment

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

done

if (nburns > 0) {
RandomPointGenerator::apply(P, p, nburns, walk_len, randPoints,
push_back_policy, rng);
randPoints.clear();
Copy link
Member

Choose a reason for hiding this comment

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

do you need to store the random points here? It looks like a waste of space (and time).

Copy link
Member Author

Choose a reason for hiding this comment

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

Please see #180

if (nburns > 0) {
RandomPointGenerator::apply(P, p, c, a, nburns, walk_len, randPoints,
push_back_policy, rng);
randPoints.clear();
Copy link
Member

Choose a reason for hiding this comment

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

similar to above

Copy link
Member Author

Choose a reason for hiding this comment

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

Please see #180

typename NT,
typename Point
>
void exponential_sampling(PointList &randPoints,
Copy link
Member

Choose a reason for hiding this comment

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

please avoid using the same code in all 4 exponential_sampling specializations

Copy link
Member Author

Choose a reason for hiding this comment

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

Please see #180.
We might include it there too.

Copy link
Member Author

Choose a reason for hiding this comment

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

Please see #180.
We might include it there too.

@vissarion vissarion added the dingo label Oct 1, 2021
@vissarion vissarion merged commit d74cdcb into GeomScale:develop Oct 3, 2021
@TolisChal TolisChal mentioned this pull request Oct 4, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants