  We demonstrate for the first time that ill-conditioned, non-smooth,
constrained distributions in very high dimension, upwards of 100,000, can be
sampled efficiently $\textit{in practice}$. Our algorithm incorporates
constraints into the Riemannian version of Hamiltonian Monte Carlo and
maintains sparsity. This allows us to achieve a mixing rate independent of
smoothness and condition numbers.
  On benchmark data sets in systems biology and linear programming, our
algorithm outperforms existing packages by orders of magnitude. In particular,
we achieve a 1,000-fold speed-up for sampling from the largest published human
metabolic network (RECON3D). Our package has been incorporated into the COBRA
toolbox.
