Novel matrix hit and run for sampling polytopes and its GPU implementation

We propose and analyze a new Markov Chain Monte Carlo algorithm that generates a uniform sample over full and non-full-dimensional polytopes. This algorithm, termed “Matrix Hit and Run” (MHAR), is a modification of the Hit and Run framework. For a polytope in $$\mathbb {R}^n$$ R n defined by m linear constraints, the regime $$n^{1+\frac{1}{3}} \ll m$$ n 1 + 1 3 ≪ m has a lower asymptotic cost per sample in terms of soft-O notation ( $$\mathcal {O}^*$$ O ∗ ) than do existing sampling algorithms after a warm start. MHAR is designed to take advantage of matrix multiplication routines that require less computational and memory resources. Our tests show this implementation to be substantially faster than the hitandrun R package, especially for higher dimensions. Finally, we provide a python library based on PyTorch and a Colab notebook with the implementation ready for deployment in architectures with GPU or just CPU.