Tau-leaped parallel discrete Hamiltonian sampler

Use this as an inference-time sampler for categorical diffusion, discrete energy-based models, or VAE latent refinement where exact Gillespie simulation is too serial. At each sampler step, evaluate U on all one-coordinate neighbors in a tensorized batch. For a configuration with d coordinates and k possible replacements, create a [batch,d,k] neighbor tensor, obtain delta_e=U(x^e)-U(x), and compute q_e=relu(p_e)*exp(-0.5*delta_e), using the persistent-momentum rate from the formula. Choose h so that the mean total event count h*sum_e q_e is between 0.1 and 1.0. Generate Bernoulli indicators B_e with probability 1-exp(-h*q_e). Resolve collisions independently per coordinate: among active replacements for coordinate i, choose one using normalized q_e, and leave the coordinate unchanged if none is active. Update momenta with p<-rho*p+sqrt(1-rho^2)*xi after each leap, or after every K leaps. The paper’s mathematical approximation is the Poisson/Bernoulli count; estimate discretization bias empirically by comparing to exact Gillespie trajectories at matched physical time. First benchmark a 20-dimensional, 32-category synthetic energy model and a small discrete diffusion model on binarized MNIST. Compare wall-clock effective samples and reverse-sampling quality against sequential Gibbs and exact Gillespie. Sweep h and report target-KL, mode coverage, GPU utilization, and samples per second. Success is at least 4x wall-clock throughput at a tau-leap setting whose KL or generated-sample metrics remain within 1% of exact sampling.