|
| 1 | +#' Sample points from a convex polytope |
| 2 | +#' |
| 3 | +#' Sample perfect uniformly distributed points from well known convex bodies |
| 4 | +#' or run geometric random walks to approximate arbitrary distributions. |
| 5 | +#' |
| 6 | +#' The \eqn{d}-dimensional unit simplex is the set of points \eqn{\vec{x} \in \R^d}, |
| 7 | +#' s.t.: \eqn{\sum_i x_i \leq 1}, \eqn{x_i \geq 0}. The \eqn{d}-dimensional canonical |
| 8 | +#' simplex is the set of points \eqn{\vec{x} \in \R^d}, s.t.: \eqn{\sum_i x_i = 1}, |
| 9 | +#' \eqn{x_i \geq 0}. |
| 10 | +#' |
| 11 | +#' @param P A convex polytope object. It is an object from class (a) Hpolytope or |
| 12 | +#' (b) Vpolytope or (c) Zonotope or (d) VpolytopeIntersection. |
| 13 | +#' @param n The number of points that the function is going to sample. |
| 14 | +#' @param random_walk Optional. A list that declares the random walk method and |
| 15 | +#' related parameters as follows: |
| 16 | +#' \describe{ |
| 17 | +#' \item{\code{walk}}{A string: (a) \code{'CDHR'} for coordinate direction hit-and-run, |
| 18 | +#' b) \code{'RDHR'} for random direction hit-and-run, c) \code{'BaW'} for the ball walk, |
| 19 | +#' d) \code{'BiW'} for the billiard walk, e) \code{'BiRCDHR'} for the bicriteria random walk, |
| 20 | +#' f) \code{'Dikin'} for dikin walk, g) \code{'Vaidya'} for vaidya walk, |
| 21 | +#' h) \code{'John'} for john walk, i) \code{'BRDHR'} for boundary hit-and-run, |
| 22 | +#' j) \code{'Hamiltonian'} for Hamiltonian walk, k) \code{'accelerated_billiard'} for |
| 23 | +#' accelerated billiard walk, l) \code{'BilliardRef'} for billiard walk with reflections, |
| 24 | +#' m) \code{'HRlu'} for logconcave settings with H-polytope, |
| 25 | +#' n) \code{'HRnr'} for logconcave densities with V-polytope, o) \code{'HESSIANRWR'} |
| 26 | +#' for logconcave densities with H-polytope and sparse constrained problems.} |
| 27 | +#' \item{\code{walk_length}}{The number of the steps per generated point for the random walk. |
| 28 | +#' The default value is \eqn{1}.} |
| 29 | +#' \item{\code{nburns}}{The number of points to burn before start sampling. The default |
| 30 | +#' value is \eqn{1}.} |
| 31 | +#' \item{\code{starting_point}}{A \eqn{d}-dimensional numerical vector that declares a |
| 32 | +#' starting point in the interior of the polytope for the random walk. The default choice |
| 33 | +#' is the center of the ball as computed by \code{inner_ball()}.} |
| 34 | +#' \item{\code{BaW_rad}}{The radius for the ball walk.} |
| 35 | +#' \item{\code{L}}{The maximum length of the billiard trajectory or the radius for the step |
| 36 | +#' of Dikin, Vaidya or John walk.} |
| 37 | +#' \item{\code{solver}}{Specify ODE solver for logconcave sampling. Options are i) leapfrog, |
| 38 | +#' ii) euler iii) runge-kutta iv) richardson} |
| 39 | +#' \item{\code{step_size}}{Optionally chosen step size for logconcave sampling. Defaults to |
| 40 | +#' a theoretical value if not provided.} |
| 41 | +#' } |
| 42 | +#' @param distribution Optional. A list that declares the target density and related |
| 43 | +#' parameters as follows: |
| 44 | +#' \describe{ |
| 45 | +#' \item{\code{density}}{A string: (a) \code{'uniform'} for the uniform distribution, |
| 46 | +#' (b) \code{'gaussian'} for the multidimensional spherical distribution, |
| 47 | +#' (c) \code{'logconcave'} with form proportional to exp(-f(x)) where f(x) is L-smooth |
| 48 | +#' and m-strongly-convex, (d) \code{'exponential'} for the exponential distribution. |
| 49 | +#' The default target distribution is the uniform distribution.} |
| 50 | +#' \item{\code{variance}}{The variance of the multidimensional spherical Gaussian or the |
| 51 | +#' exponential distribution. The default value is 1.} |
| 52 | +#' \item{\code{mode}}{A \eqn{d}-dimensional numerical vector that declares the mode of the |
| 53 | +#' Gaussian distribution. The default choice is the center computed by \code{inner_ball()}.} |
| 54 | +#' \item{\code{bias}}{The bias vector for the exponential distribution. The default vector |
| 55 | +#' is \eqn{c_1 = 1} and \eqn{c_i = 0} for \eqn{i \neq 1}.} |
| 56 | +#' \item{\code{L_}}{Smoothness constant (for logconcave).} |
| 57 | +#' \item{\code{m}}{Strong-convexity constant (for logconcave).} |
| 58 | +#' \item{\code{negative_logprob}}{Negative log-probability (for logconcave).} |
| 59 | +#' \item{\code{negative_logprob_gradient}}{Negative log-probability gradient (for logconcave).} |
| 60 | +#' } |
| 61 | +#' @param seed Optional. A fixed seed for the number generator. |
| 62 | +#' @param ... Additional arguments passed along to methods. |
| 63 | +#' |
| 64 | +#' @references Robert L. Smith, "Efficient Monte Carlo Procedures for Generating Points |
| 65 | +#' Uniformly Distributed Over Bounded Regions," Operations Research, 1984. |
| 66 | +#' @references B.T. Polyak, E.N. Gryazina, "Billiard walk - a new sampling algorithm for control |
| 67 | +#' and optimization," IFAC Proceedings Volumes, 2014. |
| 68 | +#' @references Y. Chen, R. Dwivedi, M. J. Wainwright and B. Yu, "Fast MCMC Sampling Algorithms on |
| 69 | +#' Polytopes," Journal of Machine Learning Research, 2018. |
| 70 | +#' @references Lee, Yin Tat, Ruoqi Shen, and Kevin Tian, "Logsmooth Gradient Concentration and |
| 71 | +#' Tighter Runtimes for Metropolized Hamiltonian Monte Carlo," arXiv:2002.04121, 2020. |
| 72 | +#' @references Shen, Ruoqi, and Yin Tat Lee, "The randomized midpoint method for log-concave |
| 73 | +#' sampling," Advances in Neural Information Processing Systems, 2019. |
| 74 | +#' @references Augustin Chevallier, Sylvain Pion, Frederic Cazals, "Hamiltonian Monte Carlo with |
| 75 | +#' boundary reflections, and application to polytope volume calculations," Research Report |
| 76 | +#' preprint hal-01919855, 2018. |
| 77 | +#' |
| 78 | +#' @return A \eqn{d \times n} matrix that contains, column-wise, the sampled points from the |
| 79 | +#' convex polytope P. |
| 80 | +#' @examples |
| 81 | +#' # uniform distribution from the 3d unit cube in H-representation using ball walk |
| 82 | +#' P = gen_cube(3, 'H') |
| 83 | +#' points = sample_points(P, n = 100, random_walk = list("walk" = "BaW", "walk_length" = 5)) |
| 84 | +#' |
| 85 | +#' # gaussian distribution from the 2d unit simplex in H-representation with variance = 2 |
| 86 | +#' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) |
| 87 | +#' b = c(0,0,1) |
| 88 | +#' P = Hpolytope(A=A,b=b) |
| 89 | +#' points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2)) |
| 90 | +#' |
| 91 | +#' # uniform points from the boundary of a 2-dimensional random H-polytope |
| 92 | +#' P = gen_rand_hpoly(2,20) |
| 93 | +#' points = sample_points(P, n = 100, random_walk = list("walk" = "BRDHR")) |
| 94 | +#' |
| 95 | +#' # For sampling from logconcave densities see the examples directory |
| 96 | +#' |
| 97 | +#' @export |
| 98 | +sample_points <- function(P, n, random_walk = NULL, distribution = NULL, seed = NULL, ...) { |
| 99 | + UseMethod("sample_points") |
| 100 | +} |
| 101 | + |
| 102 | +#' @rdname sample_points |
| 103 | +#' @export |
| 104 | +sample_points.default <- function(P, n, random_walk = NULL, distribution = NULL, seed = NULL, ...) { |
| 105 | + sample_points_internal(P, n, random_walk, distribution, seed) |
| 106 | +} |
| 107 | + |
| 108 | +#' @rdname sample_points |
| 109 | +#' @export |
| 110 | +sample_points.Spectrahedron <- function(P, n, random_walk = NULL, distribution = NULL, seed = NULL, ...) { |
| 111 | + stop("Sampling from Spectrahedron is not implemented yet in Rvolesti.") |
| 112 | +} |
0 commit comments