Uniformly sampling high-dimensional polytopes

The function rlim() provides access to implementations of two MCMC methods for uniformly sampling high-dimensional polytopes: the so-called mirror and billard walks (MiW and BiW).

rlim(
  lim,
  Hpol = NULL,
  type = "MiW",
  jmp = NULL,
  scale = 10,
  thin = 1,
  burn = 0,
  nsamp = 3000,
  nsim = NULL,
  starting_point = NULL,
  tol = sqrt(.Machine$double.eps),
  seed = NULL
)

Arguments

lim: A list with four components A, B, G and H representing the polytope to be sampled; see the section Details below.
Hpol: A Hpolytope object representing the reduced polytope to be sampled. If NULL (the default), it is computed from lim. Hpol is used only if lim is set to NULL.
type: A character string specifying the MCMC algorithm to be used for sampling,
jmp: A numeric value (if type is "BiW") or a numeric vector (if type is "MiW") giving the jump rate parameter of the MCMC algorithm to be used. If NULL (the default), it is automatically computed from lim; see the section Details below.
scale: A numeric value. See the section Details below for details.
thin: An integer giving the thinning parameter for the MCMC algorithm. Default is 1 (no thinning).
burn: An integer giving the burning period for the MCMC algorithm. Default is 0 (no burn-in).
nsamp: An integer; the length of the output sample of points into the polytope. Default is 3000.
nsim: An integer value giving the number of sampled points in the polytope. If not NULL (the default value), then nsamp is set to (nsim-burn)/thin. See Details Section below.
starting_point: A numeric vector giving the coordinates of a point inside the polytope, used as starting point in the MCMC algorithm.
tol: A numeric value specifying the tolerance for numeric computations.
seed: An integer used to set the seed of the PRNG.

Value

A nsamp*p matrix whose rows are the coordinates of the nsamp points sampled in the polytope. The dimension p is:

the dimension of the inner space the polytope relies into (i.e., the number of flows), if Hpol is NULL;
the dimension of Hpol is Hpol is an object of class Hpolytope.

Details

A polytope \(\mathcal{P}\) is a convex subset or \(\mathbb{R}^n\), \(n \in \mathbb{N^*}\), defined as the intersection of hyper-planes and half-spaces. More precisely, \(\mathcal{P} = \{ x \in \mathbb{R}^n: Ax = B, Gx \geq H \}\), where \(A\) is an \(m\times n\) matrix, with \(m \leq n\), \(B \in \mathbb{R}^m\), \(G\) is a \(k \times n\) matrix and \(H \in \mathbb{R}^k\). For polytopes involved in Linear Inverse Models (LIM), \(n\) is the number of flows in the metabolic network. The \(k\) inequality constraints insures that \(\mathcal{P}\) is bounded. The polytope to be sampled is given to rlim() whether by:

creating a list with four elements A, B, G, H, giving the set of equality and inequality constraints for the polytope, and using it for lim argument. This list can be created manually or by importing a LIM declaration file thanks to the function df2lim(). Or by
creating an object of class Hypolytope (originating from package volesti) thanks to the function Hpolytope(); see Hpolytope-class. This Hypolytope object is passed to argument Hpol, and requires argument lim to be set to NULL for being processed.

The main practical difference between lim and Hpolytope objects is that Hypolytope allows for inequality constraints only, hence the polytope is of full dimension \(n\).

Two MCMC methods for sampling a polytope are implemented. The MiW has been introduced by Van Oevelen et al. (2010), while the BiW has been introduced by Polyak and Gryazina (2014). Both are reflective Hamiltonian algorithms: each new point of the sample is generated from the previous one by choosing randomly a direction and moving forward this direction for a random length. If the end point of this segment belongs to the polytope, then it is added to the sample; if not, then, the trajectory is reflected on borders of polytope until the end point belongs to the polytope. BiW and MiW differ in the probability distributions used for generating the length of the trajectory; see Girardin at al. (2023) for details. Still, these distributions depend on a parameter called (in both cases) the jump rate of the trajectory.

The jump rate has a strong influence on the performance of the algorithm: too small, a very large number of points nsim would be needed to efficiently sample the polytope. Too large, the trajectories reflect an important number of times before resulting in a point inside the polytope, making the process very slow. If jmp is NULL (default), it is computed automatically, depending the type (MiW or BiW). Precisely,

if type is "MiW", then jmp is a random vector of length equal to the dimension of the reduced polytope. It is equal to the ranges of the reduced polytope divided by scale (default for scale is 10);
if type is "BiW", the jmp is a single value equal to the radius of the largest ball included in the reduced polytope.

References

D. Van Oevelen, K. Van den Meersche, F. J. R. Meysman, K. Soetaert, J. J. Middelburg and A. F. Vézina, Quantifying Food Web Flows Using Linear Inverse Models, Ecosystems 13, 32-45 (2010).B.T. Polyak and E.N. Gryazina, Billiard walk - a new sampling algorithm for control and optimization, IFAC Proceedings Volumes, 47(3), 6123-6128 (2014).V. Girardin, T. Grente, N. Niquil and P. Regnault, Comparing and updating R packages of MCMC Algorithms for Linear Inverse Modeling of Metabolic Networks, hal: (2023)

Examples