Efficient sampling of constraint spaces in practice

Efficient sampling from constraint spaces is a critical task in various scientific and industrial fields, especially when dealing with high-dimensional, complex distributions. In practice, the goal is to generate samples that respect the underlying constraints of the problem, which could involve specific conditions like boundaries, relationships, or optimization criteria.

1. Efficient sampling of constraint spaces in
practice
Vissarion Fisikopoulos
joint withApostolos Chalkis, Elias Tsigaridas
PGMODAYS 2024

2. (Truncated) distributions
▶Multivariate probability distribution with density functionπ(x)
▶Truncateπto a polytopeP:={Ax≤b}we obtain p.d.f.πP
πP(x) =
f(x)π(x)
R
P
π(x)dx
,f(x) =
æ
1,ifx∈P
0,ifx/∈P
The support is the polytopePandπPis the uniform ditribution overP.
In general the support is a convex bodyK⊂R
n

3. Sampling from (truncated) distributions
Problem
Sample (efficiently?) from a (truncated) distribution with density
πK?
Interesting directions
▶Algorithms and (efficient) implementations
▶Complexity bounds(gaps between theoretical bounds and
practical performance, provide (statistical guarantees))
▶Take advantage of the geometry of supportK⊂R
n
(polytope, non-linear, e.g., spectrahedron, basic semi-algebraic set)
▶Applications
(Volume, integration, Bayesian inference, optimization, ...)

4. Simple cases and simplistic approaches
▶Fundamental shapes (hypercube, hypershpere, simplex) admit
efficient methods
▶Acception/rejection sampling does not scale to high
dimensions

5. Geometric Random Walks
▶AGeometric Random Walkstarts at some interior point and
at each step, chosen according
to some.
▶Goal: analyse geometric random walks as Marcov Chains with
some target distributionB p q
Steps of a ball walk.Uniform target distribution

6. Useful questions and terminology
▶Does the random walk converges asymptotically to
the target distribution? (Correctness
▶How fast does it converge?
(Equivalently) How many steps do we have to perform until
we get a point that isϵ-close to a point draw from the target
distribution? (mixing time)
▶Does the initial point of the walk affects the efficiency?
(warm start)
▶What is the
▶Do we assume anything aboutK? (isotropic position
rounded)

7. Target probability distributions
Definition
Letπ(x)∝e
−f(x)
, wheref:R
d
→Ris a convex function.
π(x) is calledlog-concave (LC) probability density.
▶Letπ(x) be restricted toconvex bodyK⊂R
d
.
▶Special cases: Uniform, Gaussian, Exponential/Boltzmann.

8. Ball walk
Ball Walk(K,p, δ,f): convexK⊂R
d
,p∈P, radiusδ,f:
R
d
→R+
1. xinB(p, δ).
2.returnxwith probability min
æ
1,
f(x)
f(p)
œ
;
returnpwith the remaining probability.B p q
If the density is not restricted inK, then it is theMetropolis-Hastings
algorithm.

9. Hit-and-Run
Hit and Run(K,p,f): convexK⊂R
d
, pointp∈P,f:R
d
→
R+
1. ℓthroughp.
2.returna random point on the chordℓ∩Kchosen from
the distributionπℓ,frestricted inK∩ℓ.` p q ` p q
▶Q: How do we computeℓ∩K? Can we do itexactly?

10. Billiard walk - Uniform case
BW(K,pi, τ,R)
1. L=−τlnη,η∼U(0,1).
2. vto define the trajectory. then the
direction becomesv←v−2⟨v,s⟩.
3.
s,||s||= 1,
4.returnthe end of the trajectory aspi+1.
If the number of reflections exceedsR, thenreturnpi+1=pi.

11. Billiard walk - Uniform case
BW(K,pi, τ,R)
1. L=−τlnη,η∼U(0,1).
2. vto define the trajectory. then the
direction becomesv←v−2⟨v,s⟩.
3.
s,||s||= 1,
4.returnthe end of the trajectory aspi+1.
If the number of reflections exceedsR, thenreturnpi+1=pi.

12. Billiard walk - Uniform case
BW(K,pi, τ,R)
1. L=−τlnη,η∼U(0,1).
2. vto define the trajectory. then the
direction becomesv←v−2⟨v,s⟩.
3.
s,||s||= 1,
4.returnthe end of the trajectory aspi+1.
If the number of reflections exceedsR, thenreturnpi+1=pi.

13. Billiard walk - Uniform case
BW(K,pi, τ,R)
1. L=−τlnη,η∼U(0,1).
2. vto define the trajectory. then the
direction becomesv←v−2⟨v,s⟩.
3.
s,||s||= 1,
4.returnthe end of the trajectory aspi+1.
If the number of reflections exceedsR, thenreturnpi+1=pi.

14. Billiard walk - Uniform case
BW(K,pi, τ,R)
1. L=−τlnη,η∼U(0,1).
2. vto define the trajectory. then the
direction becomesv←v−2⟨v,s⟩.
3.
s,||s||= 1,
4.returnthe end of the trajectory aspi+1.
If the number of reflections exceedsR, thenreturnpi+1=pi.

15. Hamiltonian Monte Carlo
▶Similar to billiard walk but with non-linear trajectory
▶Trajectory is defined by Hamiltonian dynamics simulated using
a time-reversible and volume-preserving numerical integrator
(typically the leapfrog integrator)
▶Reflected: Kby using boundary
reflections.
▶Riemannian: Kthe trajectory is always
insideK.

16. Mixing time experiment (uniform case)
▶Uniform sampling from the hypercube [−1,1]
200
and
projection toR
3
.
▶Rows:,,
Random Directions Hit and Run,.
▶Columns: walk length,{1, 50, 100, 150, 200}

17. Complexity bounds
Year & Authors Random walk Mixing time
∗
Distribution
[Smith: 1986]
e
O(d
3
) any LC
[Berbee, Smith: 1987]
e
O(d
10
) any LC
[Lovasz,Simonovits’90]
e
O(d
3
) any LC
[Kannan,Narayanan’12]
e
O(d
2
) uniform (H-polytope)
[Polyak,Dabbene’14] Billiard walk ?? uniform
[Afshar,Domke’15] Reflective HMC ?? any LC (polytopes)
[Lee,Vempala’16] O(md
3/4
) uniform (H-polytope)
[Lee,Vempala’17]
e
O(md
2/3
) any LC (H-polytopes)
[Chen,Dwivedi,Wainwright,Yu’19]
e
O(d
5/2
) uniform (H-polytope)
[Chen,Dwivedi,Wainwright,Yu’19] O(m
1/2
d
3/2
) uniform (H-polytope)
▶Cost per sample:cost per step×mixing time (#steps).
▶Thecost per stepdepends on the convex body.
▶Hit-and-Run (HR): widely used & well studied.
▶Coordinate Hit-and-Run (CDHR): seems more efficient than HR in
practice.
▶Most existing software uses either CDHR or HR (H-polytopes).

18. Convex bodies
▶H-polytopes: system of linear inequalities
▶V-polytopes: convex hull of point sets
▶Minkowski sums of polytopes
▶Spectrahedra: feasible sets of linear matrix inequalities

19. Geometric and algebraic oracles
▶Membership oracle (Ball walk)
▶Boundary (intersection) oracle (HnR)
▶Reflection oracle (Billiard, ReHMC)
▶Optimization oracle (Minkowski sums, Secondary polytopes)

20. Volume Computation
Computing the exact volume ofP,
▶is #P-hard for all the representations
▶is open if both H- and V- representations available
▶is APX-hard (oracle model)
▶Randomized approximation algorithms
Theorem
[Dyer, Frieze, Kannan’91] Pand any
0≤ϵ, δ≤1, there is a randomized algorithm which computes an
estimateVs.t. with probability 1−δwe have
(1−ϵ)vol(P)≤V≤(1 +ϵ)vol(P), and the number of oracle calls
ispoly(d,1/ϵ,log(1/δ)).

21. MultiphaseMonteCarlo
Let a sequence of functions{f0, . . . ,fm},fi:R
d
→R. Then,
vol(P) =
Z
P
dx=
Z
P
fm(x)dx
R
P
fm−1(x)dx
R
P
fm(x)dx
· · ·
R
P
f0(x)dx
R
P
f1(x)dx
R
P
dx
R
P
f0(x)dx
Then selectfis.t.,
▶The number of phases,m, is as small as possible.
▶Each integral ratio can be efficiently estimated by sampling
fromπ∝firestricted toP(using geometric random walks).
▶There is a closed formula for
R
P
fm(x)dx.
complexity = #phases×#points per phase×cost per point

22. State-of-the-art
Authors-Year
Complexity
fi random walk
(oracle calls)
[Dyer, Frieze, Kannan’91] e
O(d
23
)
Indicator function
grid walk
of a ball
[Kannan, Lovasz, Simonovits’97]e
O(d
5
)
Indicator function
ball walk
of a ball
[Lovasz, Vempala’03]
e
O(d
4
) Exponential hit-and-run
[Cousins, Vempala’15]
e
O(d
3
) Spherical Gaussians ball walk
▶Cannot be implemented as they are due to large constants in
the complexity and pessimistic theoretical bounds.
Practical algorithms:
adjustments (experimental)
▶[Emiris, Fisikopoulos’14]
▶[Cousins, Vempala’16]
▶[Chevallier, Cazals, Fearnhead’22]
▶[Emiris, Chalkis, F:’23]
latter: most efficient, scales to 1000s dims, supports: V-polytopes,
Zonotopes, extends to spectahedra

23. Applications in finance
Portfolio analysis
▶The set of
a simplex.
▶Constraints on investments yield a general polytope.
▶Portfolios with same
trading price series over time) lie on an ellipsoid.
Randomized geometric tools for anomaly detection in stock markets
[Bachelard,Chalkis,F,Tsigaridas’23]

24. Applications in structural biology
[Chalkis,F, Tsigaridas, Zafeiropoulos]
▶Metabolic networks model the reactions of metabolites in an
organim or system.
▶Each reaction has a flow or rate called.
▶The set of states of the network where fluxes are in balance
(rate of production = rate of consumption) is a convex
polytope.
▶Sampling from polytope yield probability densities for reaction
fluxes (example: thioredoxin)

25. GeomScale org
https://geomscale.github.io
C++ volume approximation & sampling from convex
bodies
Python interface with extra tools for metabolic
network analysis (FBA, copulas, visualization)
R interface with extra tools for finance (portfolio
analysis)
————————————————————————————
NumFOCUS Affiliated Project.Support from an open community.More than 15 000 lines of code.
Thank you!