Support sparse sampling with crhmc and edit examples

R/HpolytopeSparseClass.R

@@ -0,0 +1,51 @@
+#' An R class to represent an H-polytope
+#'
+#' An H-polytope is a convex polytope defined by a set of linear inequalities or equivalently a \eqn{d}-dimensional H-polytope with \eqn{m} facets is defined by a \eqn{m\times d} matrix A and a \eqn{m}-dimensional vector b, s.t.: \eqn{Ax\leq b}.
+#'
+#' \describe{
+#'    \item{A}{An \eqn{m\times d} numerical matrix.}
+#'
+#'    \item{b}{An \eqn{m}-dimensional vector b.}
+#'
+#'    \item{volume}{The volume of the polytope if it is known, \eqn{NaN} otherwise by default.}
+#'
+#'    \item{type}{A character with default value 'Hpolytope', to declare the representation of the polytope.}
+#' }
+#'
+#' @examples
+#' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE)
+#' b = c(0,0,1)
+#' P = Hpolytope(A = A, b = b)
+#'
+#' @name HpolytopeSparse-class
+#' @rdname HpolytopeSparse-class
+#' @exportClass HpolytopeSparse
+library(Matrix)
+HpolytopeSparse <- setClass (
+  # Class name
+  "HpolytopeSparse",
+
+  # Defining slot type
+  representation (
+    Aineq = "sparseMatrix",
+    bineq = "numeric",
+    Aeq = "sparseMatrix",
+    beq = "numeric",
+    lb = "numeric",
+    ub = "numeric",
+    dimension = "numeric",
+    type = "character"
+    ),
+
+  # Initializing slots
+  prototype = list(
+    Aineq = as(0, "CsparseMatrix"),
+    bineq = as.numeric(NULL),
+    Aeq = as(0, "CsparseMatrix"),
+    beq = as.numeric(NULL),
+    lb = as.numeric(NULL),
+    ub = as.numeric(NULL),
+    dimension = as.numeric(NaN),
+    type = "HpolytopeSparse"
+  )
+)

R/RcppExports.R

@@ -311,7 +311,17 @@ rounding <- function(P, method = NULL, seed = NULL) {
 #' @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) CRHMC for Riemannian HMC with H-polytope constraints (uniform and general logconcave densities),
+#' xiii) \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 and \code{'CRHMC'}
+#' for logconcave densities with H-polytope and sparse constrainted problems.}
 #' \item{\code{walk_length} }{ The number of the steps per generated point for the random walk. The default value is \eqn{1}.}
 #' \item{\code{nburns} }{ The number of points to burn before start sampling. The default value is \eqn{1}.}
 #' \item{\code{starting_point} }{ A \eqn{d}-dimensional numerical vector that declares a starting point in the interior of the polytope for the random walk. The default choice is the center of the ball as that one computed by the function \code{inner_ball()}.}
@@ -360,7 +370,7 @@ rounding <- function(P, method = NULL, seed = NULL) {
 #' # gaussian distribution from the 2d unit simplex in H-representation with variance = 2
 #' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE)
 #' b = c(0,0,1)
-#' P = Hpolytope(A = A, b = b)
+#' P = Hpolytope(A=A,b=b)
 #' points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2))
 #'
 #' # uniform points from the boundary of a 2-dimensional random H-polytope

man/examples/nuts_rand_poly.R

@@ -25,8 +25,8 @@ f <- function(x) (norm_vec(x)^2 + sum(x))
 # Negative log-probability gradient oracle
 grad_f <- function(x) (2 * x + 1)
 
-dimension <- 50
-facets <- 200
+dimension <- 5
+facets <- 20
 
 # Create domain of truncation
 H <- gen_rand_hpoly(dimension, facets)
@@ -39,11 +39,11 @@ P <- Hpolytope(A = Tr$Mat[1:nrow(Tr$Mat), 2:ncol(Tr$Mat)], b = Tr$Mat[,1])
 x_min = matrix(0, dimension, 1)
 
 # Warm start point from truncated Gaussian
-warm_start <- sample_points(P, n = 1, random_walk = list("nburns" = 5000),
+warm_start <- sample_points(P, n = 1, random_walk = list("nburns" = 500),
                             distribution = list("density" = "gaussian", "variance" = 1/2, "mode" = x_min))
 
 # Sample points
-n_samples <- 20000
+n_samples <- 1000
 
 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))
@@ -53,4 +53,6 @@ hist(samples[1,], probability=TRUE, breaks = 100)
 
 psrfs <- psrf_univariate(samples)
 n_ess <- ess(samples)
-
+cat("ESS=", n_ess, "\n")
+cat("\n")
+cat("PSRF=", psrfs, "\n")

man/examples/simple_hmc_rand_poly.R

@@ -25,8 +25,8 @@ f <- function(x) (norm_vec(x)^2 + sum(x))
 # Negative log-probability gradient oracle
 grad_f <- function(x) (2 * x + 1)
 
-dimension <- 50
-facets <- 200
+dimension <- 5
+facets <- 20
 
 # Create domain of truncation
 H <- gen_rand_hpoly(dimension, facets)

man/examples/sparse_crhmc.R

@@ -31,7 +31,7 @@ A = matrix(c(1, -1), ncol=1, nrow=2, byrow=TRUE)
 b = c(2,1)
 
 # Create domain of truncation
-P <- volesti::Hpolytope$new(A, b)
+P <- volesti::Hpolytope(A=A, b=b)
 
 # Mode of logconcave density
 x_min <- c(-0.5)
@@ -41,7 +41,7 @@ L <- 2
 m <- 2
 
 # Sample points
-n_samples <- 80000
+n_samples <- 1000
 n_burns <- n_samples / 2
 cat("---Sampling without hessian\n")
 pts <- sample_points(P, n = n_samples, random_walk = list("walk" = "CRHMC", "step_size" = 0.3, "nburns" = n_burns, "walk_length" = 1, "solver" = "implicit_midpoint"), distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f, "L_" = L, "m" = m))
@@ -79,16 +79,16 @@ invisible(capture.output(dev.off()))
 
 walk="CRHMC"
 library(Matrix)
-bineq=matrix(c(10,10,10,10,10), nrow=5, ncol=1, byrow=TRUE)
+bineq=c(10,10,10,10,10)
 Aineq = matrix(c(1,0,-0.25,-1,2.5,1,0.4,-1,-0.9,0.5), nrow=5, ncol=2, byrow = TRUE)
 Aineq = as( Aineq, 'dgCMatrix' )
-beq=matrix(,nrow=0, ncol=1, byrow=TRUE)
+beq=(0)
 Aeq = matrix(, nrow=0, ncol=2, byrow = TRUE)
 Aeq=as( Aeq, 'dgCMatrix' )
 lb=-100000*c(1,1);
 ub=100000*c(1,1);
 cat("---Sampling the normal distribution in a pentagon\n")
-P <- volesti::sparse_constraint_problem$new(Aineq, bineq,Aeq, beq, lb, ub)
+P <- volesti::HpolytopeSparse(Aineq=Aineq, bineq=bineq, Aeq=Aeq, beq=beq, lb=lb, ub=ub)
 points <- sample_points(P, n = n_samples, random_walk = list("walk" = "CRHMC", "step_size" = 0.3, "nburns" = n_burns, "walk_length" = 1, "solver" = "implicit_midpoint"), distribution = list("density" = "logconcave", "variance" = 8))
 jpeg("pentagon.jpg")
 plot(ggplot(data.frame( x=points[1,], y=points[2,] )) +