  Hamiltonian Monte Carlo (HMC) is a popular sampling method in Bayesian
inference. Recently, Heng & Jacob (2019) studied Metropolis HMC with couplings
for unbiased Monte Carlo estimation, establishing a generic parallelizable
scheme for HMC. However, in practice a different HMC method, multinomial HMC,
is considered as the go-to method, e.g. as part of the no-U-turn sampler. In
multinomial HMC, proposed states are not limited to end-points as in Metropolis
HMC; instead points along the entire trajectory can be proposed. In this paper,
we establish couplings for multinomial HMC, based on optimal transport for
multinomial sampling in its transition. We prove an upper bound for the meeting
time - the time it takes for the coupled chains to meet - based on the notion
of local contractivity. We evaluate our methods using three targets: 1,000
dimensional Gaussians, logistic regression and log-Gaussian Cox point
processes. Compared to Heng & Jacob (2019), coupled multinomial HMC generally
attains a smaller meeting time, and is more robust to choices of step sizes and
trajectory lengths, which allows re-use of existing adaptation methods for HMC.
These improvements together paves the way for a wider and more practical use of
coupled HMC methods.
