polytope — Uniform sampling from convex polytopes

This module provides functions to uniformly sample points subject to a system of linear inequality constraints, \(Ax <= b\) (convex polytope), and linear equality constraints, \(Ax = b\) (affine projection).

A comparison of MCMC algorithms to generate uniform samples over a convex polytope is given in [Chen2018]. Here, we use the Hit & Run algorithm described in [Smith1984].

The R-package hitandrun provides similar functionality to this module.

References

[Chen2018]Chen Y., Dwivedi, R., Wainwright, M., Yu B. (2018) Fast MCMC Sampling Algorithms on Polytopes. JMLR, 19(55):1−86 https://arxiv.org/abs/1710.08165
[Smith1984]Smith, R. (1984). Efficient Monte Carlo Procedures for Generating Points Uniformly Distributed Over Bounded Regions. Operations Research, 32(6), 1296-1308. www.jstor.org/stable/170949
diversipy.polytope.sample(n_points, lower, upper, A1=None, b1=None, A2=None, b2=None, thin=1)

Sample a number of points from a convex polytope A1 x <= b1 using the Hit & Run algorithm.

Lower and upper bounds need to be provided to ensure that the polytope is bounded. Equality constraints A2 x = b2 may be optionally provided.

Parameters:
  • n_points (int) – Number of samples to generate.
  • lower (1d-array of shape (dimension)) – Lower bound in each dimension. If not wanted set to -np.inf.
  • upper (1d-array of shape (dimension)) – Upper bound in each dimension. If not wanted set to np.inf.
  • A1 (2d-array of shape (n_constraints, dimension)) – Left-hand-side of A1 x <= b1.
  • b1 (1d-array of shape (n_constraints)) – Right-hand-side of A1 x <= b1.
  • A2 (2d-array of shape (n_constraints, dimensions), optional) – Left-hand-side of A2 x = b2.
  • b2 (1d-array of shape (n_constraints), optional) – Right-hand-side of A2 x = b2.
  • thin (int, optional) – The thinning factor of the generated samples. A thinning of 10 means a sample is taken every 10 steps.
Returns:

Points sampled from the polytope.
Return type:

2d-array of shape (n_points)
diversipy.polytope.hitandrun(A, b, x0)

Generator for uniform sampling from the convex polytope Ax <= b using the Hit & Run algorithm described in [Smith1984].

Parameters:
  • A (2d-array of shape (n_constraints, dimension)) – Left-hand-side of Ax <= b.
  • b (1d-array of shape (n_constraints)) – Right-hand-side of Ax <= b.
  • x0 (1d-array of shape (dimension)) – Initial point that satisfies A x0 <= b.
Yields:

1d-array of shape (dimension) – Point sampled from the polytope.
diversipy.polytope.chebyshev_center(A, b)

Find the center of the polytope Ax <= b.

Parameters:
  • A (2d-array of shape (n_constraints, dimension)) – Left-hand-side of Ax <= b.
  • b (1d-array of shape (n_constraints)) – Right-hand-side of Ax <= b.
Returns:

Chebyshev center of the polytope
Return type:

1d-array of shape (dimension)
diversipy.polytope.constraints_from_bounds(lower, upper)

Construct the inequality constraints Ax <= b that correspond to the given lower and upper bounds.

Parameters:
  • lower (array-like) – lower bound in each dimension
  • upper (array-like) – upper bound in each dimension
Returns:

  • A (2d-array of shape (2 * dimension, dimension)) – Left-hand-side of Ax <= b.
  • b (1d-array of shape (2 * dimension)) – Right-hand-side of Ax <= b.
diversipy.polytope.solve_equality(A, b)

Compute a basis of the nullspace of A, and a particular solution to Ax = b. This allows to to construct arbitrary solutions as the sum of any vector in the nullspace, plus the particular solution.

Parameters:
  • A (2d-array of shape (n_constraints, dimension)) – Left-hand-side of Ax <= b.
  • b (1d-array of shape (n_constraints)) – Right-hand-side of Ax <= b.
Returns:

  • N (2d-array of shape (dimension, dimension)) – Orthonormal basis of the nullspace of A.
  • xp (1d-array of shape (dimension)) – Particular solution to Ax = b.
diversipy.polytope.check_Ab(A, b)

Check if matrix equation Ax=b is well defined.

Parameters:
  • A (2d-array of shape (n_constraints, dimension)) – Left-hand-side of Ax <= b.
  • b (1d-array of shape (n_constraints)) – Right-hand-side of Ax <= b.