No-U-Turn Sampler in HMC #216

TolisChal · 2022-03-10T13:11:11Z

This PR implements the No-U-Turn Sampler (NUTS) in HMC:

-- Implements a new structure for the random walk.
-- Optimizes the number of gradient computations in both HMC and NUTS; it reduces this number to half per step.
-- Implements a new C++ test for NUTS in ./test/logconcave_sampling_test.cpp, namely benchmark_nuts_hmc (for both truncated and non-truncated cases).
-- Extends the R interface to expose the NUTS sampler in R.
-- Implements a new R example in ./R-proj/examples/logconcave/nuts_rand_poly.R.

This PR resolves the issue #123

papachristoumarios · 2022-03-10T19:59:54Z

R-proj/examples/logconcave/simple_hmc_rand_poly.R

@@ -29,10 +29,10 @@ dimension <- 50
 facets <- 200

 # Create domain of truncation
-H <- gen_rand_hpoly(dimension, facets)
+H <- gen_rand_hpoly(dimension, facets, seed = 15)


Why a specific seed here?

For some seeds, the polytope is not bounded and the script fails.

So, I fixed the seed to test a specific instance and check if something changes/fails
after a commit.

papachristoumarios · 2022-03-10T20:00:46Z

R-proj/src/sample_points.cpp

@@ -242,7 +263,7 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo
 //' @param n The number of points that the function is going to sample from the convex polytope.
 //' @param random_walk Optional. A list that declares the random walk and some related parameters as follows:
 //' \itemize{
-//' \item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, v) \code{'dikin'} for dikin walk, vi) \code{'vaidya'} for vaidya walk, vii) \code{'john'} for john walk, viii) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or ix) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR x) \code{'HMC'} for Hamiltonian Monte Carlo (logconcave densities) xi) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method xii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (spherical Gaussian or exponential distribution). The default walk is \code{'aBiW'} for the uniform distribution or \code{'CDHR'} for the Gaussian distribution and H-polytopes and \code{'BiW'} or \code{'RDHR'} for the same distributions and V-polytopes and zonotopes.}
+//' \item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, v) \code{'dikin'} for dikin walk, vi) \code{'vaidya'} for vaidya walk, vii) \code{'john'} for john walk, viii) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or ix) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR x) \code{'NUTS'} for NUTS Hamiltonian Monte Carlo sampler (logconcave densities) xi) \code{'HMC'} for Hamiltonian Monte Carlo (logconcave densities) xii) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method (logconcave densities) xiii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (spherical Gaussian or exponential distribution). The default walk is \code{'aBiW'} for the uniform distribution, \code{'CDHR'} for the Gaussian distribution and H-polytopes and \code{'BiW'} or \code{'RDHR'} for the same distributions and V-polytopes and zonotopes. \code{'NUTS'} is the default sampler for logconcave densities.}


This line is too long. Can we split it?

Yes, I will do so. Thank you.

papachristoumarios · 2022-03-10T20:07:17Z

include/ode_solvers/leapfrog.hpp

@@ -56,6 +56,8 @@ struct LeapfrogODESolver {
  pts xs;
  pts xs_prev;

+  Point grad_x;


There is probably something going on with the solver. I run one of the examples at examples/logconcave/simple_hmc.cpp (I also added NUTS example which I am going to push in a PR).

While NUTS seemed to return the correct marginals (in 2D), HMC returned the following wrong marginal.
The density is π(x) \propto exp(-f(x)) with f(x) = 2 x^T x + x.sum() for x belonging to a 2D-cube.

Also, have your changes affected how eta is set on vanilla hmc?

Thank you for this comment.
I think the current commit works properly considering HMC.
I implemented the same optimization as in NUTS for the number of evaluations of the gradient.
In particular, the second evaluation in the leapfrog's step is used for the first momenta update in the next leapfrog step.

I did not change anything in HMC parameterization.

papachristoumarios · 2022-03-10T20:09:00Z

include/preprocess/estimate_L_smooth_parameter.hpp

+        randPoints[0] = p;
+
+        listOfPoints.push_back(randPoints);
+        //std::cout<<(listOfPoints[i])[0].getCoefficients().transpose()<<std::endl;


Please remove unused lines.

Thank you. I did so.

papachristoumarios · 2022-03-10T20:13:37Z

include/preprocess/estimate_L_smooth_parameter.hpp

+    //std::cout<<"length = "<<listOfPoints.size()<<std::endl;
+    NT L = std::numeric_limits<NT>::lowest(), Ltemp;
+
+    for (int i=0; i<rnum-1; i++)


Can't you just try subsequent points, ignoring the times the sampler stays at a place. This way you would also not need to store the points.

Yes, but what about not subsequent points?

papachristoumarios · 2022-03-10T20:15:14Z

include/random_walks/hamiltonian_monte_carlo_walk.hpp

@@ -125,8 +125,10 @@ struct HamiltonianMonteCarloWalk {
      // Pick a random velocity
      v = GetDirection<Point>::apply(dim, rng, false);

-      solver->set_state(0, x);


See my comment above. HMC returns wrong marginals.

I think the current commit fixes this issue.

papachristoumarios · 2022-03-10T20:16:33Z

include/random_walks/nuts_hmc_walk.hpp

@@ -0,0 +1,384 @@
+// VolEsti (volume computation and sampling library)


Have you also tried to involve the average number of reflections on the burn-in method?

No, I did not. I added a TODO comment to generalize Nesterov's algorithm in the truncated setting.

papachristoumarios

@TolisChal Thank you for the PR! NUTS seems functional from the simple examples I have tried it on (we should also try more examples)!

However, the HMC walk seems broken and returns the wrong marginals (see more on the corresponding comment), which I think is due to the changes in the leapfrog integrator.

Thank you!

papachristoumarios

Hi all,

The functionality of vanilla HMC seems to be restored.

Thank you!

vissarion

Thanks for this PR. That's a really cool feature! I have a few comments on the code.

vissarion · 2022-04-04T14:27:11Z

R-proj/examples/logconcave/nuts_rand_poly.R

+# Sample points
+n_samples <- 20000
+
+samples <- sample_points(P, n = n_samples, random_walk = list("walk" = "NUTS", "solver" = "leapfrog", "starting_point" = warm_start[,1]), distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f))


Could you reduce the length here?

done. Thanks

vissarion · 2022-04-04T14:36:42Z

include/cartesian_geom/point.h

@@ -79,6 +79,7 @@ class point
    void set_dimension(const unsigned int dim)
    {
        d = dim;
+        coeffs.setZero(d);


now this is not just setting a new dimension but also reset the whole vector/point. Should this be renamed to resize_point or something?

vissarion · 2022-04-04T14:37:20Z

include/ode_solvers/leapfrog.hpp

@@ -101,24 +103,25 @@ struct LeapfrogODESolver {

  void step(int k, bool accepted) {
    num_steps++;
-
+    


please remove trailing spaces

include/ode_solvers/leapfrog.hpp

vissarion · 2022-04-05T13:50:02Z

include/random_walks/nuts_hmc_walk.hpp

+      NT epsilon_=2)
+    {
+      epsilon = epsilon_;
+      if (F.params.L > 0){


or simpler eta = F.params.L > 0 ? 10.0 / (dim * sqrt(F.params.L)) : 0.005;

vissarion · 2022-04-05T13:53:12Z

include/random_walks/nuts_hmc_walk.hpp

+      accepted = false;
+
+      // Initialize solver
+      solver = new Solver(0, params.eta, pts{x, x}, F, bounds{P, NULL});


please use a smart pointer here to avoid memory leaks

vissarion · 2022-04-05T14:18:54Z

include/random_walks/nuts_hmc_walk.hpp

+
+      NT uu = std::log(rng.sample_urdist()) - h1;
+      int j = -1;
+      bool s = true;


what is s, does it make sense to give a more descriptive name ?

this comes from the paper. In general, I use the variable names from the nuts paper

vissarion · 2022-04-05T14:23:25Z

include/sampling/random_point_generators.hpp

    >
-    static void apply(Polytope &P,


I guess this is an fix, unrelated to NUTS.

vissarion · 2022-04-05T14:24:54Z

test/CMakeLists.txt

@@ -298,6 +298,11 @@ else ()
          COMMAND logconcave_sampling_test -tc=uld)
  add_test(NAME logconcave_sampling_test_exponential_biomass_sampling
          COMMAND logconcave_sampling_test -tc=exponential_biomass_sampling)
+  add_test(NAME logconcave_sampling_test_nuts_hmc_truncated


I think it would be helpful to also add a C++ example of NUTS.

@vissarion @TolisChal I have already opened this PR to Tolis' fork with a NUTS example: TolisChal#29

Great! It seems that you've resolved the mkl linking issue in examples ;-)
e.g. as described here #214

vissarion · 2022-11-02T10:37:24Z

Hi @TolisChal what is the status of this PR? Should you merge it after fixing the conflicts?

* Enable github actions to build examples. Avoid passing a polytope as a const reference. * Fix ambiguous call to fix function by renaming volesti's diagnostic function. (GeomScale#263) * Updating documentation (GeomScale#261) Adding WSL and MKL build instructions. * disable an R sampling test for windows * Fix the warning message in R Mac's cran test (GeomScale#285) * copy and replace lp_rlp.h * remove re-initialization of eta * update ubuntu version from 18 to 20 in R cran tests * minor improvements (explanatory comments) --------- Co-authored-by: Apostolos Chalkis <[email protected]> * delete commented out code * No-U-Turn Sampler in HMC (GeomScale#216) * initialize nuts sampler class * implement the burnin of nuts sampler * add tests and resolve bugs * implement e0 estimation in burnin of nuts * optimize leapfrog * integrate nuts into the R interface * document random walk in sample_points.cpp in R interface * fix burnin for the non-truncated case * resolve bugs in hmc and nuts pipelines * improve the preprocess in burin step of nuts * split large lines in sample_points.cpp * Add NUTS C++ example and update CMake (GeomScale#29) * resolve PR comments * fix minor bug * fix compiler bug * fix error in building the C++ examples * resolve warnings in sample_points * fix lpsolve cran warning * fix cran warning on mac * improve lpsolve cmake for cran check * fix R warning in mac test * remove lpsolve header * resolve PR comments --------- Co-authored-by: Marios Papachristou <[email protected]> Co-authored-by: Apostolos Chalkis <[email protected]> --------- Co-authored-by: Vissarion Fisikopoulos <[email protected]> Co-authored-by: Soumya Tarafder <[email protected]> Co-authored-by: Apostolos Chalkis <[email protected]> Co-authored-by: Marios Papachristou <[email protected]>

* initialize nuts sampler class * implement the burnin of nuts sampler * add tests and resolve bugs * implement e0 estimation in burnin of nuts * optimize leapfrog * integrate nuts into the R interface * document random walk in sample_points.cpp in R interface * fix burnin for the non-truncated case * resolve bugs in hmc and nuts pipelines * improve the preprocess in burin step of nuts * split large lines in sample_points.cpp * Add NUTS C++ example and update CMake (#29) * resolve PR comments * fix minor bug * fix compiler bug * fix error in building the C++ examples * resolve warnings in sample_points * fix lpsolve cran warning * fix cran warning on mac * improve lpsolve cmake for cran check * fix R warning in mac test * remove lpsolve header * resolve PR comments --------- Co-authored-by: Marios Papachristou <[email protected]> Co-authored-by: Apostolos Chalkis (TolisChal) <[email protected]>

* initialize nuts sampler class * implement the burnin of nuts sampler * add tests and resolve bugs * implement e0 estimation in burnin of nuts * optimize leapfrog * integrate nuts into the R interface * document random walk in sample_points.cpp in R interface * fix burnin for the non-truncated case * resolve bugs in hmc and nuts pipelines * improve the preprocess in burin step of nuts * split large lines in sample_points.cpp * Add NUTS C++ example and update CMake (#29) * resolve PR comments * fix minor bug * fix compiler bug * fix error in building the C++ examples * resolve warnings in sample_points * fix lpsolve cran warning * fix cran warning on mac * improve lpsolve cmake for cran check * fix R warning in mac test * remove lpsolve header * resolve PR comments --------- Co-authored-by: Marios Papachristou <[email protected]> Co-authored-by: Apostolos Chalkis <[email protected]>

* initialize nuts sampler class * implement the burnin of nuts sampler * add tests and resolve bugs * implement e0 estimation in burnin of nuts * optimize leapfrog * integrate nuts into the R interface * document random walk in sample_points.cpp in R interface * fix burnin for the non-truncated case * resolve bugs in hmc and nuts pipelines * improve the preprocess in burin step of nuts * split large lines in sample_points.cpp * Add NUTS C++ example and update CMake (#29) * resolve PR comments * fix minor bug * fix compiler bug * fix error in building the C++ examples * resolve warnings in sample_points * fix lpsolve cran warning * fix cran warning on mac * improve lpsolve cmake for cran check * fix R warning in mac test * remove lpsolve header * resolve PR comments --------- Co-authored-by: Marios Papachristou <[email protected]> Co-authored-by: Apostolos Chalkis (TolisChal) <[email protected]>

TolisChal added 7 commits March 9, 2022 00:24

initialize nuts sampler class

92fe314

implement the burnin of nuts sampler

9fc6567

add tests and resolve bugs

e7f4675

implement e0 estimation in burnin of nuts

158dc39

optimize leapfrog

351936d

integrate nuts into the R interface

dfdf2e3

document random walk in sample_points.cpp in R interface

4076760

vissarion requested review from vissarion and papachristoumarios March 10, 2022 13:51

vissarion added the enhancement label Mar 10, 2022

vissarion linked an issue Mar 10, 2022 that may be closed by this pull request

Implement NUTS #123

Closed

fix burnin for the non-truncated case

d05af8f

papachristoumarios reviewed Mar 10, 2022

View reviewed changes

papachristoumarios suggested changes Mar 10, 2022

View reviewed changes

TolisChal added 2 commits March 11, 2022 16:08

resolve bugs in hmc and nuts pipelines

97bca4a

improve the preprocess in burin step of nuts

ae80164

papachristoumarios approved these changes Mar 13, 2022

View reviewed changes

split large lines in sample_points.cpp

578b0b5

papachristoumarios approved these changes Mar 14, 2022

View reviewed changes

vissarion reviewed Apr 5, 2022

View reviewed changes

papachristoumarios and others added 4 commits June 4, 2022 12:15

Add NUTS C++ example and update CMake (#29)

0f1bb2f

Merge branch 'develop' into nuts_hmc

adcf382

resolve PR comments

7bdf929

fix minor bug

51ae29a

fix compiler bug

44d7d17

Apostolos Chalkis added 10 commits October 13, 2023 23:12

resolve conflicts

bd0fbb5

fix error in building the C++ examples

5752b36

resolve warnings in sample_points

a70f978

fix lpsolve cran warning

39d72ad

fix cran warning on mac

40cf11e

improve lpsolve cmake for cran check

5e9fad6

fix R warning in mac test

274c02f

remove lpsolve header

30db634

resolve conflicts

c248779

resolve PR comments

12061f2

TolisChal merged commit 074a562 into GeomScale:develop Oct 17, 2023
27 checks passed

TolisChal mentioned this pull request Nov 1, 2023

Reflective Hamiltonian Monte Carlo with Leafrog steps for exponential sampling #143

Closed

TolisChal deleted the nuts_hmc branch June 21, 2024 19:14

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

No-U-Turn Sampler in HMC #216

No-U-Turn Sampler in HMC #216

TolisChal commented Mar 10, 2022 •

edited

Loading

papachristoumarios Mar 10, 2022

TolisChal Mar 14, 2022

papachristoumarios Mar 10, 2022

TolisChal Mar 14, 2022

papachristoumarios Mar 10, 2022 •

edited

Loading

papachristoumarios Mar 10, 2022

TolisChal Mar 14, 2022

papachristoumarios Mar 10, 2022

TolisChal Mar 14, 2022

papachristoumarios Mar 10, 2022

TolisChal Mar 14, 2022

papachristoumarios Mar 10, 2022

TolisChal Mar 14, 2022

papachristoumarios Mar 10, 2022

TolisChal Mar 14, 2022

papachristoumarios left a comment

papachristoumarios left a comment

vissarion left a comment

vissarion Apr 4, 2022

TolisChal Jul 20, 2022

vissarion Apr 4, 2022

vissarion Apr 4, 2022

vissarion Apr 5, 2022

TolisChal Oct 17, 2023

vissarion Apr 5, 2022

vissarion Apr 5, 2022

TolisChal Oct 17, 2023 •

edited

Loading

vissarion Apr 5, 2022

vissarion Apr 5, 2022

papachristoumarios Apr 5, 2022

vissarion Apr 8, 2022 •

edited

Loading

vissarion commented Nov 2, 2022

		@@ -0,0 +1,384 @@
		// VolEsti (volume computation and sampling library)

		@@ -101,24 +103,25 @@ struct LeapfrogODESolver {

		void step(int k, bool accepted) {
		num_steps++;

No-U-Turn Sampler in HMC #216

No-U-Turn Sampler in HMC #216

Conversation

TolisChal commented Mar 10, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

papachristoumarios Mar 10, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

papachristoumarios left a comment

Choose a reason for hiding this comment

papachristoumarios left a comment

Choose a reason for hiding this comment

vissarion left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

TolisChal Oct 17, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

vissarion Apr 8, 2022 • edited Loading

Choose a reason for hiding this comment

vissarion commented Nov 2, 2022

TolisChal commented Mar 10, 2022 •

edited

Loading

papachristoumarios Mar 10, 2022 •

edited

Loading

TolisChal Oct 17, 2023 •

edited

Loading

vissarion Apr 8, 2022 •

edited

Loading