Gaussian Accelerated Billiard Walk #320
Conversation
lucaperju
commented
Jul 12, 2024
lucaperju
commented
- New Gaussian Accelerated Billiard Walk Sampling
- Fix bug causing Uniform Accelerated Billiard Walk to sometimes go outside of the polytope
- Remove unused variable
* 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 <apostolos.chalkis@quantagonia.com>
* 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 <papachristoumarios@gmail.com> Co-authored-by: Apostolos Chalkis <apostolos.chalkis@quantagonia.com>
…ce parameter, add define lp_solve guards to orderpolytope
Remove volesti4dingo branch
…_source Build only with lp-solve source code, no need for liblpsolve55.so files
* Enable c++17 support in tests and fix saxpy calls * Fix structure, code style and tests in autodiff
* correcting the includes to point the boost library directory * modifs to include only the include/ and external/ directories * modifs * Update CMakeLists.txt --------- Co-authored-by: Vissarion Fisikopoulos <fisikop@gmail.com> Co-authored-by: vfisikop <160006479+vfisikop@users.noreply.github.com>
Remove R interface
Refactor trigonometric_positive_intersect function for hpolytopes
update_position complexity improvement
* improve the implementation of maximum ellipsoid computation * minor fix in rounding unit-test * fix file copyrights * apply termination criterion after transformation in max ellipsoid rounding * resolve PR comments --------- Co-authored-by: Apostolos Chalkis <apostolos.chalkis@quantagonia.com>
* generalize rounding loop * support sparse cholesky operator * complete sparse support in max_inscribed_ball * complete sparse support in preprocesing * add sparse tests * change main rounding function name * improve explaining comments * resolve PR comments * changing the dates in copyrights * use if constexpr instead of SNIFAE * update the examples to cpp17 * update to cpp17 order polytope example * fix templating in mat_computational_operators * fix templating errors and change header file to mat_computational_operators * first implementation of the volumetric barrier ellipsoid * add criterion for step_iter * restructure code that computes barriers' centers * remove unused comments * resolve PR comments * remove NT typename from max_step() --------- Co-authored-by: Apostolos Chalkis <apostolos.chalkis@quantagonia.com>
include/convex_bodies/hpolytope.h
Outdated
|
|
||
| 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)) |
There was a problem hiding this comment.
4.0 * params.inner_vi_ak can be computed only once
|
|
||
|
|
||
| template | ||
| < |
There was a problem hiding this comment.
could please fix the indentation here?
| initialize(P, p, E, rng); | ||
| } | ||
|
|
||
| template <typename GenericPolytope> |
There was a problem hiding this comment.
Could you please comment on the parts of the code that are different from uniform billiard walk and which ones are identical? Does it makes sense to reuse parts of the code instead of repeating ?
There was a problem hiding this comment.
Mainly, the only different parts are the way the vector is chosen each time, and the reflection function, there are also a few added lines in the initialization for computing the AE and AEA matrices, and a few lines after each line_positive_intersect which uses _AA to adjust for the scaling coefficient. I don't really know how we could reuse parts of the code from uniform billiard walk, other than if we were to replace it. Since, in theory, this should be the same as uniform abw if the given covariance matrix is the identity, but that would change how the walk is called, and might introduce other complications
| it = 0; | ||
| std::pair<NT, int> pbpair; | ||
| if(!was_reset) | ||
| pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, |
| RandomNumberGenerator &rng) | ||
| { | ||
| _update_parameters = update_parameters(); | ||
| _Len = compute_diameter<GenericPolytope> |
| @@ -104,8 +104,16 @@ struct AcceleratedBilliardWalk | |||
| Point p0 = _p; | |||
|
|
|||
| it = 0; | |||
| std::pair<NT, int> pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, | |||
| std::pair<NT, int> pbpair; | |||
| if(!was_reset) | |||
There was a problem hiding this comment.
Good point, thanks! Could you please write it following the style of the rest of the file? i.e.
if (!was_reset) {
...
}
| T = -std::log(rng.sample_urdist()) * _L; | ||
|
|
||
|
|
||
| Eigen::LLT<MT> lltOfE(E.inverse()); // compute the Cholesky decomposition of inv(E) |
There was a problem hiding this comment.
Why do we need to compute the cholesky decomposition of E in every iteration since it should be considered fixed?
We could compute _L_cov in the constructor and initialize a member variable _E. Then, E shouldn't be an input in apply() neither in initialize().
| T = -std::log(rng.sample_urdist()) * _L; | ||
|
|
||
|
|
||
| Eigen::LLT<MT> lltOfE(E.inverse()); // compute the Cholesky decomposition of inv(E) |
There was a problem hiding this comment.
plz replace Eigen::LLT<MT> lltOfE(E.inverse()); with
Eigen::LLT<MT> lltOfE(E.llt().solve(MT::Identity(E.cols(), E.cols())));
| for(int i = 0; i < P.num_of_hyperplanes(); ++i) | ||
| { | ||
| _AEA(i) = _AE.row(i).dot(P.get_mat().row(i)); | ||
| } |
There was a problem hiding this comment.
This probably is equivalent to
_AEA = _AE.cwiseProduct(P.get_mat()).rowwise().sum();
Plz check.
|
|
||
|
|
|
|
||
|
|
* Order Polytopes generation * minor changes for PR * remove space and comment * remove bug * Unit test for Random Order Polytope, and minor changes * remove comment * Order Polytopes generation * minor changes for PR * remove space and comment * remove bug * Unit test for Random Order Polytope, and minor changes * remove comment * fix rebase bugs * remove accidental line
